BZOJ1951 [Sdoi2010]古代猪文

本文版权归ljh2000和博客园共有，欢迎转载，但须保留此声明，并给出原文链接，谢谢合作。

本文作者：ljh2000
作者博客：http://www.cnblogs.com/ljh2000-jump/
转载请注明出处，侵权必究，保留最终解释权！

题目链接：BZOJ1951

正解：中国剩余定理+费马小定理+$Lucas$定理

解题报告：

　　这道题可以说是数论+组合大综合题…

　　题目求的是$sum_{d|n} g^{C_n^d}$，模数$p$是质数（特判$g=p$），根据费马小定理，我们可以把幂上的数对$p-1$取模。

　　然后转化为求$sum_{d|n} C_n^d$ $mod$ $(p-1)$

 　  $n$高达$10^9$，而且模数不是个质数...

　　考虑用$Lucas$定理来实现大组合数取模：$C(n,m)=C(n/p,m/p)*C(n$ $mod$ $p,m$ $mod$ $p)$ $mod$ $p$

 　 然而模数不是质数，接下来就是这类题目的惯用做法了：把$p-1$质因数分解，对于每个质因子做一遍$Lucas$，最后用中国剩余定理合并起来得到最终$ans$。

　　概括起来就是：暴枚$n$的因数 ，把模数质因数分解之后对于每个质因子用$Lucas$定理实现大组合数取模，最后用$CRT$合并。

//It is made by ljh2000
//有志者，事竟成，破釜沉舟，百二秦关终属楚；苦心人，天不负，卧薪尝胆，三千越甲可吞吴。
#include <iostream>
#include <cstdlib>
#include <cstring>
#include <cstdio>
#include <cmath>
#include <algorithm>
#include <ctime>
#include <vector>
#include <queue>
#include <map>
#include <set>
#include <string>
#include <complex>
#include <bitset>
using namespace std;
typedef long long LL;
typedef long double LB;
const double pi = acos(-1);
const int MAXN = 50011;
const int mo = 999911659;
int prime[4]={2,3,4679,35617};
LL jie[4][MAXN],nj[4][MAXN],n,g,m[4],r[4],ans;
//费马小定理之后，转化为求sum {C_n^d} mod p-1

inline int getint(){
    int w=0,q=0; char c=getchar(); while((c<'0'||c>'9') && c!='-') c=getchar();
    if(c=='-') q=1,c=getchar(); while (c>='0'&&c<='9') w=w*10+c-'0',c=getchar(); return q?-w:w;
}

inline LL fast_pow(LL x,LL y,int mod){ LL r=1; while(y>0) { if(y&1) r*=x,r%=mod; x*=x; x%=mod; y>>=1; } return r; }
inline LL C(int n,int m,int k){ if(n<m) return 0; return jie[k][n]*nj[k][m]%prime[k]*nj[k][n-m]%prime[k]; }
inline void exgcd(LL a,LL b,LL &d,LL &x,LL &y){ 
	if(b==0) { d=a; x=1; y=0; return ; }
	exgcd(b,a%b,d,y,x);
	y-=a/b*x;
}

inline LL Lucas(int n,int m,int k){
	if(n<m) return 0; if(m==0) return 1;
	if(n<prime[k] && m<prime[k]) return C(n,m,k);
	return Lucas(n/prime[k],m/prime[k],k) * Lucas(n%prime[k],m%prime[k],k) %prime[k];
}

inline LL CRT(){
	LL M=m[0],R=r[0],d,x,y,z,bb;
	for(int i=1;i<4;i++) {
		exgcd(M,m[i],d,x,y); z=r[i]-R;
		bb=m[i]; bb/=d; z/=d;
		x*=z; x%=bb; x+=bb; x%=bb;
		R+=x*M;//!!!
		M=M*m[i]/d;
	}
	return R;
}

inline void solve(int x){
	for(int i=0;i<4;i++)
		r[i]+=Lucas(n,x,i),r[i]%=prime[i];
	//ans*=fast_pow(g,CRT(),mo); ans%=mo;可以最后一起CRT
}

inline void work(){
	n=getint(); g=getint(); if(g==mo) { printf("0"); return ; }//!!!
	g%=mo;//!!!
	for(int j=0;j<4;j++) {
		jie[j][0]=nj[j][0]=1; m[j]=prime[j];
		for(int i=1;i<prime[j];i++)
			jie[j][i]=jie[j][i-1]*i,jie[j][i]%=prime[j];
		nj[j][ prime[j]-1 ] = fast_pow(jie[j][ prime[j]-1 ],prime[j]-2,prime[j]);
		for(int i=prime[j]-2;i>=1;i--)
			nj[j][i]=nj[j][i+1]*(i+1),nj[j][i]%=prime[j];
	}

	for(int i=1;i*i<=n;i++) {
		if(n%i!=0) continue;//!!!!!!
		solve(i);
		if(i*i!=n) solve(n/i);
	}
	ans=fast_pow(g,CRT(),mo);
	printf("%lld",ans);
}

int main()
{
    work();
    return 0;
}
//有志者，事竟成，破釜沉舟，百二秦关终属楚；苦心人，天不负，卧薪尝胆，三千越甲可吞吴。