CodeForces 235 E.Number Challenge(モビウス反転+数論)

Description
Σi=1 aΣj=1 bΣk=1 cd(ijk)、a,b,c≦2000を求めます
Input
3つの整数a,b,c(1≦a,b,c≦2000)
Output
出力結果モード1073741824
Sample Input
2 2 2
Sample Output
20
Solution
まず2つの結論を証明します.
1 . d(mn)=∑i|m∑j|n[(i,j)=1]
m=pa 11 pa 22とする...paxx,n=pb11pb22...pbxx
mnの任一正因子lに対して、l=pc 11 pc 22と仮定する...pcxx、0≦ci≦ai+bi
s=pd 11 pd 22...pdxx , t=pe11pe22...pexx、di={ci 0 ci≦aici>ai、ei={0 ci−aici≦aici>ai
(s,t)=1,s|m,t|nがあり,この分解が一意であることは明らかである.
逆に、s,tが(s,t)=1,s|m,t|nを満たす場合、ci=pcxx|mn
明らかに唯一を表すべきで,結論をまとめると成立する.
2 . d(rst)=∑i|r∑j|s∑k|t[(i,j)=1][(i,k)=1][(j,k)=1]
d(rst)=∑i|r∑j|st[(i,j)=1]=∑i|s∑j|s∑k|t[(i,jk)=1][(j,k)=1]=∑i|r∑j|s∑k|t[(i,j)=1][(i,k)=1][(j,k)=1]
さらに
∑i=1a∑j=1b∑k=1cd(ijk) ====∑r=1a∑s=1b∑t=1c∑i|r∑j|s∑k|t[(i,j)=1][(i,k)=1][(j,k)=1]∑i=1a∑j=1b∑k=1c⌊ai⌋⌊bj⌋⌊ck⌋[(i,j)=1][(i,k)=1][(j,k)=1]∑i=1a∑j=1b∑k=1c⌊ai⌋⌊bj⌋⌊ck⌋[(i,j)=1][(i,k)=1]∑d|(j,k)μ(d)∑i=1a⌊ai⌋∑d=1min(b,c)μ(d)∑j=1⌊bd⌋⌊bjd⌋[(i,jd)=1]∑k=1⌊cd⌋⌊ckd⌋[(i,kd)=1]
プリプロセッシング
gcd(i,j)と
mu(d)、列挙
a,d,時間複雑度
O(n2log n)
Code

#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
using namespace std;
typedef long long ll;
typedef pair<int,int>P;
const int INF=0x3f3f3f3f,maxn=2005;
#define mod 1073741824 
int gcd[maxn][maxn],mu[maxn],p[maxn],vis[maxn],res;
void init(int n=2000)
{
    res=0;
    mu[1]=1;
    for(int i=2;i<=n;i++)
    {
        if(!vis[i])p[res++]=i,mu[i]=-1;
        for(int j=0;j1;
            if(i%p[j])mu[i*p[j]]=-mu[i];
            else
            {
                mu[i*p[j]]=0;
                break;
            }
        }
    }
    for(int i=0;i<=n;i++)gcd[i][0]=gcd[0][i]=i;
    for(int i=1;i<=n;i++)
        for(int j=1;j<=i;j++)
            gcd[j][i]=gcd[i][j]=gcd[j][i%j];
}
int deal(int n,int d,int i)
{
    int ans=0;
    for(int j=1;j<=n;j++)
        if(gcd[i][j*d]==1)ans+=n/j;
    return ans;
}
int main()
{
    init();
    int a,b,c;
    while(~scanf("%d%d%d",&a,&b,&c))
    {
        int ans=0;
        for(int i=1;i<=a;i++)
        {
            int res=0;
            for(int d=1;d<=min(b,c);d++) 
            {
                res+=((ll)mu[d]*deal(b/d,d,i)*deal(c/d,d,i)%mod+mod)%mod;
                if(res>=mod)res-=mod;
            }
            ans+=(ll)(a/i)*res%mod;
            if(ans>=mod)ans-=mod;
        }
        printf("%d
",ans);
    }
    return 0;
}