ACWING 213 古代猪文

https://www.acwing.com/problem/content/215/

题意

给定整数$q,n(1\le q,n\le 10^9)$，计算

$$q^{{\sum }_{d|n}{C}_n^d}mod999911659$$

思路

首先，进行一个特判。当$q=999911659$，显然答案为0。
否则由于999911659是一个质数，故其与q互质，满足欧拉定理的条件。
之后，由欧拉定理的推论

若正整数$a,n$互质，则对于任意整数$a$，有$a^b\equiv a^{(bmod\phi (n))}(modn)$

可以得到

$$q^{{\sum }_{d|n}{C}_n^d}\equiv q^{{\sum }_{d|n}{C}_n^{d}mod999911658}(mod999911659)$$

注意到999911658显然是一个合数，我们尝试将它分解质因数得到

$$999911658=2\times 3\times 4679\times 35617$$

它们的次数全都是1。

---

接下来，由于问题要求我们求$d|n$，因此我们可以根据这些质因数枚举999911658的约数。对每一个约数，可以利用Lucas定理

$$C_n^m\equiv C_{nmodp}^{mmodp}\ast C_{\frac{n}{p}}^{\frac{m}{p}}(modp)$$

求出组合数$C_n^d$，并计算出$\sum d|n{C}_n^d$对四个质因数取模的结果，记为$a_1\sim a_4$

---

最后，利用中国剩余定理

设$m_1\sim m_n$是两两互质的整数，$m={\prod }_{i=1}^n{m}_i,M_i=\frac{m}{m_i}$，$t_i$是线性同余方程$M_i{t}_i\equiv 1(mod{m}_i)$。对于任意的$n$个整数$a_1\sim a_n$，方程组

$$\{\begin{array}{l}x\equiv a_1(mod{m}_1)\\ x\equiv a_2(mod{m}_2)\\ \dots \\ x\equiv a_n(mod{m}_n)\end{array}$$

有整数解为$x={\sum }_{i=1}^n{a}_i{M}_i{t}_i$

---

解出线性方程组

$$\{\begin{array}{l}xmod2=a_1\\ xmod3=a_2\\ xmod4679=a_3\\ xmod35617=a_4\end{array}$$

• 这一步是可以转化成CRT的目标形式，即$x\equiv a_i(mod{m}_i)$的。这是因为$xmod{m}_i=a_i$代表$a_{i}mod{m}_i\equiv a_i$，又$xmod{m}_i=a_i$。

得到${\sum }_{d|n}^{C_n^d}mod999911658$的最小非负整数解$x$，然后利用快速幂求$q^x$即可。

#include <iostream>
#include <cstring>
#include <algorithm>
using namespace std;
typedef long long ll;
const int a[4] = {2, 3, 4679, 35617};
int factorials[36000], b[4], n, g, ans, x, y, mod = 999911658; // b中存储了CRT右侧的数字

int power(int a, int b)
{
    int res = 1;
    while (b)
    {
        if (b & 1)
            res = (ll)res * a % mod;
        b >>= 1;
        a = (ll)a * a % mod;
    }
    return res;
}

void exgcd(int a, int b, int &x, int &y)
{
    if (!b)
    {
        x = 1, y = 0;
        return;
    }
    exgcd(b, a % b, x, y);
    int z = x;
    x = y;
    y = z - y * (a / b);
}

int inv(int a, int p) // 求逆元
{
    int x, y;
    exgcd(a, p, x, y);
    return (x % p + p) % p;
}

int calc(int x, int mod) // lucas定理求组合数
{
    int ans = 1, y, a, b;
    for (y = n; x; x /= mod, y /= mod)
    {
        a = x % mod, b = y % mod;                                                                                            // 余数
        ans = (ll)ans * factorials[b] % mod * inv(factorials[a], mod) % mod * inv(b < a ? 0 : factorials[b - a], mod) % mod; // 求m!/n!(m-n)!
    }
    return ans;
}

int main()
{
    cin >> n >> g;
    g %= mod + 1;
    if (!g)
    { // 特判999911659
        cout << 0 << endl;
        return 0;
    }
    for (int i = factorials[0] = 1; i <= a[3]; i++)
    {
        // 阶乘预处理
        factorials[i] = (ll)factorials[i - 1] * i % mod;
    }
    for (int i = 1; i * i <= n; i++) // 枚举质因数可构成的约数
    {
        if (n % i == 0) // i是约数
        {
            for (int j = 0; j < 4; j++)
                b[j] = (b[j] + calc(i, a[j])) % a[j]; // 求C_n^i
            if (i * i != n)
            {
                for (int j = 0; j < 4; j++)
                    b[j] = (b[j] + calc(n / i, a[j])) % a[j]; // 枚举约数 并计算相应的组合数
            }
        }
    }
    for (int i = 0; i < 4; i++)
    {
        exgcd(mod / a[i], a[i], x, y); // CRT
        ans = (ans + (ll)x * (mod / a[i]) % mod * b[i]) % mod;
    }
    ans = (ans + mod) % mod, mod++; // 映射到最小非负值
    ans = power(g, ans);            // q^x
    cout << ans << endl;
    return 0;
}