BZOJ 1951: [Sdoi2010]古代猪文

好困啊刷点水题醒醒脑

题目要求G^(sigma(k|n,C(n,k))%P的值

根据费马小定理我们可以把指数对P-1求余

但是问题是P-1不是质数

所以将P-1分解为4个质数的乘积，分别求解之后用CRT合并（终于知道CRT合并是啥了）

指数求出来之后快速幂搞一搞就可以了

#include<cstdio>
#include<iostream>
#include<cstring>
#include<cmath>
#include<queue>
#include<vector>
#include<algorithm>
#include<map>
#include<set>
#include<stack>
#include<cstdlib>
#include<ctime>
#define rep(i,l,r) for(int i=l;i<=r;i++)
#define per(i,r,l) for(int i=r;i>=l;i--)
#define mmt(a,v) memset(a,v,sizeof(a))
#define tra(i,u) for(int i=head[u];i;i=e[i].next)
using namespace std;
typedef long long ll;
const int P=999911659;
ll qmul(ll a,ll b,ll p){
	ll ans=1;
	for(;b;b>>=1,a=a*a%p)if(b&1)ans=ans*a%p;
	return ans;
}
int t[]={2,3,4679,35617};
ll f[4][40000],g[4][40000],inv[4][40000];
ll c(int n,int m,int i){
	if(n<m)return 0;
	if(n==m)return 1;
	return f[i][n]*g[i][m]%t[i]*g[i][n-m]%t[i];
}
ll lucas(ll n,ll m,int i){
	if(!m)return 1;
	return c(n%t[i],m%t[i],i)*lucas(n/t[i],m/t[i],i)%t[i];
}
void exgcd(ll a,ll b,ll &x,ll &y){
	if(b){exgcd(b,a%b,y,x);y-=x*(a/b);}
	else{x=1;y=0;}
}
ll solve(ll a,int i){
	ll m=(P-1)/t[i],x,y;
	exgcd(m,t[i],x,y);
	x=(x+t[i])%t[i];
	return a*m*x%(P-1);
}
ll ans[4];
int main(){
	//freopen("a.in","r",stdin);
	rep(i,0,3){
		f[i][0]=g[i][1]=inv[i][1]=g[i][0]=1;
		rep(j,1,t[i])f[i][j]=f[i][j-1]*j%t[i];
		rep(j,2,t[i]-1)inv[i][j]=(t[i]-(t[i]/j)*inv[i][t[i]%j]%t[i]);
		rep(j,2,t[i]-1)g[i][j]=g[i][j-1]*inv[i][j]%t[i];
	}
	int n,g;scanf("%d%d",&n,&g);
	if(g==P){puts("0");return 0;}
	g%=P;
	rep(i,1,n){
		if(i*i>n)break;
		if(n%i==0){
			int tmp=n/i;
			rep(j,0,3){
				if(tmp!=i)ans[j]=(ans[j]+lucas(n,i,j))%t[j];
				ans[j]=(ans[j]+lucas(n,tmp,j))%t[j];
			}
		}
	}
	int res=0;
	rep(i,0,3)res=(res+solve(ans[i],i))%(P-1);
	printf("%lld\n",qmul(g,res,P));
	return 0;
}