spoj 3871 gcd extreme

题目大意给出一个n，求sum(gcd(i,j),0<i<j<=n);
可以明显的看出来s[n]=s[n-1]+f[n];
f[n]=sum(gcd(i,n),0<i<n);
现在麻烦的是求f[n]
gcd(x,n)的值都是n的约数，则f[n]=
sum{i*g(n,i),i是n的约数}，注意到gcd(x,n)=i的
充要条件是gcd(x/i,n/i)=1,因此满足条件的
x/i有phi(n/i)个，说明gcd(n,i)=phi(n/i).
f[n]=sum{i*phi(n/i),1=<i<n}
因此可以搞个for循环对i循环，只要i<n,f[n]+=i*phi(n/i);


#include <iostream>
#include <cstring>
using namespace std;
#define Max 1000000

long long phi[Max+5],ans[Max+5];
int prime[Max/3];
bool flag[Max+5];

void init()
{
    long long  i,j,num=0;
    memset(flag,1,sizeof(flag));
    phi[1]=0;
    for(i=2;i<=Max;i++)//欧拉筛选
    {
        if(flag[i])
        {
            prime[num++]=i;
            phi[i]=i-1;
        }
        for(j=0;j<num && prime[j]*i<=Max;j++)
        {
            flag[i*prime[j]]=false;
            if(i%prime[j]==0)
            {
                phi[i*prime[j]]=phi[i]*prime[j];
                break;
            }
            else phi[i*prime[j]]=phi[i]*(prime[j]-1);
        }
    }
    //for(i=1;i<=10;i++)
    //    cout<<phi[i]<<endl;

    ans[0] =0;
    for(i=1;i*i<=Max;i++){
        ans[i*i] += i*phi[i];
        for(j =i+1;j*i<=Max;j++)
            ans[i*j] += i*phi[j]+j*phi[i];
    }
    for(int i=1;i<=Max;i++)
        ans[i] += ans[i-1];
}

int main(){
    init();
    long long n;
    while(cin>>n&&n){

        cout<<ans[n]<<endl;
    }
}