「SDOI2010」古代猪文（bzoj1951）

题目写了一大堆背景。

一句话题意就是求 $q^{\sum_{d|n}C_{n}^{d}} \mod 999911659$。

 

因为$n$是质数，只有当$q$是$n$的倍数时（此题数据范围原因，最多$q=n$），两数不互质，无法使用一些同余方程公式。这种特殊情况下可直接观察原式，发现答案就是$0$。

下面考虑的就都是$q≠n$的倍数的情况了。

 

根据$φ$的性质，由于这里$n(999911659)$是个质数，所以欧拉函数$φ(n)=n-1=999911658$；且当前考虑的是$q≠n$的情况，所以$q,n$必定互质，可以用欧拉定理。

先放此题结论：$q^{\sum_{d|n}C_{n}^{d}} \mod 999911659 \equiv q^{\sum_{d|n}C_{n}^{d}\mod 999911658} (\mod 999911659\space)$

但原欧拉定理不是这样的，这里推一下吧。

原欧拉定理：$q^{φ(M)} \equiv 1 (\mod M\space)$（设模数为$M$）

由于同余满足加、减、乘运算，两边同时翻$x(x为非负整数)$次方，得

$q^{φ(M)*x} \equiv 1 (\mod M\space)$

两边再同时乘以$a^y(y为非负整数)$，得

$q^{φ(M)*x+y} \equiv q^y (\mod M\space)$

设 $k=φ(M)*x+y$，则 $q^k \equiv q^y (\mod M\space)$

由于 $y=k \mod φ(M)$。

所以 $q^k \equiv q^{k\mod φ(M)} (\mod M\space)$

由于这题的值都是正整数，所以可以用非负的$x,y$得出的$k$表示任意非负整数，即可以表示$\sum_{d|n}C_{n}^{d}$这种非负整数。

因此上述结论对 $q^{\sum_{d|n}C_{n}^{d}} \mod 999911659$ 成立。

 

于是本题的关键就是求 $\sum_{d|n}C_{n}^{d}\mod 999911658$ 了。

但$999911658$不是质数，没法直接套同余公式。我们需要把它分解质因数。

通过另写一个暴力程序可知把这个数分解质因数得 $999911658=2*3*4679*35617$

分解出来的质数都非常小，设其中一个为$M$，我们直接$O(\sqrt n)$地枚举$n$的约数$d$，用 $lucas$ 求组合数$C_{n}^{d}\mod M$。

求和即可解出$\sum_{d|n}C_{n}^{d}\mod M$。

然后因为四个模数都是质数，所以可以用中国剩余定理合并答案。

$x\mod 2=a_1$

$x\mod 3=a_2$

$x\mod 4769=a_3$

$x\mod 35617=a_4$

即可得到 $\sum_{d|n}C_{n}^{d}\mod 999911658$ 的最小非负整数解$x$。再用快速幂求 $q^x$ 即可得到原问题的答案。

 

#include<bits/stdc++.h>
#define ll long long
using namespace std;
inline ll read(){
    ll x=0; bool f=1; char c=getchar();
    for(;!isdigit(c);c=getchar()) if(c=='-') f=0;
    for(; isdigit(c);c=getchar()) x=(x<<3)+(x<<1)+(c^'0');
    if(f) return x;
    return 0-x;
}
const ll a[4]={2,3,4679,35617};
ll n,g,p[4][35620],b[4],inv[4][35620];
ll ans,mod=999911659;
inline ll C(ll x,ll y,ll mod){
    if(x>y) return 0;
    return p[mod][y]*inv[mod][p[mod][y-x]*p[mod][x]%a[mod]]%a[mod];
}
ll lucas(ll x,ll y,ll mod){ //求组合数C(x,n)%a[mod] 
    if(x==0) return 1;
    return C(x%a[mod],y%a[mod],mod)*lucas(x/a[mod],y/a[mod],mod);
}
void exgcd(ll a,ll b,ll &x,ll &y){
    if(b==0){x=1,y=0; return;}
    exgcd(b,a%b,x,y);
    ll t=x; x=y, y=t-(a/b)*y;
}
ll Pow(ll a,ll b){
    ll ret=1;
    while(b>0){
        if(b&1) (ret*=a)%=mod;
        (a*=a)%=mod;
        b>>=1;
    }
    return ret;
}
int main(){
    n=read(),g=read();
    g%=mod;
    if(!g){printf("0\n"); return 0;}
    ll i,j;
    for(i=0;i^4;++i){
        for(j=p[i][0]=1;j<=a[i];++j) p[i][j]=p[i][j-1]*j%a[i];
    }
    for(i=0;i^4;++i){
        inv[i][0]=inv[i][1]=1;
        for(j=2;j^a[i];++j){
            inv[i][j]=-(a[i]/j)*inv[i][a[i]%j],
            inv[i][j]=(inv[i][j]%a[i]+a[i])%a[i];
            //prllf("%d %d %d %d %d",a[i],j,-(a[i]/j),inv[i][a[i]%j],inv[i][j]); system("pause");
        }
    }
    ll to=sqrt(n);
    for(i=1;i<=to;++i){
        for(j=0;j^4;++j){
            if(n%i==0){
                b[j]=(b[j]+lucas(i,n,j))%a[j];
                if(i*i!=n)
                    b[j]=(b[j]+lucas(n/i,n,j))%a[j];
            }
        }
    }
    --mod;
    ll x,y;
    for(i=0;i^4;++i){
        exgcd(mod/a[i],a[i],x,y);
        x=(x%a[i]+a[i])%a[i];
        ans=(ans+x*(mod/a[i])%mod*b[i])%mod;
    }
    ++mod;
    printf("%lld\n",Pow(g,ans));
    return 0;
}