题意
求下列表达式的值
\(\large{g^{\sum_{d|n}{\binom{n}{d}}} \pmod{999911659} }\)
其中，\(n, d \leqslant 10^9.\)
Solution
由欧拉定理可知，
\(\large{ 原式 = g^{\sum_{d|n}{\binom{n}{d}} \pmod{999911658}} }\)
显然只需要考虑分子，考虑到 \(999911658\) 范围下的组合数无法解决，考虑到 \(999911658=2*3*4679*35617\) ，为square-free number，则根据中国剩余定理，考虑分别求出在质数因子下的结果，
$$\begin
x \equiv \sum_{d|n}{\binom} \pmod{2} \
x \equiv \sum_{d|n}{\binom} \pmod{3} \
x \equiv \sum_{d|n}{\binom} \pmod{4679} \
x \equiv \sum_{d|n}{\binom} \pmod{35617} \
\end$$
特解 \(x\) 即为所求。
而对于里面每一个同余方程右侧的结果，Lucas定理可以比较优秀的解决，只需要预处理出阶乘和阶乘的逆元，利用Lucas定理即可求解。
Tips
在Lucas定理运算过程中，可能会出现 $d > n$的情况，需要特判出该种情况下的解 (为 \(0\) )，才能避免RE解出答案。
此外，当 \(g\) 为 \(999911659\) 的倍数时，答案为 \(0\)，不能使用欧拉定理进行计算。
Code
点击查看代码
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
const int N = 4e4 + 50;
const int mod = 999911659;
inline int read() {
register int w = 0, f = 1;
register char c = getchar();
while (c > '9' || c < '0') {
if (c == '-') f = -1;
c = getchar();
}
while (c >= '0' && c <= '9') {
w = w * 10 + c - '0';
c = getchar();
}
return w * f;
}
int n, g;
int p;
inline int qpow(register int a, register int b) {
register int base = 1;
while (b) {
if (b & 1) base = 1ll * base * a % p;
a = 1ll * a * a % p;
b >>= 1;
}
return base;
}
int mul[N], inv[N];
inline int Comb(register int n, register int m) {
if (n < m) return 0;
if (m == 0) return 1;
if (m < p && n < p) return 1ll * mul[n] * inv[m] % p * inv[n - m] % p;
return 1ll * Comb(n / p, m / p) * Comb(n % p, m % p) % p;
}
int a[5], m[5];
inline int exgcd(register int a, register int b, int &x, int &y) {
if (!b) {
x = 1;
y = 0;
return a;
}
register int d = exgcd(b, a % b, x, y);
register int tmp = x;
x = y;
y = tmp - a / b * y;
return d;
}
int a1, m1;
int main() {
n = read(), g = read();
if (g % mod == 0) {
printf("0\n");
return 0;
}
mul[0] = inv[0] = 1;
m[1] = p = 2;
for (register int i = 1; i <= p; ++i) {
mul[i] = 1ll * mul[i - 1] * i % p;
inv[i] = qpow(mul[i], p - 2);
}
for (register int d = 1; d * d <= n; ++d)
if (n % d == 0) {
a[1] = (a[1] + Comb(n, d)) % p;
if (d * d != n) a[1] = (a[1] + Comb(n, n / d)) % p;
}
m[2] = p = 3;
for (register int i = 1; i <= p; ++i) {
mul[i] = 1ll * mul[i - 1] * i % p;
inv[i] = qpow(mul[i], p - 2);
}
for (register int d = 1; d * d <= n; ++d)
if (n % d == 0) {
a[2] = (a[2] + Comb(n, d)) % p;
if (d * d != n) a[2] = (a[2] + Comb(n, n / d)) % p;
}
m[3] = p = 4679;
for (register int i = 1; i <= p; ++i) {
mul[i] = 1ll * mul[i - 1] * i % p;
inv[i] = qpow(mul[i], p - 2);
}
for (register int d = 1; d * d <= n; ++d)
if (n % d == 0) {
a[3] = (a[3] + Comb(n, d)) % p;
if (d * d != n) a[3] = (a[3] + Comb(n, n / d)) % p;
}
m[4] = p = 35617;
for (register int i = 1; i <= p; ++i) {
mul[i] = 1ll * mul[i - 1] * i % p;
inv[i] = qpow(mul[i], p - 2);
}
for (register int d = 1; d * d <= n; ++d)
if (n % d == 0) {
a[4] = (a[4] + Comb(n, d)) % p;
if (d * d != n) a[4] = (a[4] + Comb(n, n / d)) % p;
}
a1 = a[1], m1 = m[1];
for (register int i = 2; i <= 4; ++i) {
register int x, y;
register int g = exgcd(m1, m[i], x, y);
register int t = m[i] / g;
x = (x % t + t) % t;
register int tmp = m1 / g * m[i];
a1 = (a1 + 1ll * (a[i] - a1) / g * x % tmp * m1 % tmp + tmp) % tmp;
m1 = tmp;
}
a1 = (a1 % (mod - 1) + (mod - 1)) % (mod - 1);
p = mod;
printf("%d\n", qpow(g, a1));
return 0;
}