[SDOI2010] 古代猪文

题意

求下列表达式的值

$$g^{\sum _{d|n}(\genfrac{}{}{0ex}{}{n}{d})}(\mathrm{mod}999911659)$$

其中，$n,d⩽{10}^9.$

Solution

由欧拉定理可知，

$$\mathrm{原}\mathrm{式}=g^{\sum _{d|n}(\genfrac{}{}{0ex}{}{n}{d})(\mathrm{mod}999911658)}$$

显然只需要考虑分子，考虑到 $999911658$ 范围下的组合数无法解决，考虑到 $999911658=2\ast 3\ast 4679\ast 35617$ ，为square-free number，则根据中国剩余定理，考虑分别求出在质数因子下的结果，

$$\{\begin{array}{l}x\equiv \sum _{d|n}(\genfrac{}{}{0ex}{}{n}{d})(\mathrm{mod}2)\\ x\equiv \sum _{d|n}(\genfrac{}{}{0ex}{}{n}{d})(\mathrm{mod}3)\\ x\equiv \sum _{d|n}(\genfrac{}{}{0ex}{}{n}{d})(\mathrm{mod}4679)\\ x\equiv \sum _{d|n}(\genfrac{}{}{0ex}{}{n}{d})(\mathrm{mod}35617)\end{array}$$

特解 $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;
}