BZOJ1951: [Sdoi2010]古代猪文

题目：http://www.lydsy.com/JudgeOnline/problem.php?id=1951

题意：

题解：我们发现模数是个质数，所以我们只需要求出指数 mod 999911659-1 的结果即可。

          这是个合数=2*3*4679*35617  所以我们就可以分开来做然后CRT合并。

          分开直接lucas即可。

代码：

#include<cstdio>

#include<cstdlib>

#include<cmath>

#include<cstring>

#include<algorithm>

#include<iostream>

#include<vector>

#include<map>

#include<set>

#include<queue>

#include<string>

#define inf 1000000000

#define maxn 200000+5

#define maxm 200000+5

#define eps 1e-10

#define ll long long

#define pa pair<int,int>

#define for0(i,n) for(int i=0;i<=(n);i++)

#define for1(i,n) for(int i=1;i<=(n);i++)

#define for2(i,x,y) for(int i=(x);i<=(y);i++)

#define for3(i,x,y) for(int i=(x);i>=(y);i--)

#define for4(i,x) for(int i=head[x],y=e[i].go;i;i=e[i].next,y=e[i].go)

#define for5(n,m) for(int i=1;i<=n;i++)for(int j=1;j<=m;j++)

using namespace std;

inline int read()

{

    int x=0,f=1;char ch=getchar();

    while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}

    while(ch>='0'&&ch<='9'){x=10*x+ch-'0';ch=getchar();}

    return x*f;

}
const int p[5]={2,3,4679,35617,999911658};
int n,m,fac[maxn],inv[maxn];
inline ll power(ll x,ll y,ll p)
{
    ll t=1;
    for(;y;y>>=1,x=x*x%p)
        if(y&1)t=t*x%p;
    return t;
}
inline void exgcd(ll a,ll b,ll &x,ll &y)
{
    if(!b){x=1;y=0;return;}
    exgcd(b,a%b,x,y);
    int t=x;x=y;y=t-a/b*y;
}
inline ll invinv(ll a,ll b)
{
    ll x,y;
    exgcd(a,b,x,y);
    return (x%b+b)%b;
}
inline int c(int n,int m,int p)
{
    if(n<m)return 0;
    if(n<p&&m<p)return fac[n]*inv[m]%p*inv[n-m]%p;
    return c(n/p,m/p,p)*c(n%p,m%p,p)%p;
}

int main()

{

    freopen("input.txt","r",stdin);

    freopen("output.txt","w",stdout);

    n=read();m=read();ll ans=0;
    m%=p[4]+1;
    if(!m){printf("0
");return 0;}
    for0(i,3)
    {
        fac[0]=1;
        for1(j,p[i]-1)fac[j]=fac[j-1]*j%p[i];
        inv[0]=inv[1]=1;
        for2(j,2,p[i]-1)inv[j]=(p[i]/j+1)*inv[j-p[i]%j]%p[i];
        for2(j,2,p[i]-1)inv[j]=inv[j]*inv[j-1]%p[i];
        ll t1=0,t2=invinv(p[4]/p[i],p[i]);
        for(int j=1;j*j<=n;j++)if(n%j==0)
        {
            (t1+=c(n,j,p[i]))%=p[i];
            if(j*j==n)continue;
            (t1+=c(n,n/j,p[i]))%=p[i];
        }
        (ans+=(ll)p[4]/p[i]%p[4]*t2%p[4]*t1%p[4])%=p[4];
    }
    (ans+=p[4])%=p[4];
    cout<<power(m,ans,p[4]+1)<<endl;

    return 0;

}