题目链接:http://codeforces.com/problemset/problem/235/E
题意:d(n)是n的约数个数函数,求值:sigma(i=1..a,j=1..b,k=1...c,d(i*j*k)),a,b,c<=2000.
236B作为一个simple的Div2B题,是这一题的弱化版,a,b,c<=100,235E的程序也能过236B.
题解:数论函数相关,就去考虑莫比乌斯辣!
某岛国大佬rng58给出了一个非常妙的式子:
sigma(i=1...a,j=1...b,k=1....c,d(i*j*k))=sigma(gcd(i,j)=gcd(j,k)=gcd(i,k)=1,(a/i)*(b/j)*(c/k)) 证明:http://rng.gozaru.jp/b.pdf
菜鸡博主再大致转述一下讲一下自己的想法:
首先,对于二维的情况,d(a*b)=sigma(i|a,j|b,[gcd(i,j)==1]),
这个结论因为 [SDOI2015]约数个数和 这道题而变得十分有名,
证明:考虑每个质数p对于答案的贡献,
设a=a'*p^k1,b=b'*p^k2,
显然质数p对于答案的贡献是(k1+k2+1)*d(a'*b').
而在右式中,只考虑与p有关的部分,共有(k1+k2+1)对,
而每一对两个数分别乘上d(a'*b')所对应的两个数均符合gcd(a,b)=1的条件。
所以上面这个式子得证。
考虑扩展到三维情况,d(a*b*c)=sigma(i|a,j|b,k|c,[gcd(i,j)=1][gcd(j,k)=1][gcd(i,k)=1]).
这个式子事实上也是正确的(该公式对于任意维数都成立),证明类似上面的证明。
然后进行转化,左式=sigma(i=1..a,j=1...b,k=1...c,p|i,q|j,r|k,[gcd(p,q)==gcd(q,r)==gcd(p,r)==1])
考虑每一个p被(a/p)个数整除,每一个q被(b/q)个数整除,每一个r被(c/r)个数整除,
计算每一组(i,j,k)对于答案的贡献,
得到sigma(i=1..a,j=1...b,k=1...c,(a/i)*(b/j)*(c/k)*[gcd(i,j)=1]*[gcd(i,k)=1]*[gcd(j,k)=1])
于是就可以开始莫比乌斯反演了,
利用[n=1]=miu*1化开一个等式: 原式=sigma(i=1...a,(a/i)*(sigma(d=1...min(b,c),miu(d)*sigma(d|j,gcd(i,j)=1,b/j)*sigma(d|k,gcd(i,k)=1,c/k))))
=sigma(i=1...a,(a/i)*(sigma(d=1...min(b,c),miu(d)*sigma(gcd(i,d*j)=1,b/(d*j))*sigma(gcd(i,d*k)=1,c/(d*k))))
因为gcd(i,d*j)=1,所以首先要保证gcd(i,d)=1,然后再保证gcd(i,j)=1,其中j=1....b/d,对于k的部分可以用相同的办法处理,
枚举i的时间是O(a),枚举d的时间是O(min(b,c)),然后枚举j,k的时间分别是sigma(b/i)和sigma(c/i),总时间复杂度是O(a*b*log(b)).
Code:
#include <bits/stdc++.h>
#define ll long long
#define mod 1073741824ll
using namespace std;
int gcd[2005][2005];
ll p[3],a,b,c;
int miu[2005],pri[505],cnt=0;
bool isp[2005];
inline void gmiu()
{miu[1]=1;
for (int i=2;i<=2000;i++)
{if (!isp[i]) {pri[++cnt]=i;miu[i]=-1;}
for (int j=1;j<=cnt&&pri[j]*i<=2000;j++)
{isp[i*pri[j]]=1;
if (i%pri[j]==0) {miu[i*pri[j]]=0;break;}
miu[i*pri[j]]=miu[i]*(-1);
}
}
}
int main (){
ll i,j,d,k;
gmiu();
scanf ("%lld%lld%lld",&p[0],&p[1],&p[2]);
sort(p,p+3);a=p[0];b=p[1];c=p[2];
for (i=1;i<=2000;i++)
{for (j=1;j<=i;j++)
{if (i%j==0) {gcd[i][j]=j;}
else {gcd[i][j]=gcd[j][i%j];}
gcd[j][i]=gcd[i][j];
}
}
ll ans=0;
for (i=1;i<=a;i++)
{ll s1=a/i,s2=0;
for (d=1;d<=b;d++)
{if (gcd[i][d]!=1) {continue;}
ll a1=0,a2=0;
for (j=1;j<=b/d;j++)
{if (gcd[i][j]==1)
{a1+=b/(d*j);a1%=mod;}
}
for (k=1;k<=c/d;k++)
{if (gcd[i][k]==1)
{a2+=c/(d*k);a2%=mod;}
}
s2+=miu[d]*a1*a2;s2%=mod;
if (s2<0) {s2+=mod;}
}
ans+=s1*s2;ans%=mod;
}
printf ("%lld\n",ans);
return 0;
}