[SDOI2018]旧试题解题报告

解题报告

[SDOI2018]旧试题解题报告：

题意：求$\sum_{i=1}^A\sum_{j=1}^B\sum_{k=1}^C d(i\cdot j\cdot k)\bmod (10^9+7)$其中$d(x)=\sum_{w\mid x}1$即$x$的因数个数

做这道题前，请保证你已经学会莫比乌斯反演并且做出了一部分题目，否则可以到我博珂里康康。

我们先证明一个毒瘤式子$d(x\cdot y\cdot y)=\sum_{i\mid x}\sum_{j\mid y}\sum_{k\mid z}[\gcd(i,j,k)==1]$

当$x=\prod_{i=1}^s p_i^{k_i}$因为$d(x)=\prod_{i=1}^s(k_i+1)$为$x$的因数个数，因此考虑将$x\cdot y\cdot z$的每个因数对应一组互质的因数

对于每个质数$p\mid x\cdot y\cdot z$，且$x$中$p$的幂次为$a$，$y$中$p$的幂次为$b$，$z$中$p$的幂次为$c$，则$x\cdot y\cdot z$中$p$的幂次为$a+b+c$。

根据$d$的积性与$\gcd(p^{a+b+c},\frac{x\cdot y\cdot z}{p^{a+b+c}})==1$，得：$d(x\cdot y\cdot z)=d(\frac{x\cdot y\cdot z}{p^{a+b+c}})\cdot d(p^{a+b+c})=d(\frac{x\cdot y\cdot z}{p^{a+b+c}})\cdot(a+b+c+1)$，即$p$对左式的贡献为$a+b+c+1$

令$i\mid x,j\mid y,k\mid z$且$\gcd(i,j,k)=1$，则$i$中$p$的幂次$u$不超过$a$，$j$中$p$的幂次$v$不超过$b$，$k$中$p$的幂次$w$不超过$c$，且$u$，$v$，$w$中起码有两个为$0$

• 当$v=w=0,u\ne 0$时，$u=1,2,\cdots a$，即有$a$种取值
• 当$u=w=0,v\ne 0$时，$v=1,2,\cdots b$，即有$b$种取值
• 当$u=v=0,w\ne 0$时，$w=1,2,\cdots c$，即有$c$种取值
• 当$u=v=w=0$时，只有$1$种取值

即每个$p$可以对右式同样产生$a+b+c+1$的贡献，因此对于每个$x\cdot y\cdot z$的质因子$p$，其会对左右两式均产生$a+b+c+1$即$p$的指数$+1$的贡献

因此我们可以把原式化为：$\sum_{i=1}\sum_{j=1}^B\sum_{k=1}^C\sum_{x\mid i}\sum_{y\mid j}\sum_{z\mid k}[\gcd(x,y,z)==1]$
$=\sum_{i=1}\sum_{j=1}^B\sum_{k=1}^C\sum_{x\mid i}\sum_{y\mid j}\sum_{z\mid k}[\gcd(x,y)==1][\gcd(y,z)==1][\gcd(x,z)==1]$

改变枚举顺序可得：$\sum_{x=1}^A\lfloor\frac{A}{x}\rfloor\sum_{y=1}^B\lfloor\frac{B}{y}\rfloor\sum_{z=1}^C\lfloor\frac{C}{z}\rfloor[\gcd(x,y)==1][\gcd(y,z)==1][\gcd(x,z)==1]$

使用莫比乌斯反演得原式可得：$\sum_{x=1}^A\lfloor\frac{A}{x}\rfloor\sum_{y=1}^B\lfloor\frac{B}{y}\rfloor\sum_{z=1}^C\lfloor\frac{C}{z}\rfloor\sum_{u\mid\gcd(x,y)}\mu(u)\sum_{v\mid\gcd(y,z)}\mu(v)\sum_{w\mid\gcd(x,z)}\mu(w)$

继续交换枚举顺序得：$\sum_{u=1}^{min(A,B)}\mu(u)\sum_{v=1}^{min(B,C)}\mu(v)\sum_{w=1}^{min(A,C)}\mu(w)\sum_{i=1}^{\lfloor\frac{A}{lcm(u,v)}\rfloor}\lfloor\frac{A}{lcm(u,v)\cdot i}\rfloor\sum_{j=1}^{\lfloor\frac{B}{lcm(v,w)}\rfloor}\lfloor\frac{B}{lcm(v,w)\cdot j}\rfloor\sum_{k=1}^{\lfloor\frac{C}{lcm(u,w)}\rfloor}\lfloor\frac{C}{lcm(u,w)\cdot k}\rfloor$

