数学计数原理（Pólya，高精度）：SGU 294 He's Circles

  1 #include <iostream>
  2 #include <cstring>
  3 #include <cstdio>
  4 using namespace std;
  5 const int maxn=200010;
  6 int n,tot;
  7 int phi[maxn],pri[maxn];
  8 
  9 void Linear_Shaker(){
 10     phi[1]=1;
 11     for(int i=2;i<=n;i++){
 12         if(!phi[i]){
 13             phi[i]=i-1;
 14             pri[++tot]=i;
 15         }
 16         for(int j=1;j<=n;j++){
 17             if(i*pri[j]>n)break;
 18             if(i%pri[j]!=0)
 19                 phi[i*pri[j]]=phi[i]*(pri[j]-1);
 20             else{
 21                 phi[i*pri[j]]=phi[i]*pri[j];
 22                 break;
 23             }    
 24         }
 25     }
 26 }
 27 
 28 int POW[4]={1,10,100,1000};
 29 const int mod=10000;
 30 struct ExtInt{
 31     int num[20010],len;
 32     ExtInt(int x){len=0;
 33         memset(num,0,sizeof(num));
 34         do{
 35             num[++len]=x%mod;
 36             x/=mod;
 37         }while(x);
 38     }
 39     
 40     void Scf(){
 41         char s[10010];scanf("%s",s+1);
 42         memset(num,0,sizeof(num));len=1;
 43         for(int i=strlen(s+1),cnt=0;i>=1;i--){
 44             num[len]+=POW[cnt++]*(s[i]-‘0‘);
 45             if(cnt==4)cnt=0,len+=1;
 46         }
 47     }
 48     
 49     int operator [](int x){
 50         return num[x];
 51     }    
 52     
 53     void Prf(){
 54         printf("%d",num[len]);
 55         for(int i=len-1;i>=1;i--)
 56             printf("%04d",num[i]);
 57         printf("\n");    
 58     }
 59 };
 60 
 61 ExtInt operator +(ExtInt a,int b){
 62     ExtInt ret(0);
 63     ret.len=a.len;
 64     for(int i=1,in=0;i<=ret.len||in;i++){
 65         ret.num[i]=a[i]+b+in;in=ret[i]/mod;
 66         ret.num[i]%=mod;ret.len=max(ret.len,i);
 67     }
 68     return ret;
 69 }
 70 
 71 ExtInt operator +(ExtInt a,ExtInt b){
 72     ExtInt ret(0);
 73     ret.len=max(a.len,b.len);
 74     for(int i=1,in=0;i<=ret.len||in;i++){
 75         ret.num[i]=a[i]+b[i]+in;in=ret[i]/mod;
 76         ret.num[i]%=mod;ret.len=max(ret.len,i);
 77     }
 78     return ret;
 79 }
 80 
 81 ExtInt operator *(ExtInt a,ExtInt b){
 82     ExtInt ret(0);
 83     for(int i=1;i<=a.len;i++){
 84         for(int j=1,in=0;j<=b.len||in;j++){    
 85             ret.num[i+j-1]+=a[i]*b[j]+in;in=ret[i+j-1]/mod;
 86             ret.num[i+j-1]%=mod;ret.len=max(ret.len,i+j-1);
 87         }
 88         while(!ret[ret.len])
 89             ret.num[ret.len--]=0;
 90         ret.len=max(ret.len,1);    
 91     }
 92     return ret;
 93 }
 94 
 95 ExtInt operator ^(ExtInt a,int k){
 96     ExtInt ret(1);
 97     while(k){
 98         if(k&1)ret=ret*a;
 99         k>>=1;a=a*a;    
100     }
101     return ret;
102 }
103 
104 ExtInt operator /(ExtInt a,int k){
105     for(int i=a.len,tot=0;i>=1;i--){
106         tot=tot*10000+a[i];
107         a.num[i]=tot/k;
108         tot%=k;
109     }
110     while(!a[a.len])
111         a.num[a.len--]=0;
112     return a;
113 }
114 
115 int GCD(int a,int b){
116     return b?GCD(b,a%b):a;
117 }
118 
119 int main(){
120     scanf("%d",&n);
121     Linear_Shaker();
122     ExtInt ans(0);
123     for(int d=1;d<=n;d++)
124         if(n%d==0){
125             ExtInt x(2);
126             ans=ans+(x^d)*phi[n/d];
127         }
128     ans=ans/n;    
129     ans.Prf();
130     return 0;
131 }