扩展卢卡斯学习笔记

前几天做了一道计数题，本来挺水的，非得出成模数不是质数，于是我就来学扩展卢卡斯了。
这东西感觉还不难，比较好理解。


我们要求的就是(C_{n} ^ {m} \% p)。因为(p)不一定是质数，所以可以把(p)质因数分解，后求出(C_{n} ^ {m} \% p _ {i} ^ {k})的解，这样用中国剩余定理解线性同余方程组就行了。
所以，


1.求(C_{n} ^ {m} \% p ^ {k})
我们把组合数写成阶乘的形式，然后提出里面的(p)，这样(m!, (n - m)!)就和(p ^ k)互质，就可以用扩展欧几里得就逆元了。也就是变成了这样的形式：$$frac{frac{n!}{p ^ {k1}}}{frac{m!}{p ^ {k2}} * frac{(n - m)!}{p ^ {k3}}} * p ^ {k1 - k2 - k3} % p ^ k$$于是现在我们算出(n! \% p ^ k)，然后再乘上(p ^ {k1 -k2 - k3})即可。


2.求(n! \% p ^ k)
这里放出一篇很好的博客（咕咕咕）：【知识总结】扩展卢卡斯定理（exLucas）。
不过这篇博客的代码跑的很慢，原因是在求阶乘的时候没有预处理前缀积，即(f[i] = prod _ {j = 1, (j, p ^ k) = 1} ^ {i} j)，这样递归的每一层都是(O(logn))（快速幂复杂度）的了：$$fac(n) = fac(frac{n}{p}) * f[p ^ k - 1] ^ {frac{n}{p ^ k}} * f[n % p ^ k]$$值得注意的是，这个式子把(n!)中的(p)都提出去了，所以我们求出来的相当于(frac{n}{p ^ {k1}} \% p ^ k)。


到这就没了。
放代码。(luogu P4720 【模板】扩展卢卡斯)

#include<cstdio>
#include<iostream>
#include<cmath>
#include<algorithm>
#include<cstring>
#include<cstdlib>
#include<cctype>
#include<vector>
#include<stack>
#include<queue>
#include<assert.h>
using namespace std;
#define enter puts("") 
#define space putchar(' ')
#define Mem(a, x) memset(a, x, sizeof(a))
#define In inline
typedef long long ll;
typedef double db;
const int INF = 0x3f3f3f3f;
const db eps = 1e-8;
const int maxn = 1e6 + 5;
In ll read()
{
	ll ans = 0;
  	char ch = getchar(), last = ' ';
  	while(!isdigit(ch)) last = ch, ch = getchar();
  	while(isdigit(ch)) ans = (ans << 1) + (ans << 3) + ch - '0', ch = getchar();
  	if(last == '-') ans = -ans;
  	return ans;
}
In void write(ll x)
{
  	if(x < 0) x = -x, putchar('-');
  	if(x >= 10) write(x / 10);
  	putchar(x % 10 + '0');
}
In void MYFILE()
{
#ifndef mrclr
  	freopen(".in", "r", stdin);
  	freopen(".out", "w", stdout);
#endif
}

In ll inc(ll a, ll b, ll mod) {return a + b < mod ? a + b : a + b - mod;}
In ll quickpow(ll a, ll b, ll mod)
{
  	ll ret = 1;
  	for(; b; b >>= 1, a = a * a % mod)
    	if(b & 1) ret = ret * a % mod;
  	return ret;
}
In void exgcd(ll a, ll b, ll& x, ll& y)
{
  	if(!b) x = 1, y = 0;
  	else exgcd(b, a % b, y, x), y -= a / b * x;
}
In ll inv(ll a, ll p)
{
  	ll x, y; exgcd(a, p, x, y);
  	return inc(x % p, p, p);
}

struct ExL
{
  	ll f[maxn];
  	In ll Fac(ll n, ll p, ll pk)
  	{
    	if(!n) return 1;
    	return Fac(n / p, p, pk) * quickpow(f[pk - 1], n / pk, pk) % pk * f[n % pk] % pk;
  	}
  	In ll C(ll n, ll m, ll p, ll pk)
  	{
    	if(n < m) return 0;
    	f[0] = 1;
    	for(int i = 1; i < pk; ++i) f[i] = f[i - 1] * (i % p ? i : 1) % pk;
    	ll f1 = Fac(n, p, pk), f2 = Fac(m, p, pk), f3 = Fac(n - m, p, pk), cnt = 0;
    	//以下是求p ^ (k1 - k2 - k3),其实就是仿照上面求Fac的过程.
    	for(ll i = n; i; i /= p) cnt += i / p;
    	for(ll i = m; i; i /= p) cnt -= i / p;
    	for(ll i = n - m; i; i /= p) cnt -= i / p;
    	return f1 * inv(f2, pk) % pk * inv(f3, pk) % pk * quickpow(p, cnt, pk) % pk;
  	}
  	int cnt = 0;
  	ll a[maxn], c[maxn];
  	In ll CRT()
  	{
    	ll M = 1, ans = 0;
    	for(int i = 1; i <= cnt; ++i) M *= c[i];
    	for(int i = 1; i <= cnt; ++i)
      		ans = inc(ans, a[i] * (M / c[i]) % M * inv(M / c[i], c[i]) % M, M);
    	return ans;
  	}
  	In ll exlucas(ll n, ll m, ll p)
  	{
    	cnt = 0;
    	for(int i = 2; i * i <= p; ++i)
      	{
    		ll tp = 1;
    		while(p % i == 0) p /= i, tp *= i;
    		if(tp > 1) a[++cnt] = C(n, m, i, tp), c[cnt] = tp;
      	}
    	if(p > 1) a[++cnt] = C(n, m, p, p), c[cnt] = p;
    	return CRT();
  }
}E;

int main()
{
  	MYFILE();
  	ll n = read(), m = read(), p = read();
  	write(E.exlucas(n, m, p)), enter;
  	return 0;
}