CF235E 题解(转)
题意
- 给定\(a,b,c\)，求\(\sum_{i=1}^a\sum_{j=1}^b\sum_{k=1}^c d(i\times j\times k)\)
- \(1\leqslant a,b,c\leqslant 2000\)
分析
第一眼就看出来好像见过这道题，然后想起来它是[P4619 SDOI2018]旧试题的弱化版到底是CCF抄CF的还是CF抄CCF的呢：
我们先证明一个毒瘤式子\(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的博客：不常用的黑科技——「三元环」
时间复杂度：\(O(m\sqrt{m})\)。
代码
直接拿的旧试题的代码，加了亿点常数优化，凑合着看吧。
#include<stdio.h>
#include<vector>
#include<string.h>
#define int long long
using namespace std;
const int maxn=1000005,maxm=1000005,mod=1073741824;
int i,j,k,a,b,c,e,mx,mn,T,cnt,tot,ans;
int 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<int>v[maxn],w[maxn];
inline int min(int a,int b){
return a<b? a:b;
}
inline int max(int a,int b){
return a>b? a:b;
}
inline int gcd(int a,int b){
return b==0? a:gcd(b,a%b);
}
signed main(){
p[1]=mu[1]=1;
for(int i=2;i<maxn;i++){
if(p[i]==0)
lst[++cnt]=i,mu[i]=-1;
for(int 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(int i=1;i<maxn;i++)
if(mu[i]!=0)
ok[++tot]=i,ord[i]=tot;
for(int i=1;i<maxn;i++){
for(int j=1;i*j<maxn;j++)
d[i*j]++;
f[i]=(f[i-1]+d[i])%mod;
}
scanf("%lld%lld%lld",&a,&b,&c);
mn=min(min(a,b),c),mx=max(max(a,b),c);
e=ans=0;
for(int i=1;i<=tot;i++){
if(ok[i]>mx)
break;
for(int j=1;j<=tot;j++){
if(ok[i]*ok[j]>mx)
break;
if(mu[ok[i]*ok[j]]==0)
continue;
for(int 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(int 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(int 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(int i=1;i<=tot;i++){
if(ok[i]>min(a,b))
break;
for(int j=0;j<v[i].size();j++){
int 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(int 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(int i=1;i<=tot;i++){
if(ok[i]>mx)
break;
for(int j=0;j<v[i].size();j++)
mrk[v[i][j]]=w[i][j];
for(int j=0;j<v[i].size();j++){
int x=v[i][j];
for(int k=0;k<v[x].size();k++){
int y=v[x][k],p=mrk[y],q=w[i][j],r=w[x][k];
if(mrk[y]==0)
continue;
int 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(int j=0;j<v[i].size();j++)
mrk[v[i][j]]=0;
}
printf("%lld\n",(ans+mod)%mod);
return 0;
}