由引理$\lfloor\frac{a}{b\cdot c}\rfloor=\lfloor\frac{\lfloor\frac{a}{b}\rfloor}{c}\rfloor$得：$\sum_{i=1}^{\lfloor\frac{x}{y}\rfloor}\lfloor\frac{x}{y\cdot i}\rfloor=\sum_{i=1}^{\lfloor\frac{x}{y}\rfloor}\lfloor\frac{\lfloor\frac{x}{y}\rfloor}{i}\rfloor$，原式化为：$\sum_{u=1}^{min(A,B)}\mu(u)\sum_{v=1}^{min(B,C)}\mu(v)\sum_{w=1}^{min(A,C)}\mu(w)\sum_{i=1}^{\lfloor\frac{A}{lcm(u,v)}\rfloor}\lfloor\frac{\lfloor\frac{A}{lcm(u,v)}\rfloor}{i}\rfloor\sum_{j=1}^{\lfloor\frac{B}{lcm(v,w)}\rfloor}\lfloor\frac{\lfloor\frac{B}{lcm(v,w)}\rfloor}{j}\rfloor\sum_{k=1}^{\lfloor\frac{C}{lcm(u,w)}\rfloor}\lfloor\frac{\lfloor\frac{C}{lcm(u,w)}\rfloor}{k}\rfloor$

设$f(x)=\sum_{i=1}^x\lfloor\frac{x}{i}\rfloor$，考虑如何预处理$f(x)$：对于每一个$i\leqslant x$，$x$都会得到$i$的贡献$\lfloor\frac{x}{i}\rfloor$，即在$[1,x]$的范围内$i$的倍数个数。

因此$f(x)=\sigma(x)=\sum_{d\mid x}d$即约数和函数，这个函数通过可以枚举每个数的倍数在$O(n\log n)$的时间内预处理。

原式可化为：$\sum_{u=1}^{min(A,B)}\sum_{v=1}^{min(B,C)}\sum_{w=1}^{min(A,C)}\mu(u)\cdot\mu(v)\cdot\mu(w)\cdot\sigma(\frac{A}{lcm(u,v)})\cdot\sigma(\frac{B}{lcm(v,w)})\cdot\sigma(\frac{C}{lcm(u,w)})$

但是最后的式子还是$O(n^3)$的，就比暴力快一点点。

首先，我们可以把所有$\mu(x)=0$的$x$删去，这种$x$有不少，因此可以快很多。

其次，因为$\sigma(0)=0$，因此所有的$lcm(u,v)>A,lcm(v,w)>B,lcm(u,w)>C$可以剪掉。（其实不剪掉也没有关系，对答案没有影响，因此我们将$A,B,C$统一改为$max(A,B,C)$）

为了达到更低的复杂度，我们可以使用一种毒瘤算法：三元环计数。感性理解一下，枚举三个数$u,v,w$，且满足$lcm(u,v)\leqslant A,lcm(v,w)\leqslant B,lcm(u,w)\leqslant C$就相当于在一张只要任意两个点满足其最小公倍数小于等于$max(A,B,C)$边会有一条边的图中枚举三元环。

三元环计数算法可以在$O(m\sqrt{m})$的时间内求出所有三元环，具体见nekko的博客：不常用的黑科技——「三元环」

但是还是会TLE，因此我们要进行常数优化：

去除所有的自环，枚举三个点相同和两个点相同的情况，并用$6$种排列计算三个点都不相同的情况。

改链式前向星变vector，因为vector的内存访问是连续的，因此可以加快一点速度。

快读，O2优化和火车头优化

这样，我们就可以拿到$100pts$的好成绩！

又臭又长的代码：

