BZOJ1951 [Sdoi2010]古代猪文
题目大意：给出n,G,求G^P%M，其中M=999911659，P=∑C(n,i) (i=1~n,n%i=0)。
题解：此神题也。注意到M是个素数，则G^P%M=G^(P%phi(M)+phi(M))%M=G^(P%(M-1)+M-1)%M，所以求P%(M-1)即可。错误！前式应满足条件P>=phi(M)即P>=M-1，如果P<M-1，可以随便举出一个反例。所以，判断P是否≥M-1，若不是直接求答案，若是就进行下一步。令人头疼的是M-1不是一个素数，求P%(M-1)就要用Lucas求∑C(n,i) (i=1~n,n%i=0)，其中求C需要乘法逆元，而面对模数M-1不是素数的情况不能直接求。于是，把M-1分解质因数，用Lucas(n,i,p)(p为质因数)求出几个模数后，解一波线性模方程组，用中国剩余定理即可。Lucas(a,b,p)=C(a%p,b%p,p)*Lucas(a/p,b/p,p)%p，C的求法见代码。最后的答案跑快速幂。
代码：
1 #include<cstdio> 2 #include<cstring> 3 #include<cstdlib> 4 #include<algorithm> 5 using namespace std; 6 7 #define LL long long 8 #define M 999911659 9 LL son[]={0,2,3,4679,35617}; 10 #define maxs 36000 11 LL n,g,fac[maxs][5],p,mo[5]; 12 LL pow_mod(LL a,LL b,LL p) 13 { 14 LL tmp=a%p,ans=1; 15 while (b) {if (b&1) ans=ans*tmp%p;tmp=tmp*tmp%p;b>>=1;} 16 return ans; 17 } 18 void ext_gcd(LL a,LL b,LL& x,LL& y) 19 { 20 if (!b) {x=a;y=0;} 21 else ext_gcd(b,a%b,y,x),y-=a/b*x; 22 } 23 LL C(LL a,LL b,LL pos) 24 { 25 if (b>a) return 0;if (b==0 || b==a) return 1;LL p=son[pos]; 26 return fac[a][pos]*pow_mod(fac[b][pos]*fac[a-b][pos],p-2,p)%p; 27 } 28 LL Lucas(LL a,LL b,LL pos) 29 { 30 if (!b) return 1;LL p=son[pos]; 31 return (C(a%p,b%p,pos)*Lucas(a/p,b/p,pos)%p); 32 } 33 LL china(LL n,LL* a,LL* m) 34 { 35 LL x,y,ans=0; 36 for (int i=1;i<=n;i++) 37 { 38 LL w=(M-1)/m[i]; 39 ext_gcd(w,m[i],x,y); 40 ans=(ans+w*x*a[i])%(M-1); 41 } 42 return (ans+(M-1))%(M-1); 43 } 44 LL C2(LL a,LL b) 45 { 46 if (b>a) return 0;b=min(a-b,b);if (b>25) return M; 47 LL ans=1; 48 for (int i=1;i<=b;i++) {ans=ans*(a-i+1)/i;if (ans>=M-1) break;} 49 if (ans>=M-1) return M;return ans; 50 } 51 int main() 52 { 53 scanf("%lld%lld",&n,&g); 54 for (int j=1;j<=4;j++) fac[0][j]=1; 55 for (int j=1;j<=4;j++) 56 for (int i=1;i<=son[j];i++) 57 fac[i][j]=fac[i-1][j]*i%son[j]; 58 LL tmp=0;int ok=0; 59 #define TMP if (tmp>=M-1) {ok=1;break;} 60 for (LL i=1;i*i<=n;i++) 61 if (!(n%i)) 62 { 63 tmp+=C2(n,i);TMP 64 if (i*i!=n) {tmp+=C2(n,n/i);TMP} 65 } 66 for (int j=1;j<=4;j++) 67 { 68 mo[j]=0; 69 for (LL i=1;i*i<=n;i++) 70 if (!(n%i)) 71 { 72 mo[j]=(mo[j]+Lucas(n,i,j))%son[j]; 73 if (i*i!=n) mo[j]=(mo[j]+Lucas(n,n/i,j))%son[j]; 74 } 75 } 76 p=ok?china(4,mo,son)+M-1:tmp; 77 printf("%lld",pow_mod(g,p,M)%M); 78 return 0; 79 }
