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 }
View Code

 

posted @ 2016-12-13 21:57  Blue233333  阅读(120)  评论(0编辑  收藏  举报