1951: [Sdoi2010]古代猪文

                                                       1951: [Sdoi2010]古代猪文

 

链接：Click Here~

题目：

   一道很好的组合数学题！！！！题目很长，不过就下面几段话有用。

  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就是最后的答案。

   我们可以很容易的看出问题是由两部分组成的，求解G^P和P。

   我们可以容易的得出G^P用快速幂可以容易的解决，现在的问题是如何快速的求出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\n",powmod((LL)G,ans,MOD));   //G^P
    }
    return 0;
}