bzoj1951 [Sdoi2010]古代猪文 ( 欧拉定理+CRT+lucas )

bzoj1951 [Sdoi2010]古代猪文

原题地址：http://www.lydsy.com/JudgeOnline/problem.php?id=1951

题意：
一个朝代流传的猪文文字恰好为N的k分之一，其中k是N的一个正约数（可以是1和N）。不过具体是哪k分之一，以及k是多少，由于历史过于久远，已经无从考证了。考虑到所有可能的k。显然当k等于某个定值时，该朝的猪文文字个数为N / k。然而从N个文字中保留下N / k个的情况也是相当多的。如果所有可能的k的所有情况数加起来为P的话，那么他研究古代文字的代价将会是G的P次方。 现在他想知道研究古代文字的代价除以999911659的余数是多少。
即
给出N，G，求：

$G^{\sum_{k|N}C_n^k}$

数据范围
1 <= G <= 1000000000，1 <= N <= 1000000000。

题解：
$P=\sum_{k|N}C_n^k$
求$G^P (mod\ \ 999911659)$

思路比较简单，首先求G^P显然不能直接算，
由于欧拉定理
$a^x \equiv a^{x \ mod \ φ(m) + φ(m) } (mod \ m)$
$x>φ(m)$
因此转化为求：
$P=\sum_{k|N}C_n^k\ \ \ mod \ φ(m)$
$G^P (mod\ \ 999911659)$

由于φ(m)=999911658，并不是质数，所以不能直接用lucas，
因为 999911658=2*3*4679*35617，每个质因子的质数都是1，转化为分别求$\sum_{k|N}C_n^k$模四个质因子的值，得到四个模线性方程，再用CRT合并起来。

注意：
1.线性求逆元时求的那个inv[N]是N=mod-1，不然之上的fac都是0，没法计算逆元。
2.欧拉定理成立条件 gcd(G,mod)=1，因此当G=999911659特判答案为0。
欧拉定理的扩展$a^x \equiv a^{x \ mod \ φ(m) + φ(m) } (mod \ m)$并不要求a，m互质。

---

$a^x \equiv a^{x \ mod \ φ(m) + φ(m) } (mod \ m)$

欧拉定理的扩展并不要求a，m互质，
证明见Mercer的博客

代码：

#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cstring>
using namespace std;
const int mod=999911659;
const int N=40000;
int n,g,m[5]={0,2,3,4679,35617},fac[N+10][6],inv[N+10][6],a[10]; 
//999911658=2*3*4679*35617
int modpow(int A,int B,int Mod)
{
    int ans=1; int base=A;
    for(;B;B>>=1)
    {
        if(B&1) ans=(1LL*ans*base)%Mod;
        base=(1LL*base*base)%Mod;
    }
    return ans;
}

void init()
{
    for(int i=1;i<=4;i++)
    fac[0][i]=inv[0][i]=inv[m[i]][i]=1;

    for(int j=1;j<=4;j++)
    for(int i=1;i<=N;i++)
    fac[i][j]=(1LL*fac[i-1][j]*i)%m[j]; 

    for(int i=1;i<=4;i++)
    {
        inv[m[i]-1][i]=modpow(fac[m[i]-1][i],m[i]-2,m[i]);
        for(int j=m[i]-2;j>=1;j--)
        inv[j][i]=(1LL*inv[j+1][i]*(j+1))%m[i];
    }           
}
int comb(int A,int B,int i)
{
    if(A<B||A<0||B<0) return 0;
    int iv=(1LL*inv[B][i]*inv[A-B][i])%m[i];
    int ans=(1LL*fac[A][i]*iv)%m[i];
    return ans;
}
int lucas(int A,int B,int i)
{
    if(A<B) return 0;
    if(B==0) return 1;
    return (1LL*lucas(A/m[i],B/m[i],i)*comb(A%m[i],B%m[i],i))%m[i];
}
void solve(int x)
{
    for(int i=1;i<=4;i++)
    a[i]=(a[i]+lucas(n,x,i)%m[i])%m[i]; 
}
int crt()
{
    int M=1;
    for(int i=1;i<=4;i++) M=M*m[i];
    int x=0;
    for(int i=1;i<=4;i++)
    {
        int Mi=M/m[i];
        int Ri=modpow(Mi,m[i]-2,m[i]);
        int ret=(1LL*((1LL*Mi*Ri)%M)*a[i]%M)%M; 
        x=(x+ret)%M;
    }
    return x;
}
int main()
{
    init(); 
    scanf("%d%d",&n,&g);
    if(g==mod) 
    {
        printf("0\n"); return 0;
    }
    for(int i=1;i*i<=n;i++)
    {if(n%i==0) {solve(i); if(i!=n/i) solve(n/i);}} 
    int P=crt();
    printf("%d\n",modpow(g,P,mod));
    return 0;
}