『排列计数·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}\ast 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\ast 3\ast 4679\ast 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行哈哈哈哈哈