1951: [Sdoi2010]古文字猪

1951: [Sdoi2010]古代猪文

题目：

一道非常好的组合数学题。！！。题目非常长。只是就以下几段话实用。

iPig认为仅仅要符合文献，每一种能整除N的k都是有可能的。他打算考虑到全部可能的k。显然当k等于某个定值时，该朝的猪文文字个数为N / k。然而从N个文字中保留下N / k个的情况也是相当多的。iPig估计，假设全部可能的k的全部情况数加起来为P的话。那么他研究古代文字的代价将会是G的P次方。如今他想知道猪王国研究古代文字的代价是多少。

因为iPig认为这个数字可能是天文数字。所以你仅仅须要告诉他答案除以999911659的余数就能够了。

就是要求出P = C(N,i|N),然后解出G^P % 999911659就是最后的答案。

我们能够非常easy的看出问题是由两部分组成的，求解G^P和P。

我们能够easy的得出G^P用高速幂能够easy的解决，如今的问题是怎样高速的求出P呢？这是一个难题。我们在继续的分解，发现P的组成了没有。对！就是组合公式C(N,i|N)！而我们又能够知道组合公式假设数据较小的时候能够用杨辉三角的暴力公式得到。而假设是大组合数的话就要用到Lucas定理。（AekdyCoin空间给出了这方面的具体说明）而这题的组合数显然是大组合数。

而直接求就算是long long 也要益处啊！

！怎么办？我们能够想到大组合数的通常处理方法。就是取模。

模?哪来的模？而此时我们又有费马小定理能够知道当P是素数的时候:

G^P % MOD = G ^ (P % (MOD - 1)) % MOD

所以，此时MOD - 1 = 999911658 是一个偶数，如今就不符合Lucas定理了。

所以，此时我们要对这个偶数进行质数分解。999911658 = 2 * 3 * 4679 * 35617 可是，如今问题又来了。你把模给分解了，那怎样求出结果呢？

我们在好好研究一下。你发现了什么没有？对!就是有例如以下的等式关系：

x = a1 % m1

x = a2 % m2

x = a3 % m3

对！

就是中国剩余定理！

如今。经过了一步一步的分解最终把问题攻克了。你该懂了吧！

一道神题啊！

！

/*
    算法：模非素数的组合数
    题目：就是要求出P = C(N,i|N),然后解出G^P % 999911659
    1、对模质因子分解
    2、中国剩余定理
    1 <= G <= 10^9，1 <= N <= 10^9

*/
#include <iostream>
#include <algorithm>
#include <cstdio>
#include <cstring>
#include <cmath>
using namespace std;

typedef long long LL;
const int phi = 999911658;
const int MOD = 999911659;
const int MAXN = 200000;

struct Prim{
   int cnt,prim[MAXN],num[MAXN];
}pz,nz;

int n,G;
LL a[MAXN];
LL fact[5][37777];

//分解质因子
void Div(Prim& p,int x){
    p.cnt = 0;   //初始化
    int k = sqrt(x + 0.5);
    for(int i = 2;i <= k;++i){
        if(0 == x % i){
            p.num[++p.cnt] = 0;
            p.prim[p.cnt] = i;
            while(0 == x % i){
                ++p.num[p.cnt];
                x /= i;
            }
        }
    }

    if(x > 1){
        p.num[++p.cnt] = 1;
        p.prim[p.cnt] = x;
    }
}

//预处理阶乘
void init(){
   int i,j;
   for(i = 1;i <= pz.cnt;++i){
       fact[i][0] = 1;
       for(j = 1;j <= pz.prim[i];++j)
          fact[i][j] = fact[i][j-1] * j % pz.prim[i];
   }
}

//高速幂取模
LL powmod(LL a,LL b,const int& mod){
   LL res = 1;
   a %= mod;

   while(b > 0){
       if(b & 1) res = res * a % mod;
       a = a * a % mod;
       b >>= 1;
   }

   return res;
}

LL Norma_C(int n,int m,int i){
    int p = pz.prim[i];
    return fact[i][n] * powmod(fact[i][m]*fact[i][n-m]%p,p-2,p) % p;
}

//Lucas
LL Lucas(int n,int m,int i){
    if(!m) return 1;
    int p = pz.prim[i];
    if(n%p < m%p) return 0;
    LL tmp = Norma_C(n%p,m%p,i);
    return tmp * Lucas(n / p,m / p,i) % p;
}

//计算C(N,i)
void deal(int sum){
    for(int i = 1;i <= pz.cnt;++i){
        a[i] = (a[i] + Lucas(n,sum,i)) % pz.prim[i];
    }
}

//查找 i|N
void dfs(int dep,int sum){
     if(dep > nz.cnt){
        deal(sum);
        return;
     }

     dfs(dep+1,sum); //不要该数
     for(int i = 1;i <= nz.num[dep];++i){ //要该数的情况
        sum *= nz.prim[dep];
        dfs(dep+1,sum);
     }
}

//扩展欧几里得
void extgcd(LL a,LL b,LL& d,LL& x,LL& y){
    if(!b){ d = a; x = 1; y = 0; }
    else { extgcd(b,a % b,d,y,x); y -= x * (a / b); }
}

//中国剩余定理
LL china(){
   LL M = phi,d,y,x = 0;
   for(int i = 1;i <= pz.cnt;++i){
      LL w = M / pz.prim[i];
      extgcd(pz.prim[i],w,d,d,y);
      x = (x + y * w * a[i]) % M;
   }
   x = (x + M) % M;
   return x;
}

int main()
{
    while(~scanf("%d%d",&n,&G)){
        memset(a,0,sizeof(a));
        G %= MOD;
        if(!G){
            puts("0");
            continue;
        }
        
        ////////////////////////////////
        
        Div(pz,phi);    //对模进行质因子分解        
        Div(nz,n);      //对数进行质因子分解
        init();         //预处理阶乘
        dfs(1,1);       
        LL ans = china();  //中国剩余定理
        
        //////////////////////////////////
        
        printf("%lld
",powmod((LL)G,ans,MOD));   //G^P
    }
    return 0;
}