[SDOI2010] 古代猪文

Description

猪王国的文明源远流长，博大精深。

iPig在大肥猪学校图书馆中查阅资料，得知远古时期猪文文字总个数为N。当然，一种语言如果字数很多，字典也相应会很大。当时的猪王国国王考虑到如果修一本字典，规模有可能远远超过康熙字典，花费的猪力、物力将难以估量。故考虑再三没有进行这一项劳猪伤财之举。当然，猪王国的文字后来随着历史变迁逐渐进行了简化，去掉了一些不常用的字。

iPig打算研究古时某个朝代的猪文文字。根据相关文献记载，那个朝代流传的猪文文字恰好为远古时期的k分之一，其中k是N的一个正约数（可以是1和N）。不过具体是哪k分之一，以及k是多少，由于历史过于久远，已经无从考证了。

iPig觉得只要符合文献，每一种能整除N的k都是有可能的。他打算考虑到所有可能的k。显然当k等于某个定值时，该朝的猪文文字个数为N / k。然而从N个文字中保留下N / k个的情况也是相当多的。iPig预计，如果所有可能的k的所有情况数加起来为P的话，那么他研究古代文字的代价将会是G的P次方。

现在他想知道猪王国研究古代文字的代价是多少。由于iPig觉得这个数字可能是天文数字，所以你只需要告诉他答案除以999911659的余数就可以了。

Input：

有且仅有一行：两个数N、G，用一个空格分开。

Output

有且仅有一行：一个数，表示答案除以999911659的余数。

Range

10%的数据中，$1 <= N <= 50$；

20%的数据中，$1 <= N <= 1000$；

40%的数据中，$1 <= N <= 100000$；

100%的数据中，$1 <= G <= 1000000000，1 <= N <= 1000000000$。

Solution

先一句话题意：求 $G^{\Sigma _{d \mid n} C_n^d} \; \% 999911659$ 的值。

若 G =999911659，则上式结果为 0。否则，因为999911659是质数，所以 G，n 互质。由欧拉定理推论得：

$$G^{\Sigma _{d \mid n} C_n^d} \equiv G^{\sum _{d \mid n} C_n^d \; \% \; 999911658} \quad (mod\;999911659)$$

尝试分解质因数，可以发现 $999911658 = 2 \times 3 \times 4679 \times 35617$。因为所以质因子的指数都为 1，所以我们枚举 $n$ 的约数 $d$，然后运用 $Lucas$ 求出组合数 $C_n^d$ ，分别计算出 $\Sigma _{d \mid n} C_n^d$ 对 2,3,4679,35617 四个质数取模的结果。然后中国剩余定理合并即可。

快速幂非递归版本 1A，递归版本无限 wa。玄学错误=.=

Code

 

#include<cstdio>
#define mod 999911659
#define int long long

int n,g;
int m[6],r[6];
int factorcnt;
int inv[100005];
int fac[100005];
int factor[10005];

void init_factor(){
    for(int i=1;i*i<=n;i++){
        if(n%i) continue;
        factor[++factorcnt]=i;
        if(i*i!=n) factor[++factorcnt]=n/i;
    }
}

int lucas(int x,int y,int p){
    if(x<y) return 0;
    if(x>=p) return (lucas(x%p,y%p,p)*lucas(x/p,y/p,p))%p;
    return (fac[x]*inv[fac[y]]*inv[fac[x-y]])%p;
}

int calc(int x){
    fac[0]=fac[1]=1;
    inv[0]=inv[1]=1;
    for(int i=2;i<=m[x];i++) fac[i]=(fac[i-1]*i)%m[x];
    for(int i=2;i<=m[x];i++) inv[i]=((m[x]-m[x]/i)*inv[m[x]%i])%m[x];
    for(int i=1;i<=factorcnt;i++)
        (r[x]+=lucas(n,factor[i],m[x]))%=m[x];
}

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

int innv(int a,int b){
    int x,y;
    exgcd(a,b,x,y);
    return (x%b+b)%b;
}

int CRT(){
    int M=1,ans=0;
    for(int i=1;i<=4;i++) M*=m[i];
    for(int i=1;i<=4;i++)
        (ans+=(M/m[i])*innv(M/m[i],m[i])*r[i])%=M;
    return (ans%M+M)%M;
}

int ksm(int x,int y){
    /*if(y==0) return 1;
    if(y==1) return x%mod;
    int c=ksm(x,y>>1);
    if(y&1) return (x*c*c)%mod;
    return (c*c)%mod;*/
    int r=1;
    for(;y;y>>=1){
        if(y&1) (r*=x)%=mod;
        (x*=x)%=mod;
    }
    return r;
}

signed main(){
    scanf("%lld%lld",&n,&g);
    if(g==mod){
        puts("0");
        return 0;
    }
    m[1]=2; m[2]=3; m[3]=4679; m[4]=35617;
    init_factor();
    for(int i=1;i<=4;i++) calc(i);
    printf("%lld\n",ksm(g,CRT())%mod);
    return 0;
}