题目大意:给出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 #include2 #include 3 #include 4 #include 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 3600011 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);TMP64 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 }