#include<stdio.h>
#include<vector>
#include<string.h>
using namespace std;
const long long maxn=1000005,maxm=1000005,mod=1000000007;
long long i,j,k,a,b,c,e,mx,mn,T,cnt,tot,ans;
long long lst[maxn],d[maxn],p[maxn],mu[maxn],ok[maxn],ord[maxn],deg[maxn],f[maxn],from[maxm],to[maxm],lcm[maxm],mrk[maxn];
vector<long long>v[maxn],w[maxn];
inline long long min(long long a,long long b){
    return a<b? a:b;
}
inline long long max(long long a,long long b){
    return a>b? a:b;
}
inline long long gcd(long long a,long long b){
    return b==0? a:gcd(b,a%b);
}
int main(){
    p[1]=mu[1]=1;
    for(i=2;i<maxn;i++){
        if(p[i]==0)
            lst[++cnt]=i,mu[i]=-1;
        for(j=1;j<=cnt;j++){
            if(i*lst[j]>=maxn)
                break;
            p[i*lst[j]]=1;
            if(i%lst[j]==0){
                mu[i*lst[j]]=0;
                break;
            }
            mu[i*lst[j]]=-mu[i];
        }
    }
    for(i=1;i<maxn;i++)
        if(mu[i]!=0)
            ok[++tot]=i,ord[i]=tot;
    for(i=1;i<maxn;i++){
        for(j=1;i*j<maxn;j++)
            d[i*j]++;
        f[i]=(f[i-1]+d[i])%mod;
    }
    scanf("%lld",&T);
    while(T--){
        memset(deg,0,sizeof(deg));
        memset(v,0,sizeof(v));
        memset(w,0,sizeof(w));
        scanf("%lld%lld%lld",&a,&b,&c);
        mn=min(min(a,b),c),mx=max(max(a,b),c);
        e=ans=0;
        for(i=1;i<=tot;i++){
            if(ok[i]>mx)
                break;
            for(j=1;j<=tot;j++){
                if(ok[i]*ok[j]>mx)
                    break;
                if(mu[ok[i]*ok[j]]==0)
                    continue;
                for(k=j+1;k<=tot;k++){
                    if(ok[i]*ok[j]*ok[k]>mx)
                        break;
                    if(mu[ok[i]*ok[k]]==0||gcd(ok[j],ok[k])>1)
                        continue;
                    from[++e]=ord[ok[i]*ok[j]],to[e]=ord[ok[i]*ok[k]],lcm[e]=ok[i]*ok[j]*ok[k];
                    deg[from[e]]++,deg[to[e]]++;
                } 
            }
        }
        for(i=1;i<=tot;i++){
            if(ok[i]>mn)
                break;
            ans+=mu[ok[i]]*mu[ok[i]]*mu[ok[i]]*f[a/ok[i]]*f[b/ok[i]]*f[c/ok[i]];
        }
        for(i=1;i<=e;i++){
            v[from[i]].push_back(to[i]),w[from[i]].push_back(lcm[i]);
            v[to[i]].push_back(from[i]),w[to[i]].push_back(lcm[i]);
        }
        for(i=1;i<=tot;i++){
            if(ok[i]>min(a,b))
                break;
            for(j=0;j<v[i].size();j++){
                long long x=ok[i],y=ok[v[i][j]],z=w[i][j];
                ans=(ans+mu[x]*mu[y]*mu[y]*f[a/z]*f[b/z]*f[c/y]%mod+mod)%mod;
                ans=(ans+mu[x]*mu[x]*mu[y]*f[a/x]*f[b/z]*f[c/z]%mod+mod)%mod;
                ans=(ans+mu[x]*mu[x]*mu[y]*f[a/z]*f[b/x]*f[c/z]%mod+mod)%mod;
            }
        }
        memset(v,0,sizeof(v));
        memset(w,0,sizeof(w));
        for(i=1;i<=e;i++){
            if(deg[from[i]]>=deg[to[i]])
                v[from[i]].push_back(to[i]),w[from[i]].push_back(lcm[i]);
            else v[to[i]].push_back(from[i]),w[to[i]].push_back(lcm[i]);
        }
        for(i=1;i<=tot;i++){
            if(ok[i]>mx)
                break;
            for(j=0;j<v[i].size();j++)
                mrk[v[i][j]]=w[i][j];
            for(j=0;j<v[i].size();j++){
                long long x=v[i][j];
                for(k=0;k<v[x].size();k++){
                    long long y=v[x][k],p=mrk[y],q=w[i][j],r=w[x][k];
                    if(mrk[y]==0)
                        continue;
                    long long st1,st2,st3,st4,st5,st6;
                    st1=f[a/p]*f[b/q]*f[c/r]%mod;
                    st2=f[a/p]*f[b/r]*f[c/q]%mod;
                    st3=f[a/q]*f[b/p]*f[c/r]%mod;
                    st4=f[a/q]*f[b/r]*f[c/p]%mod;
                    st5=f[a/r]*f[b/p]*f[c/q]%mod;
                    st6=f[a/r]*f[b/q]*f[c/p]%mod;
                    ans=(ans+mu[ok[i]]*mu[ok[x]]*mu[ok[y]]*(st1+st2+st3+st4+st5+st6)%mod+mod)%mod;
                }
            }
            for(j=0;j<v[i].size();j++)
                mrk[v[i][j]]=0;
        }
        printf("%lld\n",(ans+mod)%mod);
    }
    return 0;
}