「SDOI2010」古代猪文 - 数论 + Lucas 定理

发表于 更新于 分类于 ， 阅读次数： Disqus：

题意描述

洛谷链接

LibreOJ 链接

给定整数 \(g\)，\(n\)（\(1 \leq q, n \leq 10^9\)），计算：

\[ g^{\sum_{d|n}{C_{n}^{d}}} \operatorname{mod} 999911659 \]

解题思路

\(999911659\) 为质数，显然与 \(g\) 互质。

对于该式，我们可以用欧拉定理推论转化为：

\[ g^{\sum_{d|n}{C_{n}^{d}} \bmod \varphi(999911659)} \bmod 999911659 \\ = g^{\sum_{d|n}{C_{n}^{d}} \bmod 999911658} \bmod 999911659 \]

于是该题可转化为求出 \(\sum_{d|n}{C_{n}^{d}} \operatorname{mod} 999911658\)。

我们将 \(999911658\) 分解质因数，可得 \(888811658 = 2 \times 3 \times 4679 \times 35617\)，其中没有相同的质因数。

于是对于方程 \(x \equiv \sum_{d|n}{C_{n}^{d}} \pmod {999911658}\)，显然对于下列同余方程组：

\[ \begin{cases} x \equiv \sum_{d|n}{C_{n}^{d}} \pmod {2} \\ x \equiv \sum_{d|n}{C_{n}^{d}} \pmod {3} \\ x \equiv \sum_{d|n}{C_{n}^{d}} \pmod {4679} \\ x \equiv \sum_{d|n}{C_{n}^{d}} \pmod {35617} \end{cases} \]

其中 \(x\) 的解即为方程 \(x \equiv \sum_{d|n}{C_{n}^{d}} \pmod {999911658}\) 的解。

我们先对 \(n\) 分解质因数，再用中国剩余定理求出该方程组的解，用 Lucas 定理和预处理阶乘及逆元优化求解 \(C_{n}^{d} \bmod p\)。最终利用快速幂求出 \(g^x \bmod 999911659\) 即得出答案。

代码演示

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
#include <iostream>

const int MAXN = 35617, MOD = 999911659;

static long long fact[MAXN + 1], inv[MAXN + 1];

long long power(long long x, long long y, long long p) {
    long long re = 1;
    while (y) {
        if (y & 1) re = (re * x) % p;
        x = (x * x) % p;
        y >>= 1;
    }
    return re;
}

long long C(long long x, long long y, long long p) {
    if (x < y) return 0;
    return ((fact[x] * inv[y]) % p * inv[x - y]) % p;
}

long long lucas(long long n, long long m, long long p) {
    if (n < m) return 0;
    if (!n) return 1;
    return (lucas(n / p, m / p, p) * C(n % p, m % p, p)) % p;
}

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

long long crt(int k, long long* a, long long* r) {
    long long n = 1, ans = 0;
    for (int i = 1; i <= k; i++) n = n * r[i];
    for (int i = 1; i <= k; i++) {
        long long m = n / r[i], b, y;
        exgcd(m, r[i], b, y);
        ans = (ans + a[i] * m * b % n) % n;
    }
    return (ans % n + n) % n;
}

void init(int p) {
    fact[0] = 1;
    for (long long i = 1; i <= p; i++) {
        fact[i] = (fact[i - 1] * i) % p;
    }

    inv[p - 1] = power(fact[p - 1], p - 2, p);
    for (long long i = p - 2; i >= 0; i--) {
        inv[i] = (inv[i + 1] * (i + 1)) % p;
    }
}

int main() {
    long long n, g;

    std::cin >> n >> g;

    if (g == MOD) {
        std::cout << "0" << std::endl;
        return 0;
    }

    static long long f[5] = { -1, 2, 3, 4679, 35617 }, a[5];

    for (int i = 1; i <= 4; i++) {
        init(f[i]);

        for (long long d = 1; d * d <= n; d++) {
            if (n % d == 0) {
                a[i] = (a[i] + lucas(n, d, f[i])) % f[i];
                if (d * d != n) {
                    a[i] = (a[i] + lucas(n, n / d, f[i])) % f[i];
                }
            }
        }
    }

    std::cout << power(g, crt(4, a, f), MOD) << std::endl;

    return 0;
}