输入文件只有一行，用空格隔开的三个整数n, m, p 。保证有m <=n,p为素数。
#include<cstdio> typedef unsigned long long u64; typedef unsigned int u32; int n,m; u32 P; int gcd(int a,int b){ for(int c;b;c=a,a=b,b=c%b); return a; } int phi(int n){ int v=n; for(int i=2;i*i<=n;++i)if(n%i==0){ do n/=i;while(n%i==0); v=v/i*(i-1); } if(n>1)v=v/n*(n-1); return v; } inline u32 fix(int a){ return a+(a>>31&P); } struct num{ u32 x; num(u32 a=0):x(a){} num operator+(num w){return fix(x+w.x-P);} num operator*(num w){return u64(x)*w.x%P;} void operator+=(num w){x=fix(x+w.x-P);} }; num s[5][1007],gs[11],iv[117],f0[57][1007],f1[57][1007],ans; void cal(int m,int n){ int g=gcd(n,m); num v=0; for(int d=1;d<=g;++d)if(g%d==0)v+=f0[m/d][n/d]*phi(d); v+=f1[m][n]*m; // printf("%d,%d:%d\n",m,n,v*iv[m*2]); ans+=v*iv[m*2]; } int main(){ scanf("%d%d%u",&n,&m,&P); if(m>n)m=n; s[0][0]=iv[1]=1; for(int i=2;i<=115;++i)iv[i]=iv[P%i]*(P-P/i); for(int i=1;i<=n;++i){ f0[1][i]=f1[1][i]=s[0][i-1]+s[1][i-1]+s[2][i-1]; gs[1]=f0[1][i]+s[3][i-1]; // printf("[%d:%d]\n",i,gs[1]); for(int j=2;j<=3;++j)gs[j]=gs[j-1]*(gs[1]+(j-1))*iv[j]; for(int j=3;j;--j){ for(int k=n;k>=i;--k){ for(int t=1;t<=j;++t){ int w=k-t*i; if(w>=0)s[j][k]+=gs[t]*s[j-t][w]; } } } } for(int i=2;i<=m;++i){ for(int j=i;j<=n;++j){ for(int k=1;k<j;++k)f0[i][j]+=f0[i-1][j-k]*f0[1][k]; if(i&1){ for(int k=2;k<j;k+=2)f1[i][j]+=f0[i>>1][k>>1]*f0[1][j-k]; }else{ for(int k=1;k<j;++k)f1[i][j]+=f1[i-1][j-k]*f0[1][k]; if(~j&1)f1[i][j]+=f0[i>>1][j>>1]; f1[i][j]=f1[i][j]*iv[2]; } // printf("%d,%d %d %d\n",i,j,f0[i][j],f1[i][j]); } } for(int i=3;i<=m;++i)cal(i,n); printf("%d\n",ans.x); return 0; }
原文：http://www.cnblogs.com/ccz181078/p/7147947.html