BZOJ1951: [Sdoi2010]古代猪文

Posted 2020-08-31

传送门

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

很容易得到答案$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;
}