51nod 1040 最大公约数之和

分析:只要求出每个最大公约数出现的次数就可以了,而最大公约数必然是n的因子,考虑n的任意一个因子m,设t满足gcd(t,n)=m,等价于gcd(t/m,n/m)=1,由t<=n的t/m<=t/m,也就是说,这样的t的个数恰好是φ(n/m)个,故结果就是(m|n)∑m*φ(n/m).复杂度为O((logn)^2).

  1 #include<iostream>
  2 #include<cstring>
  3 #include<cstdio>
  4 #include<vector>
  5 using namespace std;
  6 typedef long long ll;
  7 int prinum[10000],len=0;
  8 ll Fai(int n0){
  9     int n=n0;
 10     if(n==1)return 1;
 11     for(int i=0;i<len;i++){
 12         if(n%prinum[i]==0){
 13             int t=1;
 14             while(n%prinum[i]==0){
 15                 n/=prinum[i];
 16                 t*=prinum[i];
 17             }
 18             if(n==1){
 19                 return (ll)n0-t/prinum[i];
 20             }
 21             return Fai(t)*Fai(n);
 22         }
 23     }
 24     return n0-1;
 25 }
 26 struct Pair{
 27     int x,y;
 28 };
 29 void CalPri(){
 30     int num[50000];
 31     for(int i=0;i<50000;i++)num[i]=i;
 32     for(int i=2;i<50000;i++){
 33         if(num[i]==0)continue;
 34         prinum[len++]=i;
 35         for(int j=2*i;j<50000;j+=i)
 36             num[j]=0;
 37     }
 38 }
 39 void decompose(int n,vector<Pair> &v){
 40     v.clear();
 41     Pair p;
 42     for(int i=0;i<len&&prinum[i]<=n;i++){
 43         if(n%prinum[i]==0){
 44             p.x=prinum[i];
 45             p.y=0;
 46             while(n%prinum[i]==0){
 47                 n/=prinum[i];
 48                 p.y++;
 49             }
 50             v.push_back(p);
 51         }
 52     }
 53     if(n!=1){
 54         p.x=n;
 55         p.y=1;
 56         v.push_back(p);
 57     }
 58 }
 59 ll qpow(int a,int n){
 60     ll k=1,y=a,ans=1;
 61     while(k){
 62         if(k&n)ans*=y;
 63         y*=y;
 64         k<<=1;
 65     }
 66     return ans;
 67 }
 68 ll solve(int n){
 69     vector<Pair> v;
 70     decompose(n,v);
 71     int r[100];
 72     ll result=0,m;
 73     memset(r,0,sizeof(r));
 74     int i=v.size()-1;
 75     while(i>=0){
 76         m=1;
 77         for(int j=0;j<v.size();j++){
 78             m*=qpow(v[j].x,r[j]);
 79         }
 80         result+=m*Fai(n/m);
 81         if(r[i]>=v[i].y){
 82             while(i>=0&&r[i]>=v[i].y){
 83                 i--;
 84             }
 85             if(i==-1)break;
 86             r[i]++;
 87             memset(r+i+1,0,sizeof(r)-(i+1)*sizeof(int));
 88             i=v.size()-1;
 89         }
 90         else r[i]++;
 91     }
 92     return result;
 93 }
 94 int main(){
 95     CalPri();
 96     int n;
 97     cin>>n;
 98     if(n==1)cout<<1<<endl;
 99     else cout<<solve(n)<<endl;
100     return 0;
101 }