『排列计数·Lucas定理』「SDOI2010」古代猪文

题目大意

题目很长，这里就不说了吧。粘个链接：古代猪文
题目大意就是：给定整数n,q，计算$q^{{\sum }_{d|n}{C}_n^d}mod999911659$.

题解

这道题的$n$很大；如果暴力枚举每一个$d$，再直接求解组合数的话，由于$n$和$d$很大，在计算组合数的时候会浪费很大的时间，因此我们需要对这一个算法进行优化。

在讨论本题之前我们先来讨论一下$Lucas$定理以及其实现过程。

• 对于任意质数P，有:$$C_n^m\equiv C_{n/P}^{m/P}?C_{n\%P}^{m\%P},(modp)$$
• 证明过程很复杂，我们在这里不详细讨论。我们主要考虑如何运用到组合数的计算内。
• 对于利用$Lucas$定理求解组合数，我们可以使用递归来进行求解，以求解$C_n^m$即：
  当$n=m$时返回1；若$n$和$m$都同时小于模数p，我们直接暴力计算组合数即可;否则我们就必须使用$Lucas$定理来进行递归实现。
• 友情提醒：必须要小于质数$p$才能提前计算而不能小于某一个预先处理好的大数字；因为当这一个大数字超过$p$时就会出现大量模$p$等于$0$的情况，那么你就会发现组合数计算出来以后全部都是0.但是一般情况下因为p比较大所以不存在任何问题；基础这道题的特殊性，模数比较小，组合数的计算上就会出现问题了。
• 代码如下：你可以将$p[k]$理解一个取模的质数；$power$是一个快速幂函数，用于求解乘法逆元。

LL C(LL n,LL m,LL k)
{
    if (n == m) return 1;if (n < m) return 0;if (n<p[k] && m<p[k]){
    LL sum1 = jc[n][k];LL sum2 = jc[n-m][k]*jc[m][k]%p[k];return sum1*power(sum2,p[k]-2,p[k])%p[k];}return C(n/p[k],m/p[k],k)*C(n%p[k],m%p[k],k)%p[k];
}

OK，说完卢卡斯以后回归本题。

首先我们需要特判一种特殊的情况，当q是999911659=P的倍数时，答案为0.

若不是上面的情况，我们可以根据欧拉定理推论:$a^b\equiv a^{b\%phi(p)}(modp)$,得到：
$$q^{{\sum }_{d|n}{C}_n^d}\equiv q^{{\sum }_{d|n}{C}_n^d\%999911658}$$

因此我们的目标就是要求解这一串数的值：${\sum }_{d|n}{C}_n^d\%999911658$

然后我们对模数$P$分解质因数，得到：$999911659=2?3?4679?35617$。可以分别计算。
然后这一题就非常好玩了，我们可以暴力求$n$的每一个因数，得到对于的$C_n^d$对四个质数取模的值，分别累加在$a_0$,$a_1$,$a_2$,$a_3$上。即$a_i={\sum }_{d|n}{C}_n^d\%p_i$.为什么可以累加？因为最后是求和的，因此余数可以直接累加。

最后，你可以得到一个方程组：
$$?\begin{array}{l}xmod2=a_1\\ xmod3=a_2\\ xmod4769=a_3\\ xmod35617=a_4\end{array}$$
显然这个东西直接套用中国剩余定理的结论即可；一通线性同余方程exgcd一下就好了。

然后我们就可以求解出一组$x_0$,则通解$x$可以表示为：$x_0\equiv x(mod999911658)$.我们只要找到最小的解即可。

为了加快求组合数的速度，我们要预处理出阶乘；乘法逆元求不求都无所谓。然后就十分愉快的得到代码了：

#include <bits/stdc++.h>
typedef long long LL;
using namespace std;const LL P = 999911659;
const LL p[4] = {
    2,3,4679,35617};
LL n,q;
LL a[4];
LL jc[100000][4];void init(void)
{
    jc[0][0] = jc[0][1] = jc[0][2] = jc[0][3] = 1;for (LL i=1;i<=5e4;++i) {
    jc[i][0] = jc[i-1][0]*i%p[0];jc[i][1] = jc[i-1][1]*i%p[1];jc[i][2] = jc[i-1][2]*i%p[2];jc[i][3] = jc[i-1][3]*i%p[3];}return;
}LL power(LL a,LL b,LL c) {
    LL res = 1;while (b>0){
    if (b%2 == 1) res = res*a%c;a = a*a%c;b /= 2;}return res;
}LL C(LL n,LL m,LL k)
{
    if (n == m) return 1;if (n < m) return 0;if (n<p[k] && m<p[k]){
    LL sum1 = jc[n][k];LL sum2 = jc[n-m][k]*jc[m][k]%p[k];return sum1*power(sum2,p[k]-2,p[k])%p[k];}return C(n/p[k],m/p[k],k)*C(n%p[k],m%p[k],k)%p[k];
}void work(void)
{
    for (LL i=0;i<4;++i){
    for (LL j=1;j*j<=n;++j) if (n%j == 0){
    if (j*j != n) a[i] = (a[i]+C(n,j,i))%p[i];a[i] = (a[i]+C(n,n/j,i))%p[i];}}
}LL x,y;
LL exgcd(LL a,LL b,LL c) 
{
     if (b == 0) {
     x = c/a; y = 0; return a; } LL t = exgcd(b,a%b,c); LL X = x,Y = y; x = Y; y = X-a/b*Y; return t; 
}void china(void)
{
    LL sum = 0;LL mm = 999911658;LL M[4] = {
    };LL r[4] = {
    };for (LL i=0;i<4;++i) {
    M[i] = mm/p[i];exgcd(M[i],p[i],1);r[i] = x;sum = sum+r[i]*M[i]*a[i];}sum = (sum%mm+mm)%mm;cout<<power(q,sum,P);return;
}int main(void)
{
    freopen("test.in","r",stdin);freopen("test.out","w",stdout);init();cin>>n>>q;if (q == P) return puts("0"),0;work();china();return 0;
} 

一道数学题硬是被我写到了100行哈哈哈哈哈