「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$ 即可得到原问题的答案。
1 #include<bits/stdc++.h> 2 #define ll long long 3 using namespace std; 4 inline ll read(){ 5 ll x=0; bool f=1; char c=getchar(); 6 for(;!isdigit(c);c=getchar()) if(c=='-') f=0; 7 for(; isdigit(c);c=getchar()) x=(x<<3)+(x<<1)+(c^'0'); 8 if(f) return x; 9 return 0-x; 10 } 11 const ll a[4]={2,3,4679,35617}; 12 ll n,g,p[4][35620],b[4],inv[4][35620]; 13 ll ans,mod=999911659; 14 inline ll C(ll x,ll y,ll mod){ 15 if(x>y) return 0; 16 return p[mod][y]*inv[mod][p[mod][y-x]*p[mod][x]%a[mod]]%a[mod]; 17 } 18 ll lucas(ll x,ll y,ll mod){ //求组合数C(x,n)%a[mod] 19 if(x==0) return 1; 20 return C(x%a[mod],y%a[mod],mod)*lucas(x/a[mod],y/a[mod],mod); 21 } 22 void exgcd(ll a,ll b,ll &x,ll &y){ 23 if(b==0){x=1,y=0; return;} 24 exgcd(b,a%b,x,y); 25 ll t=x; x=y, y=t-(a/b)*y; 26 } 27 ll Pow(ll a,ll b){ 28 ll ret=1; 29 while(b>0){ 30 if(b&1) (ret*=a)%=mod; 31 (a*=a)%=mod; 32 b>>=1; 33 } 34 return ret; 35 } 36 int main(){ 37 n=read(),g=read(); 38 g%=mod; 39 if(!g){printf("0\n"); return 0;} 40 ll i,j; 41 for(i=0;i^4;++i){ 42 for(j=p[i][0]=1;j<=a[i];++j) p[i][j]=p[i][j-1]*j%a[i]; 43 } 44 for(i=0;i^4;++i){ 45 inv[i][0]=inv[i][1]=1; 46 for(j=2;j^a[i];++j){ 47 inv[i][j]=-(a[i]/j)*inv[i][a[i]%j], 48 inv[i][j]=(inv[i][j]%a[i]+a[i])%a[i]; 49 //prllf("%d %d %d %d %d",a[i],j,-(a[i]/j),inv[i][a[i]%j],inv[i][j]); system("pause"); 50 } 51 } 52 ll to=sqrt(n); 53 for(i=1;i<=to;++i){ 54 for(j=0;j^4;++j){ 55 if(n%i==0){ 56 b[j]=(b[j]+lucas(i,n,j))%a[j]; 57 if(i*i!=n) 58 b[j]=(b[j]+lucas(n/i,n,j))%a[j]; 59 } 60 } 61 } 62 --mod; 63 ll x,y; 64 for(i=0;i^4;++i){ 65 exgcd(mod/a[i],a[i],x,y); 66 x=(x%a[i]+a[i])%a[i]; 67 ans=(ans+x*(mod/a[i])%mod*b[i])%mod; 68 } 69 ++mod; 70 printf("%lld\n",Pow(g,ans)); 71 return 0; 72 }
