标签:pac href lin inline ++ 合并 print type lucas定理
考虑所求即为:\(G^{\sum\limits_{d|n}C_n^d}\%999911659\)。
发现系数很大,先用欧拉定理化简系数:\(G^{\sum\limits_{d|n}C_n^d\%999911658}\%999911659\)。
实际上我们只用求\(\sum\limits_{d|n}C_n^d\%999911658\),之后快速幂即可。
发现\(999911658\)不是个质数,没办法用Lucas定理求组合数,于是考虑拆开\(999911658\),发现为\(2,3,4679,35617\)。
于是对模意义下这四个数分别求\(\sum\limits_{d|n}C_n^d\),假设第\(i\)个求出的为\(a_i\)
发现我们得到了四个形如\(x\equiv a_i\pmod{p_i}\)的方程,用中国剩余定理合并即可得到答案(这其实就是exLucas的简化版)。
code:
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll mod1=999911659;
const ll mod2=999911658;
const int maxs=40010;
ll n,m;
ll a[5],prime[]={0,2,3,4679,35617};
ll fac[maxs][5],inv[maxs][5];
inline ll power(ll x,ll k,ll mod)
{
ll res=1;
while(k)
{
if(k&1)res=res*x%mod;
x=x*x%mod;k>>=1;
}
return res;
}
inline ll C(ll n,ll m,ll op)
{
if(m>n)return 0;
return fac[n][op]*inv[n-m][op]%prime[op]*inv[m][op]%prime[op];
}
inline ll Lucas(ll n,ll m,ll op)
{
if(!m)return 1;
return C(n%prime[op],m%prime[op],op)*Lucas(n/prime[op],m/prime[op],op)%prime[op];
}
void exgcd(ll a,ll b,ll& x,ll& y)
{
if(!b){x=1,y=0;return;}
exgcd(b,a%b,x,y);
ll z=x;x=y,y=z-(a/b)*y;
}
inline ll CRT()
{
ll res=0;
for(int i=1;i<=4;i++)
{
ll x,y,M=mod2/prime[i];
exgcd(M,prime[i],x,y);
x=(x%prime[i]+prime[i])%prime[i];
res=(res+a[i]*x%mod2*M%mod2)%mod2;
}
return res;
}
int main()
{
scanf("%lld%lld",&n,&m);
if(m==mod1){puts("0");return 0;}
for(int i=1;i<=4;i++)
{
fac[0][i]=1;
for(int j=1;j<prime[i];j++)fac[j][i]=fac[j-1][i]*j%prime[i];
inv[prime[i]-1][i]=power(fac[prime[i]-1][i],prime[i]-2,prime[i]);
for(int j=prime[i]-1;j;j--)inv[j-1][i]=inv[j][i]*j%prime[i];
}
for(ll i=1;i*i<=n;i++)
{
if(n%i)continue;
for(int j=1;j<=4;j++)
{
a[j]=(a[j]+Lucas(n,i,j))%prime[j];
if(i*i!=n)a[j]=(a[j]+Lucas(n,n/i,j))%prime[j];
}
}
printf("%lld",power(m,CRT(),mod1));
return 0;
}
标签:pac href lin inline ++ 合并 print type lucas定理
原文地址:https://www.cnblogs.com/nofind/p/11933166.html