BZOJ1951 古代猪文

题目链接：

题意：

计算\(g^{\sum_{k|n}(^n_k)}\%999911659\)

\(n\le 10^9, g\le 10^9\)

题解：

首先，根据扩展欧拉定理，\(a^b≡a^{b\%\phi(p)}\ (MOD\ p), gcd(a,p)=1\)

可以把要计算的式子降幂得到：\(g^{(\sum_{k|n}(^n_k))\%999911658}\%999911659\)

接下来我们需要计算的就是\((\sum_{k|n}(^n_k))\%999911658\)

但是要是直接计算是不现实的，一是组合数要用到\(n\)的阶乘，而\(n\)有\(10^9\)的大小，其次要算的阶乘的逆元也不一定存在

可以发现\(999911658\)可以拆分成\(2\times 3\times 4679\times 35617\)

所以我们可以计算出\(\sum_{k|n}(^n_k)\)在四个质数下的模数，然后用中国剩余定理来合并

也就是说，求出来在四个模数下的值\(r_1,r_2,r_3,r_4\)之后，得到一个线性方程组：

\[\begin{cases}
x≡r_1\ (MOD\ 2) \\
x≡r_2\ (MOD\ 3) \\
x≡r_3\ (MOD\ 4679) \\
x≡r_4\ (MOD\ 35617)
\end{cases}\]

用中国剩余定理合并一下就可以得到\(x\)在模\(999911658\)下的值了

然后用快速幂就可以计算出最后的答案

而拆分后的质数都比较小，所以可以用卢卡斯定理来算在小质数下的组合数

view code

//#pragma GCC optimize("O3")
//#pragma comment(linker, "/STACK:1024000000,1024000000")
#include<bits/stdc++.h>
using namespace std;
typedef long long int LL;
const LL MOD = 999911659LL;
const int MAXN = 4e4+7;
LL fact[MAXN];
vector<LL> P,F;
LL n, g;
void preprocess(){
    LL x = MOD - 1;
    for(int i = 2; i * i <= x; i++){
        if(x%i==0){
            P.push_back(i);
            while(x%i==0) x /= i;
        }
    }
    if(x!=1) P.push_back(x);
    for(int i = 1; i * i <= n; i++){
        if(n%i==0){
            F.push_back(i);
            if(i!=n/i) F.push_back(n/i);
        }
    }
}
LL ksm(LL a, LL b, LL p){                              // a^b%p
    LL ret = 1;
    while(b){
        if(b&1) ret = ret * a % p;
        b >>= 1;
        a = a * a % p;
    }
    return ret;
}
void exgcd(LL a, LL b, LL &x, LL &y){                   //ax + by = 1
    if(!b){ x = 1, y = 0; return; }
    exgcd(b,a%b,y,x);
    y -= a / b * x;
}
void calComb(LL p){                                     // 预处理组合数
    fact[0] = 1;
    for(int i = 1; i < p; i++) fact[i] = fact[i-1] * i % p;
}
LL inv(LL a, LL p){                                     // a在模p意义下的逆元
    LL x, y;
    exgcd(a,p,x,y);
    return (x%p)+p;
}
LL C(LL n, LL m, LL p){ return n<m?0:fact[n] * inv(fact[m],p) % p * inv(fact[n-m],p) % p; }       // 组合数C(n,m)%p
LL lucas(LL n, LL m, LL p){             // 卢卡斯计算 C(n,m) % p
    LL ret = 1;
    while(n and m){
        int nn = n % p, mm = m % p;
        if(mm>nn) return 0LL;
        ret = ret * C(nn,mm,p) % p;
        n /= p; m /= p;
    }
    return ret;
}
LL CRT(vector<LL> R, vector<LL> P){
    LL ret = 0;
    for(int i = 0; i < (int)R.size(); i++)
        ret = (ret + R[i] * (MOD-1)/P[i] % (MOD-1) * inv((MOD-1)/P[i],P[i])) % (MOD-1);
    return ret;
}
LL solve(){
    preprocess();
    vector<LL> R;
    for(int p : P){
        calComb(p);
        LL r = 0;
        for(int f : F) r = (r + lucas(n,f,p)) % p;
        R.push_back(r);
    }
    LL pw = CRT(R,P);
    return ksm(g,pw,MOD);
}
int main(){
    scanf("%lld %lld",&n,&g);
    g %= MOD;
    if(g == 0) cout << 0 << endl;
    else if(g == 1) cout << 1 << endl;
    else cout << solve() << endl;
    return 0;
}