首页 > 代码库 > BZOJ1951: [Sdoi2010]古代猪文

BZOJ1951: [Sdoi2010]古代猪文

传送门

数论的套路似乎没多少..

很容易得到答案$ans=G^{C(N,i)}$其中$i$为$N$的因数。

很显然,我们可以在$O(logN)$的时间内算出所有的因数,但是怎么算大组合数呢?大力Lucas。

然而指数应该模$P-1$,$P-1$不是质数啊,这就很痛苦啊。

通过暴力计算法我们可以发现$P-1$是4个质数的合。我们可以考虑分别计算答案然后中国剩余定理合并答案即可。

//BZOJ 1951//by Cydiater//2017.2.18#include <iostream>#include <queue>#include <map>#include <cstring>#include <string>#include <algorithm>#include <iomanip>#include <cstdlib>#include <cstdio>#include <ctime>#include <cmath>#include <bitset>#include <set>#include <vector>#include <complex>using namespace std;#define ll long long#define up(i,j,n)   for(ll i=j;i<=n;i++)#define down(i,j,n) for(ll i=j;i>=n;i--)#define cmax(a,b)   a=max(a,b)#define cmin(a,b)   a=min(a,b)const ll MAXN=2e5+5;const ll oo=1LL<<50;const ll LIM=1e5+5;const ll mod=999911659;inline ll read(){    char ch=getchar();ll x=0,f=1;    while(ch>‘9‘||ch<‘0‘){if(ch==‘-‘)f=-1;ch=getchar();}    while(ch>=‘0‘&&ch<=‘9‘){x=x*10+ch-‘0‘;ch=getchar();}    return x*f;}ll N,G,m[4]={2,3,4679,35617},inv[MAXN][4],fac[MAXN][4],a[4];namespace solution{    void exgcd(ll a,ll b,ll &x,ll &y){        if(!b){x=1;y=0;return;}        exgcd(b,a%b,y,x);        y-=a/b*x;    }    ll quick_pow(ll base,ll ind){        ll tmp=1;        while(ind){            if(ind&1)(tmp*=base)%=mod;            (base*=base)%=mod;ind>>=1;        }        return tmp;    }    ll C(ll i,ll j,ll p){        if(j>i)return 0;        if(i<m[p]&&j<m[p])return fac[i][p]*inv[j][p]%m[p]*inv[i-j][p]%m[p];        return C(i/m[p],j/m[p],p)*C(i%m[p],j%m[p],p)%m[p];    }    void Prepare(){        N=read();G=read();        fac[0][0]=fac[0][1]=fac[0][2]=fac[0][3]=1;        up(i,1,LIM)up(j,0,3)fac[i][j]=fac[i-1][j]*i%m[j];        inv[1][0]=inv[1][1]=inv[1][2]=inv[1][3]=1;        up(i,2,LIM)up(j,0,3)inv[i][j]=(m[j]-m[j]/i)*inv[m[j]%i][j]%m[j];        inv[0][0]=inv[0][1]=inv[0][2]=inv[0][3]=1;        up(i,1,LIM)up(j,0,3)(inv[i][j]*=inv[i-1][j])%=m[j];    }    ll CRT(){        ll a1=a[0],m1=m[0];        up(i,1,3){            ll a2=a[i],m2=m[i],x,y;            exgcd(m1,m2,x,y);            x=(x*(a2-a1)%m2+m2)%m2;            a1+=x*m1;            m1*=m2;            a1%=m1;        }        return a1;    }    void Solve(){        if(G==mod){            puts("0");            return;        }        for(int i=1;i*i<=N;i++)if(!(N%i))up(j,0,3){            (a[j]+=C(N,i,j))%=m[j];            if(i*i!=N)(a[j]+=C(N,N/i,j))%=m[j];         }        ll tmp=CRT();        //cout<<tmp<<endl;        printf("%lld\n",quick_pow(G,tmp));    }}int main(){    //freopen("input.in","r",stdin);    using namespace solution;    Prepare();    Solve();    return 0;}

 

BZOJ1951: [Sdoi2010]古代猪文