[Codeforces 235E/236B] Number Challenge 莫比乌斯反演

题目链接：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;
}