[Sdoi2010 古代猪文]师傅说的数论基础、谢天神牛出的题、BT的数论~、Lucas、中国剩余定理

【题目大意】
求G^(Sum(C(n,i))) % 999911659
其中i | n。
【算法分析】
我们来围观999911659这个数字。
首先他是一个素数。其次，它-1以后是一个square free number。
素数有什么性质呢？
那就是a^x = a^(x % phi(p)) (mod p)。其中p是素数。
然后phi(素数)=这个素数-1。
现在问题聚焦成求 Sum(C(n,i)) % 999911658。  其中 i | n
由于那个可爱的999911658是square free number，于是我们将他分解质因数。
发现只有4个素因子，分别是2、3、4679、35617
然后呢？
因为是分解质因数，所以它们都是素数（废话=.=）
我们单独对这4个素因子求Sum(C(n,i)) % p。
于是乎因为p是素数，所以我们来一个Lucas定理，就可以很方便的算出组合数。
当然，期间还要用到我以前说的“不知名定理”（实际上叫“欧拉定理”）。
（不知道的同学看这里：http://hi.baidu.com/edwardmj/blog/item/24c2e52c0133625c4fc2264a.html）
最后呢，我们得出4个同余方程：
X=A1 (mod 2)
X=A2 (mod 3)
X=A3 (mod 4679)
X=A4 (mod 35617)
好嘞，由于mod的都是素数，所以必然它们是两两互质的。
所以直接上中国剩余定理。
不知道的同学可以参照《算法艺术与信息学竞赛》 （黑书）的P231。
当然这里面涉及到的拓展欧几里得算法在黑书里也有，位于P219。
【其它】
这里面有个小技巧，就是由于算组合数时要用到那个阶乘，可以预处理出来。
由于使用了Lucas定理，所以每个阶乘都不会超过要mod的那个数字。
于是只用处理4个数组，它们的长度分别是2、3、4679、35617。
令人激动的是：1A。
更令人激动的是：
1 27533 edward_mj 632K 15MS G++ 2.04K 2010-06-02 21:14:09 2 27541 AekdyCoin 404K 30MS G++ 2.30K 2010-06-02 21:33:18 3 27538 AekdyCoin 456K 30MS G++ 2.78K 2010-06-02 21:27:11 4 27543 AekdyCoin 400K 45MS G++ 2.17K 2010-06-02 21:36:11 5 27546 AekdyCoin 396K 60MS G++ 2.14K 2010-06-02 21:41:21 6 26437 sonicmisora 412K 230MS G++ 1.92K 2010-05-28 10:56:15 7 26483 zxytim 384K 307MS G++ 2.50K 2010-05-28 14:05:33 8 26349 root 1344K 384MS G++ 2.67K 2010-05-27 20:19:52 9 27515 sevenk 348K 762MS G++ 3.55K 2010-06-02 20:48:10 10 26890 AekdyCoin 348K 777MS G++ 3.04K 2010-05-30 20:43:11 11 26886 AekdyCoin 352K 808MS G++ 2.90K 2010-05-30 20:32:14 12 26438 hyf042 844K 929MS G++ 2.69K 2010-05-28 11:05:21 13 26891 AekdyCoin 288K 933MS G++ 3.31K 2010-05-30 20:52:39 当然，师傅由于之前写疵，所以后面那堆是他在疵状态下写的。
前面的是在使用了那个小技巧之后的时间。。。
【CODE】
#include #include #include #include using namespace std;
typedef long long lld;
const int phi=999911658;
int mm=999911659;
const int N=200000;
struct zys{int L[N],Num[N],t;}pz,nz;
int n,G;
lld a[N];
lld JC[5][37777];

void Div(zys &p,int x){
p.t=0;
for (int i=2;i*i<=x;i++)
if (x%i==0){
p.L[++p.t]=i;
p.Num[p.t]=0;
while (x%i==0){
p.Num[p.t]++;
x/=i;
}
}
if (x!=1) {p.L[++p.t]=x; p.Num[p.t]=1;}
}

void makejc(){
int i,j;
for (i=1;i<=4;i++){
JC[i][0]=1;
for (j=1;j<=pz.L[i];j++)
JC[i][j]=JC[i][j-1]*j%pz.L[i];
}
}

lld Pow_mod(lld x,lld cf,int &mod){
if (cf==0) return 1;
lld res=Pow_mod(x,cf/2,mod);
res=res*res%mod;
if (cf%2) res=res*x%mod;
return res;
}

lld Normal_C(int m,int n,int i){
int p=pz.L[i];
return JC[i][m]*Pow_mod(JC[i][n]*JC[i][m-n]%p,p-2,p)%p;
}

lld C(int m,int n,int i){
if (!n) return 1;
int p=pz.L[i];
if (m%p return Normal_C(m%p,n%p,i)*C(m/p,n/p,i)%p;
}

void deal(int sum){
for (int i=1;i<=4;i++)
a[i]=(a[i]+C(n,sum,i))%pz.L[i];
}

void dfs(int dep,int sum){
if (dep>nz.t){
deal(sum);
return;
}
dfs(dep+1,sum);
for (int i=1;i<=nz.Num[dep];i++){
sum*=nz.L[dep];
dfs(dep+1,sum);
}
}

lld exgcd(lld a,lld b,lld &x,lld &y){
if (!b){
x=1; y=0;
return a;
}
lld res=exgcd(b,a%b,x,y),t;
t=x;
x=y;
y=t-a/b*y;
}

void solve(){
lld i,j,Mi,ans=0,p,q,M=phi;
for (i=1;i<=4;i++){
Mi=M/pz.L[i];
exgcd(Mi,pz.L[i],p,q);
ans=(ans+Mi*p*a[i])%M;
}
if (ans<0) ans+=M;
cout << Pow_mod(G,ans,mm) << "n";
}

int main(){
cin >> n >> G;
G%=mm;
if (G==0) {printf("0n"); return 0;}
Div(pz,phi);
Div(nz,n);
makejc();
dfs(1,1);
solve();
return 0;
}