标签:
数论综合。费马小定理,lucas定理,中国剩余定理,exgcd,快速幂,乘法逆元。
首先要计算出n的每个约数,简单的sqrt(n)枚举即可。
然后计算C(i,m)(m个中挑i个的组合数,ps:因为网上正反俩种都有,所以标注一下。。)
设s=sum(C(i,m))
题目要求g^(s)%mod,
由费马小定理得 g^(s)=g^(s%(mod-1))。所以这里需要特判一下g是否等于mod。
可是组合数计算因为用到除法,所以模数必须为质数。999911659-1=2×3×4679×35617。
这四个数设为d[i],对这四个数分别进行计算。
由因为组合数很大,直接计算会超时。
需要用到lucas定理。 lucas(n,m)=lucas(n/mod,m/mod)*C(n%mod,m%mod)。这里mod不是999911659,而是那4个数。
所以就得到4个答案m[i]。
组合数的计算可以根据费马小定理转换成快速幂计算。
然后再用中国剩余定理合并。
这时得到的答案可以写成如下形式:
s%d[0]=m[0],s%d[1]=m[1];
改写为 s%a1=b1,s%a2=b2 (1)
则有 a1*x+b1=a2*y+b2 (x,y分别为俩个未知数,和上面的任何数没有任何关系)
设a=a1,b=a2,c=b2-b1
则上式变为 ax+by=c (正负号并没有关系,因为不会用到y)
计算t=extended_gcd(a,b,x,y)
如果c%t !=0 的话,说明此方程是无解的。在本题中t=1。(一下都基于本题t=1的情况)
然后我们计算出的是 ax+by=t。 所以我们将每一项乘以c。
x=(x*c%b+b)%b。
当然这样的x有无数组为 x1=x+k*b。
带入(1)式有 s=a1(x+k*b)+b1。
s=a1*b*k+a1*x+b1.
所以有 s%(a1*b)=a1*x+b1。
所以 b1=b1+a1*x,a1=a1*b。
在本题中这样操作完以后a1>b1,实际上即使t!=1,b1也会小于1。
最后合并完返回b1即可。
#include<cstdio> #include<algorithm> #include<cstring> #include<cmath> #define LL long long using namespace std; const int mod = 999911659; const int d[]={2,3,4679,35617}; int fac[4][40000],m[4]; int n,g; int power(LL a,int e,int mod) { LL res=1; while(e) { if(e&1) res=res*a%mod; a=a*a%mod; e>>=1; } return res; } int C(int n,int m,int x) { if(n<m) return 0; return fac[x][n]*power(fac[x][m]*fac[x][n-m]%d[x],d[x]-2,d[x])%d[x]; } int lucas(int n,int m,int x) { if(m==0) return 1; return lucas(n/d[x],m/d[x],x)*C(n%d[x],m%d[x],x)%d[x]; } int exgcd(int a,int b,int &x,int &y) { if(b==0) { x=1; y=0; return a; } else { int t=exgcd(b,a%b,y,x); y-=(a/b)*x; return t; } } int solve() { int a1,b1,a2,b2,a,b,c,x,y; a1=d[0]; b1=m[0]; for(int i=1;i<4;i++) { a2=d[i];b2=m[i]; a=a1,b=a2,c=b2-b1; exgcd(a,b,x,y); x=(c*x%b+b)%b;//if(c%t) printf("Test\n"); b1=b1+a1*x; a1=a1*b; } return b1; } int main() { for(int i=0;i<4;i++) { fac[i][0]=1; for(int j=1;j<=40000;j++) fac[i][j]=fac[i][j-1]*j%d[i]; } scanf("%d%d",&n,&g); if(g==mod) printf("0\n"); else { g=g%mod; int n2=(int)sqrt(n+0.5); for(int i=1,j;i<=n2;i++) if(n%i==0) { j=n/i; for(int x=0;x<4;x++) { m[x]=(m[x]+lucas(n,i,x))%d[x]; if(i!=j) m[x]=(m[x]+lucas(n,j,x))%d[x]; } } printf("%d\n",power(g,solve(),mod)); } return 0; }
标签:
原文地址:http://www.cnblogs.com/invoid/p/5645307.html