[SDOI2018]旧试题

[SDOI2018]旧试题

t组询问,询问(sum_{i=1}^Asum_{j=1}^Bsum_{k=1}^Cd(ijk)mod 10^9+7),(1 ≤ T ≤ 10, 1 ≤ A, B, C ≤ 10^5, ,)(1≤∑max(A,B,C)≤2 imes10 ^5)

不难由约数拆分公式类比出

[ans=sum_{i=1}^Asum_{j=1}^Bsum_{k=1}^Csum_{x|i}sum_{y|j}sum_{z|k}epsilon(gcd(x,y))epsilon(gcd(y,z))epsilon(gcd(x,z))= ]

[sum_{x=1}^Asum_{y=1}^Bsum_{z=1}^Cepsilon(gcd(x,y))epsilon(gcd(y,z))epsilon(gcd(x,z))[A/x][B/y][C/z] ]

注意到多个元素的gcd,要保证gcd的相等,而倍数型更适用于只有gcd,约数条件的问题,于是考虑辅助公式,于是我们有

[ans=sum_{x=1}^Asum_{y=1}^Bsum_{z=1}^Csum_{i|x,i|y}mu(i)sum_{j|x,j|z}mu(j)sum_{k|y,k|z}mu(k)[A/x][B/y][C/z]= ]

[sum_{i=1}^{min(A,B)}sum_{j=1}^{min(A,C)}sum_{k=1}^{min(B,C)}mu(i)mu(j)mu(k)sum_{lcm(i,j)|x}sum_{lcm(i,k)|y}sum_{lcm(j,k)|z}[A/x][B/y][C/z] ]

于是我们发现三组数,考虑三元环mobius,注意到一组i,j,k要有效,必然它们的(mu)要有值,且小于A,B,C,而显然后式子可以(O(nlog(n)))维护,即两个数(mu)有值,范围在A,B,C内,而lcm也没有超过A,B,C,我们就将其连边,这样我们就得到一张图,而解就是其中的三元环。

建图根据(lcm(i,j)gcd(i,j)=ij),枚举gcd,再枚举互质的两个数,边判(mu)的是否有值以及是否超范围即可。

而三元环寻找显然是(n^3),考虑优化,这里采取度数大的点与度数小的点连边,这样我们的时间复杂度就只有(O(msqrt{m}))(n为点,m为边)


证明:

考虑度数大于(sqrt{m})的点,这样的点不能超过(sqrt{m})个,这里的枚举不超过(msqrt{m}),而小于的点,点不超过n所以枚举次数应为(nsqrt{m}),因此当边数较多时可近似认为时间复杂度(O(msqrt{m}))


于是我们有预处理后式的方法,又能较快建出了图,又有很快地找出三元环的方式,于是问题,注意判自环即可,至此得到很好解决。

参考代码:

#include <iostream>
#include <cstdio>
#include <cstring>
#define il inline
#define ri register
#define ll long long
#define plim 100000
#define yyb 1000000007
#define swap(x,y) x^=y^=x^=y
using namespace std;
struct point{
    int next,to,len;
}ar[plim<<3];int at;
int head[plim+1];
struct edge{
    int p1,p2,len;
}e[plim<<3];int et;
bool check[plim+1];
ll fa[plim+1],fb[plim+1],fc[plim+1];
int pt,prime[10000],mu[plim+1],ds[plim+1],
         opt[plim+1],mark[plim+1];
il int max(int,int),gcd(int,int);
il void prepare(int),link(int,int,int);
int main(){
    ll ans;
    int lsy,A,B,C,mx,i,j,k;
    scanf("%d",&lsy),prepare(plim);
    while(lsy--){
        scanf("%d%d%d",&A,&B,&C),mx=max(A,max(B,C)),ans&=0;
        for(i=1;i<=A;++i)for(j=i;j<=A;j+=i)fa[i]+=A/j;
        for(i=1;i<=B;++i)for(j=i;j<=B;j+=i)fb[i]+=B/j;
        for(i=1;i<=C;++i)for(j=i;j<=C;j+=i)fc[i]+=C/j;
        for(i=1;i<=mx;++i)ans+=mu[i]*fa[i]*fb[i]*fc[i];
        for(i=1;i<=mx;++i){
            if(!mu[i])continue;
            for(j=1;j*i<=mx;++j){
                if(!mu[j*i])continue;
                for(k=j+1;(ll)j*i*k<=mx;++k){
                    if(!mu[i*k]||gcd(j,k)>1)continue;
                    int u(i*j),v(i*k),l(i*j*k);
                    ans+=mu[u]*(fa[v]*fb[l]*fc[l]+fa[l]*fb[v]*fc[l]+fa[l]*fb[l]*fc[v]);
                    ans+=mu[v]*(fa[u]*fb[l]*fc[l]+fa[l]*fb[u]*fc[l]+fa[l]*fb[l]*fc[u]);
                    e[++et].p1=u,e[et].p2=v,e[et].len=l,++ds[u],++ds[v];
                }
            }
        }
        for(i=1;i<=et;++i)
            if(ds[e[i].p1]>ds[e[i].p2])link(e[i].p1,e[i].p2,e[i].len);
            else link(e[i].p2,e[i].p1,e[i].len);
        for(i=1;i<=mx;++i){
            for(j=head[i];j;j=ar[j].next)mark[ar[j].to]=i,opt[ar[j].to]=ar[j].len;
            for(j=head[i];j;j=ar[j].next)
                for(k=head[ar[j].to];k;k=ar[k].next){
                    if(mark[ar[k].to]!=i)continue;
                    int d1(ar[j].len),d2(ar[k].len),d3(opt[ar[k].to]);
                    ans+=mu[i]*mu[ar[j].to]*mu[ar[k].to]*
                        (fa[d1]*fb[d2]*fc[d3]+fa[d1]*fb[d3]*fc[d2]+
                         fa[d2]*fb[d1]*fc[d3]+fa[d2]*fb[d3]*fc[d1]+
                         fa[d3]*fb[d1]*fc[d2]+fa[d3]*fb[d2]*fc[d1]);
                }
        }printf("%lld
",ans%yyb);
        memset(mark,0,sizeof(mark)),memset(fa,0,sizeof(fa)),memset(fb,0,sizeof(fb));
        memset(fc,0,sizeof(fc)),memset(head,0,sizeof(head)),memset(ds,0,sizeof(ds)),
            et&=at&=0;
    }
    return 0;
}
il void link(int x,int y,int len){
    ar[++at].len=len,ar[at].to=y;
    ar[at].next=head[x],head[x]=at;
}
il int gcd(int a,int b){
    while(b)swap(a,b),b%=a;return a;
}
il int max(int a,int b){
    return a>b?a:b;
}
il void prepare(int n){
    int i,j;mu[1]|=true;
    for(i=2;i<=n;++i){
        if(!check[i])prime[++pt]=i,mu[i]=-1;
        for(j=1;j<=pt&&prime[j]*i<=n;++j){
            check[i*prime[j]]|=true;
            if(!(i%prime[j]))break;
            mu[i*prime[j]]=~mu[i]+1;
        }
    }
}

原文地址:https://www.cnblogs.com/a1b3c7d9/p/10824961.html