题意
- 给定(a,b,c)，求(sum_{i=1}^asum_{j=1}^bsum_{k=1}^c d(i imes j imes k))
- (1leqslant a,b,cleqslant 2000)
分析
第一眼就看出来好像见过这道题，然后想起来它是[P4619 SDOI2018]旧试题的弱化版到底是CCF抄CF的还是CF抄CCF的呢：
我们先证明一个毒瘤式子(d(xcdot ycdot y)=sum_{imid x}sum_{jmid y}sum_{kmid 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)的因数个数，因此考虑将(xcdot ycdot z)的每个因数对应一组互质的因数
对于每个质数(pmid xcdot ycdot z)，且(x)中(p)的幂次为(a)，(y)中(p)的幂次为(b)，(z)中(p)的幂次为(c)，则(xcdot ycdot z)中(p)的幂次为(a+b+c)。
根据(d)的积性与(gcd(p^{a+b+c},frac{xcdot ycdot z}{p^{a+b+c}})==1)，得：(d(xcdot ycdot z)=d(frac{xcdot ycdot z}{p^{a+b+c}})cdot d(p^{a+b+c})=d(frac{xcdot ycdot z}{p^{a+b+c}})cdot(a+b+c+1))，即(p)对左式的贡献为(a+b+c+1)
令(imid x,jmid y,kmid 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 e 0)时，(u=1,2,cdots a)，即有(a)种取值
- 当(u=w=0,v e 0)时，(v=1,2,cdots b)，即有(b)种取值
- 当(u=v=0,w e 0)时，(w=1,2,cdots c)，即有(c)种取值
- 当(u=v=w=0)时，只有(1)种取值
即每个(p)可以对右式同样产生(a+b+c+1)的贡献，因此对于每个(xcdot ycdot z)的质因子(p)，其会对左右两式均产生(a+b+c+1)即(p)的指数(+1)的贡献
因此我们可以把原式化为：(sum_{i=1}sum_{j=1}^Bsum_{k=1}^Csum_{xmid i}sum_{ymid j}sum_{zmid k}[gcd(x,y,z)==1])(=sum_{i=1}sum_{j=1}^Bsum_{k=1}^Csum_{xmid i}sum_{ymid j}sum_{zmid k}[gcd(x,y)==1][gcd(y,z)==1][gcd(x,z)==1])
改变枚举顺序可得：
(sum_{x=1}^Alfloorfrac{A}{x} floorsum_{y=1}^Blfloorfrac{B}{y} floorsum_{z=1}^Clfloorfrac{C}{z} floor[gcd(x,y)==1][gcd(y,z)==1][gcd(x,z)==1])
使用莫比乌斯反演得原式可得：(sum_{x=1}^Alfloorfrac{A}{x} floorsum_{y=1}^Blfloorfrac{B}{y} floorsum_{z=1}^Clfloorfrac{C}{z} floorsum_{umidgcd(x,y)}mu(u)sum_{vmidgcd(y,z)}mu(v)sum_{wmidgcd(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}^{lfloorfrac{A}{lcm(u,v)} floor}lfloorfrac{A}{lcm(u,v)cdot i} floorsum_{j=1}^{lfloorfrac{B}{lcm(v,w)} floor}lfloorfrac{B}{lcm(v,w)cdot j} floorsum_{k=1}^{lfloorfrac{C}{lcm(u,w)} floor}lfloorfrac{C}{lcm(u,w)cdot k} floor)
由引理(lfloorfrac{a}{bcdot c} floor=lfloorfrac{lfloorfrac{a}{b} floor}{c} floor)得：(sum_{i=1}^{lfloorfrac{x}{y} floor}lfloorfrac{x}{ycdot i} floor=sum_{i=1}^{lfloorfrac{x}{y} floor}lfloorfrac{lfloorfrac{x}{y} floor}{i} floor)，原式化为：(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}^{lfloorfrac{A}{lcm(u,v)} floor}lfloorfrac{lfloorfrac{A}{lcm(u,v)} floor}{i} floorsum_{j=1}^{lfloorfrac{B}{lcm(v,w)} floor}lfloorfrac{lfloorfrac{B}{lcm(v,w)} floor}{j} floorsum_{k=1}^{lfloorfrac{C}{lcm(u,w)} floor}lfloorfrac{lfloorfrac{C}{lcm(u,w)} floor}{k} floor)
设(f(x)=sum_{i=1}^xlfloorfrac{x}{i} floor)，考虑如何预处理(f(x))：对于每一个(ileqslant x)，(x)都会得到(i)的贡献(lfloorfrac{x}{i} floor)，即在([1,x])的范围内(i)的倍数个数。
因此(f(x)=sigma(x)=sum_{dmid x}d)即约数和函数，这个函数通过可以枚举每个数的倍数在(O(nlog n))的时间内预处理。
原式可化为：(sum_{u=1}^{min(A,B)}sum_{v=1}^{min(B,C)}sum_{w=1}^{min(A,C)}mu(u)cdotmu(v)cdotmu(w)cdotsigma(frac{A}{lcm(u,v)})cdotsigma(frac{B}{lcm(v,w)})cdotsigma(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(msqrt{m}))的时间内求出所有三元环，具体见nekko的博客：不常用的黑科技——「三元环」
时间复杂度：(O(msqrt{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
",(ans+mod)%mod);
return 0;
}