【笔记】 exlucas

扩展 lucas

求 (inom{n}{m} mod p)，其中 (1 le m le n le 10^{18},2 le p le 10^6)。

首先我们可以对 (p) 质因数分解 (p=prod p_i^{a_i})，对于每一个质因数我们求出答案，然后用 CRT 合并。

单独考虑每一个模数 (p^k)：

我们要求解 (inom{n}{m} mod p^k)，其中 (1 le m le n le 10^{18},2 le p le 10^6) 且 (p) 为质数。

障碍

找找是什么阻挡了我们求解的脚步？原因是 ((n-m)!m!) 有可能找不到逆元。怎么解决呢？

天真的我曾经以为，当 (n) 足够大的时候，(inom{n}{m} mod p^k=0)，原因是我觉得 (n) 会有很多 (p) 的因子，于是乎 (mod p) 完之后就成了 (0)。

且不论这个想法是否正确，至少它指明了一个方向：求出 (inom{n}{m}) 有多少个 (p) 的因子。

[inom{n}{m} = dfrac{n!}{(n-m)!m!}=dfrac{frac{n!}{p^x}}{frac{(n-m)!}{p^y} imes frac{m!}{p^z}} p^{x-y-z} ]

这样子，前面这部分就不含 (p) 的因子，就可以求出逆元了。

于是我们现在的问题分如下两部分：

求解 (n!) 有多少个 (p) 的因子；
求解 (n!) 去除 (p) 的因子后的值 (mod p^k)。

Task1

不难把 (n!) 拆解成如下部分：

[n!equiv p^{lfloor frac{n}{p} floor}prod_{i=1}^{lfloor frac{n}{p} floor}i imes left(prod_{i=1land i otequiv0pmod p}^{p^k}i ight)^{lfloor frac{n}{p^k} floor} left( prod_{i=p^k+1land i otequiv0 pmod p}^{n}i ight) pmod{p^k} ]

后面两个带括号的部分都不带 (p) 因子，前面 (p^{lfloor frac{n}{p} floor}) 贡献了 (lfloor frac{n}{p} floor)，(prod_{i=1}^{lfloor frac{n}{p} floor}i) 是一个阶乘，可以递归下去。

不难发现可以 (O(log_p n)) 计算 (p) 的指数。

Task2

按照 Task1 的式子计算。

不难发现，预处理一下 (prod_{i=1land i otequiv0pmod p}^{p^k}i) 这个前缀积之后就又可以递归计算了，复杂度 (O(log_p n log n))，多出来一个 (log n) 是因为快速幂，如果需要多次计算组合数的话，可以预处理一下次幂，实现 (O(1)) 快速幂。

计算 (dfrac{frac{n!}{p^x}}{frac{(n-m)!}{p^y} imes frac{m!}{p^z}} p^{x-y-z}) 时要用 exgcd 求逆元，或者算 (varphi) 求快速幂。

最后用 CRT 把各个质数次方的答案合并即可。

Code

#include <bits/stdc++.h>

#define debug(...) fprintf(stderr ,__VA_ARGS__)
#define __FILE(x)
	freopen(#x".in" ,"r" ,stdin);
	freopen(#x".out" ,"w" ,stdout)
#define LL long long

const int MX = 1e6 + 23;

LL read(){
	char k = getchar(); LL x = 0;
	while(k < '0' || k > '9') k = getchar();
	while(k >= '0' && k <= '9') x = x * 10 + k - '0' ,k = getchar();
	return x;
}

LL qpow(LL a ,LL b ,LL p){
	LL ans = 1;
	while(b){if(b & 1) ans = ans * a % p;
		a = a * a % p ,b >>= 1;
	}return ans;
}

void exgcd(LL a ,LL b ,LL &x ,LL &y){
	if(!b) return x = 1 ,y = 0 ,void();
	exgcd(b ,a % b ,y ,x); y -= a / b * x;
}

LL inv(LL x ,LL p){
	LL ret ,tmp;
	exgcd(x ,p ,ret ,tmp);
	return (ret + p) % p;
}

LL rfac[MX];
LL other(LL n ,LL p ,LL pk){
	if(n == 0) return 1LL;
	LL ans = other(n / p ,p ,pk);
	ans = ans * qpow(rfac[pk] ,n / pk ,pk) % pk;
	ans = ans * rfac[n % pk] % pk;
	return ans;
}

LL calcexp(LL n ,LL p){
	if(n < p) return 0LL;
	return n / p + calcexp(n / p ,p);
}

void prework(LL p ,LL pk){
	rfac[0] = 1;
	for(int i = 1 ; i <= pk ; ++i){
		if(i % p == 0) rfac[i] = rfac[i - 1];
		else rfac[i] = rfac[i - 1] * i % pk;
	}
}

LL a[MX] ,xp[MX] ,cnt;
LL b[MX];

LL exlucas(LL n ,LL m ,LL p){
	int f = p;
	for(int i = 2 ; i * i <= f ; ++i){
		if(f % i) continue;
		a[++cnt] = i;
		xp[cnt] = 1;
		while(f % i == 0){
			f /= i;
			xp[cnt] *= i;
		}
	}
	if(f != 1) a[++cnt] = f ,xp[cnt] = f;
	for(int u = 1 ; u <= cnt ; ++u){
		LL P = a[u] ,PK = xp[u];
		prework(P ,PK);
		LL ex = calcexp(n ,P) - calcexp(m ,P) - calcexp(n - m ,P);
	   	LL A = other(n ,P ,PK) * inv(other(m ,P ,PK) ,PK) % PK * inv(other(n - m ,P ,PK) ,PK) % PK;
		b[u] = qpow(P ,ex ,PK) * A % PK;
	}

	LL ans = 0;
	for(int u = 1 ; u <= cnt ; ++u){
		ans = (ans + (p / xp[u]) * inv(p / xp[u] ,xp[u]) % p * b[u]) % p;
	}
	return ans;
}

int main(){
	LL n = read() ,m = read() ,p = read();
	printf("%lld
" ,exlucas(n ,m ,p));
	return 0;
}