题解 古代猪文

题意

求出以下式子的值：

\[\large \begin{aligned} g^{\sum_{d\mid n} {n\choose d}} \mod 999911659 \end{aligned} \]

首先观察到 \(999911659\) 是一个质数，且考虑到指数过大，考虑欧拉定理：

\[\large \begin{aligned} g^{999911659-1}\equiv1(mod\:\;999911659) \end{aligned} \]

于是原式变为：

\[\large \begin{aligned} g^{\sum_{d\mid n} {n\choose d}\mod 999911658} \mod 999911659 \end{aligned} \]

现在指数已经可以接受，所以考虑计算如下式子：

\[\large \begin{aligned} \sum_{d\mid n} {n\choose d}\mod 999911658 \end{aligned} \]

lucas 定理？不行！首先因为 \(999911658\) 不是质数，其次是 \(999911658\) 这个数太大了。有一种巧妙的做法。

我们知道对于模数两两互质的同余方程组，有唯一解，所以考虑分解质因数然后用 CRT 得出解，解即为贡献，这样模数为质数且较小所以可以用 lucas 定理，单次 lucas 计算的时间复杂度是 \(O(\log p\cdot \log n)\) 的。

\[\large \begin{aligned} 999911658=2 \times 3 \times 4679 \times 35617 \end{aligned} \]

即解同余式：

\[\large \begin{aligned} x &\equiv {n \choose d} (mod\;\: 2)\\ x &\equiv {n \choose d} (mod\;\: 3)\\ x &\equiv {n \choose d} (mod\;\: 4679)\\ x &\equiv {n \choose d} (mod\;\: 35617) \end{aligned} \]

\(d\) 枚举即可，求组合数用 lucas 定理就行。

我要用 EXCRT！

注意当指数为 \(0\) 时不需要研究。

#include <bits/stdc++.h>
#define int long long
using namespace std;
const int N=10,M=35620,P=999911659,PP=999911658;
int a[N],m[N];
int fac[4][M];
unordered_map<int,int> num{{2,0},{3,1},{4679,2},{35617,3}};
unordered_map<int,int> pos{{0,2},{1,3},{2,4679},{3,35617}};
int qmi(int a,int b,int p) {
	int res=1%p;
	while(b) {
		if(b&1) res=res*a%p;
		a=a*a%p;
		b>>=1;
	}
	return res;
}
int exgcd(int a,int b,int &x,int &y) {
	if(!b) {
		x=1,y=0;
		return a;
	}
	int d=exgcd(b,a%b,y,x);
	y-=a/b*x;
	return d;
}
int excrt(int n) {
	int a1,a2,m1,m2;
	a1=a[1],m1=m[1];
	for(int i=2;i<=n;i++) {
		a2=a[i],m2=m[i];
		int k1,k2;
		int d=exgcd(m1,m2,k1,k2);
		k1*=(a2-a1)/d;
		k1=(k1%(m2/d)+(m2/d))%(m2/d);
		a1+=m1*k1;
		m1*=m2/d;
	}
	return (a1%m1+m1)%m1;
}
int C(int x,int y,int p) {
	if(x<y) return 0;
	return fac[num[p]][x]*qmi(fac[num[p]][y],p-2,p)%p*qmi(fac[num[p]][x-y],p-2,p)%p;
}
int lucas(int x,int y,int p) {
	if(!y) return 1;
	return C(x%p,y%p,p)*lucas(x/p,y/p,p)%p;
}
int n,g,ans;
int solve(int x) {
	int top=0;
	for(int i=0;i<4;i++) {
		int tmp=lucas(n,x,pos[i]);
		//if(tmp==0) continue;
		a[++top]=tmp;
		m[top]=pos[i];
	}
	return excrt(top)%PP;
}
signed main() {
	cin>>n>>g;
	for(int i=0;i<4;i++) {
		fac[i][0]=1;
		for(int j=1;j<=pos[i];j++) fac[i][j]=fac[i][j-1]*j%pos[i];
	}
	for(int i=1;i<=n/i;i++)
		if(n%i==0) {
			if(i*i!=n) ans=(ans+solve(n/i))%PP;
			ans=(ans+solve(i))%PP;
		}
	if(!ans) puts("0");
	else cout<<qmi(g,ans,P)<<endl;
	return 0;
}