題目

Luogu P2480

求值：

\[G^{\sum_{k\mid n}\binom{n}{k}} \bmod 999911659 \]

分析

若 \(999911659\mid G\) ，顯然答案為 \(0\)

若 \(999911659\not\mid G\) ，因為 \(999911659\) 是素數，由尤拉定理可得：

\[G^{\sum_{k\mid n}\binom{n}{k}} \equiv G^{\sum_{k\mid n}\binom{n}{k}\bmod 999911658}\pmod{999911659} \]

我們先考慮求出下式：

\[\sum_{k\mid n}\binom{n}{k}\bmod 999911658 \]

注意到 \(999911658=2\times3\times4679\times35617\)

所以可以分別求出模 \(2,3,4679,35617\) 的值，計算組合數時需要用到 盧卡斯定理

將結果記為 \(b_1,b_2,b_3,b_4\) 再用 中國剩餘定理 合併即可：

\[\begin{cases} x \equiv b_1 \pmod{2}\\ x \equiv b_2 \pmod{3}\\ x \equiv b_3 \pmod{4679}\\ x \equiv b_4 \pmod{35617} \end{cases} \]

得到結果後再對 \(G\) 進行一次快速冪即為答案

程式碼

#include<bits/stdc++.h>
#define ll long long
using namespace std;

const int MOD = 999911659;
int n, g;
ll fact[50000 + 5];
int num[4] = {2, 3, 4679, 35617};
ll b[5];
vector<int> d;

ll qpow(ll a, ll b, ll p)
{
    if(b == 0) 
        return 1;
    ll res = qpow(a, b / 2, p);
    if(b % 2)
        return res * res % p * a % p;
    return res * res % p;
}

ll C(ll n, ll m, ll p)
{
    if(n < 0 || m < 0 || m > n)
        return 0;
    return fact[n] * qpow(fact[m], p - 2, p) % p * qpow(fact[n - m], p - 2, p) % p;
}

ll lucas(ll n, ll m, ll p)
{
    if(n < 0 || m < 0 || m > n)
        return 0;
    if(m == 0)
        return 1;
    return lucas(n / p, m / p, p) * C(n % p, m % p, p) % p;
}

ll exgcd(ll a, ll b, ll & x, ll & y)
{
    if(b == 0) {
        x = 1;
        y = 0;
        return a;
    }
    ll d = exgcd(b, a % b, y, x);
    y -= a / b * x;
    return d;
}

ll inv(ll a, ll p)
{
    ll x, y, d;
    d = exgcd(a, p, x, y);
    ll t = p / d;
    return (x % t + t) % t;
}

ll CRT()
{
    ll M = 999911658;
    ll ans = 0;
    for(int i = 0; i < 4; i++) {
        ll m = M / num[i];
        ll t = inv(m, num[i]);
        ans = (ans + m * t % M * b[i] % M) % M;
    }
    return ans;
}

int main()
{
    cin >> n >> g;
    if(g == MOD) {
        cout << 0 << endl;
        return 0;
    }
    int i;
    for(i = 1; i * i < n; i++) {
        if(n % i == 0) {
            d.push_back(i);
            d.push_back(n / i);
        }
    }
    if(i * i == n)
        d.push_back(i);
    for(i = 0; i < 4; i++) {
        ll mod = num[i];
        fact[0] = 1;
        for(int j = 1; j <= mod; j++)
            fact[j] = fact[j - 1] * j % mod;
        b[i] = 0;
        for(int j = d.size() - 1; j >= 0; j--) {
            b[i] = (b[i] + lucas(n, d[j], mod)) % mod;
        }
    }
    cout << qpow(g, CRT(), MOD) << endl;
    return 0;
}