【题解】Luogu-P2480 [SDOI2010]古代猪文

P2480 [SDOI2010]古代猪文

Description

给定正整数 \(n,g\)，求 \(g^{\sum_{k\mid n}C_{n}^{k}}\bmod 999911659\)。

• 对于 \(100\%\) 的数据，\(1\le n,g \le 10^9\)。

Solution

前置芝士：

• 欧拉定理
• \(\rm Lucas\) 定理
• \(\rm CRT\)

由欧拉定理知 \(g^{\sum_{k\mid n}C_{n}^{k}}\equiv g^{(\sum_{k\mid n}C_{n}^{k})\bmod\ \varphi(999911659)}\equiv g^{(\sum_{k\mid n}C_{n}^{k})\bmod\ 999911658}\pmod{999911659}\)

指数可以快速幂 \(O(\log 999911658)\) 解决，找因数是 \(\Theta(\sqrt{n})\) 的。

于是问题转化为求 \(C_{n}^{k}\bmod 999911658\)。

直接算是不行的，但是直接 \(\rm Lucas\) 因为模数太大也不行。

\(999911658=2×3×4679×35617\)，令 \(x=C_{n}^{k}\)，考虑分别算出

\[\begin{cases} x\bmod 2\\ x\bmod 3\\ x\bmod 4679\\ x\bmod 35617 \end{cases} \]

这个就可以用 \(\rm Lucas\)。

再用 \(\rm CRT\) 合并即可。

Code

//18 = 9 + 9 = 18.
#include <iostream>
#include <cstdio>
typedef long long ll;
using namespace std;

const int MAXN = 35620;
const int MOD = 999911659;

int qpow(int a, int b, int p)
{
	int base = a, ans = 1;
	while (b)
	{
		if (b & 1)
		{
			ans = (ll)ans * base % p;
		}
		base = (ll)base * base % p;
		b >>= 1;
	}
	return ans;
}

int inv(int a, int p)
{
	return qpow(a, p - 2, p);
}

int fac[MAXN], inv_fac[MAXN];

void init(int p)
{
	fac[0] = 1;
	for (int i = 1; i < p; i++)
	{
		fac[i] = (ll)fac[i - 1] * i % p;
	}
	inv_fac[p - 1] = inv(fac[p - 1], p);
	for (int i = p - 2; i >= 0; i--)
	{
		inv_fac[i] = (ll)inv_fac[i + 1] * (i + 1) % p;
	}
}

int C(int n, int m, int p)
{
	if (m > n)
	{
		return 0;
	}
	return (ll)fac[n] * inv_fac[n - m] % p * inv_fac[m] % p;
}

int Lucas(int n, int m, int p)
{
	if (!m)
	{
		return 1;
	}
	return (ll)Lucas(n / p, m / p, p) * C(n % p, m % p, p) % p;
}

const int a[5] = {0, 2, 3, 4679, 35617};
int b[5];

int CRT()
{
	int m = MOD - 1, ans = 0;
	for (int i = 1; i <= 4; i++)
	{
		int mi = m / a[i];
		int Mi = inv(mi, a[i]);
		ans = (ans + (ll)b[i] * mi * Mi) % m;
	}
	return ans;
}

int main()
{
	int n, g;
	scanf("%d%d", &n, &g);
	if (!(g % MOD))
	{
		puts("0");
		return 0;
	}
	for (int i = 1; i <= 4; i++)
	{
		init(a[i]);
		for (int j = 1; j * j <= n; j++)
		{
			if (!(n % j))
			{
				b[i] = (b[i] + Lucas(n, j, a[i])) % a[i];
				if (j * j != n) //注意完全平方数
				{
					b[i] = (b[i] + Lucas(n, n / j, a[i])) % a[i];
				}
			}
		}
	}
	printf("%d\n", qpow(g, CRT(), MOD));
	return 0;
}