【洛谷P2480】古代猪文

题目

题目链接：https://www.luogu.com.cn/problem/P2480
猪王国的文明源远流长，博大精深。
iPig 在大肥猪学校图书馆中查阅资料，得知远古时期猪文文字总个数为 \(n\)。当然，一种语言如果字数很多，字典也相应会很大。当时的猪王国国王考虑到如果修一本字典，规模有可能远远超过康熙字典，花费的猪力、物力将难以估量。故考虑再三没有进行这一项劳猪伤财之举。当然，猪王国的文字后来随着历史变迁逐渐进行了简化，去掉了一些不常用的字。
iPig 打算研究古时某个朝代的猪文文字。根据相关文献记载，那个朝代流传的猪文文字恰好为远古时期的 \(1/k\)，其中 \(k\) 是 \(n\) 的一个正约数（可以是 \(1\) 或 \(n\)）。不过具体是哪 \(1/k\)，以及 \(k\) 是多少，由于历史过于久远，已经无从考证了。
iPig 觉得只要符合文献，每一种 \(k|n\) 都是有可能的。他打算考虑到所有可能的 \(k\)。显然当 \(k\) 等于某个定值时，该朝的猪文文字个数为 \(n/k\)。然而从 \(n\) 个文字中保留下 \(n/k\) 个的情况也是相当多的。iPig 预计，如果所有可能的 \(k\) 的所有情况数加起来为 \(p\) 的话，那么他研究古代文字的代价将会是 \(g^p\)。
现在他想知道猪王国研究古代文字的代价是多少。由于 iPig 觉得这个数字可能是天文数字，所以你只需要告诉他答案除以 \(999911659\) 的余数就可以了。

思路：

一句话题意：给出\(n,g\)，求\(g^{\sum_{d|n}C^d_n}\ mod\ 999911659\)。
根据欧拉定理得

\[g^{\sum_{d|n}C^d_n}\ mod\ 999911659\ =\ g^{\sum_{d|n}C^d_n\ mod\ 999911658}\ mod\ 999911659 \]

筛出\(n\)的每一个因子\(d[]\)，考虑对于其中一个因子\(d[i]\)如何求\(C^{d[i]}_n\ mod\ 999911658\)
由于\(999911658=2\times 3\times 4679\times 35617\)，考虑到质因子较小且个数较少，我们可以用\(Lucas\)先求出\(a_1=C_n^{d[i]}\ mod\ 2,a_2=C_n^{d[i]}\ mod\ 3，a_3=C_n^{d[i]}mod\ 4679，a_4=C_n^{d[i]}\ mod\ 35617\)，然后解同余方程组

\[\left\{\begin{matrix}x\equiv a_1(mod\ 2)\\ x\equiv a_2(mod\ 3)\\ x\equiv a_3(mod\ 4679)\\ x\equiv a_4(mod\ 35617)\end{matrix}\right. \]

对于每一个因数\(d[i]\)我们都可以解除一个同于方程组的解\(x_i\)。最后答案就是\(g^{\sum_{d[i]}x_i}\)

代码

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
typedef long long ll;

const int N=80010,MOD=999911659;
const int v[]={0,2,3,4679,35617};
int n,m,g,s,ans,d[N];
ll fac[N];

ll power(ll x,ll k,ll p)
{
	ll ans=1;
	for (;k;k>>=1,x=x*x%p)
		if (k&1) ans=ans*x%p;
	return ans;
}

ll C(int n,int m,int p)
{
	if (n<m) return 0;
	ll niv1=power(fac[m],p-2,p),niv2=power(fac[n-m],p-2,p);
	return fac[n]*niv1%p*niv2%p;
}

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

ll crt(int j)
{
	ll M=(MOD-1)/v[j];
	return 1LL*s*M%(MOD-1)*power(M,v[j]-2,v[j])%(MOD-1);
}

int main()
{
	scanf("%d%d",&n,&g);
	if (g==MOD) return !printf("0");
	for (int i=1;i*i<=n;i++)
		if (!(n%i))
		{
			d[++m]=i;
			if (n/i!=i) d[++m]=n/i;
		}
	for (int j=1;j<=4;j++)
	{
		fac[0]=1;
		for (int i=1;i<v[j];i++)
			fac[i]=fac[i-1]*i%v[j];
		s=0;
		for (int i=1;i<=m;i++)
			s=(s+lucas(n,d[i],v[j]))%v[j];
		ans=(ans+crt(j))%(MOD-1);
	}
	printf("%lld",power(g,ans,MOD));
	return 0;
}