exLucas 算法学习笔记

Posted 2023-05-11 qwq-qaq-tat

一、问题

给定 \\(n,m,p\\)，求 \\(C_n^m\\bmod p\\)，\\(n,m\\le 10^18,p\\le 10^6\\)。
\\(p\\) 现在不一定是质数了，该怎么办？

二、解法

首先，数论题一个常见的做法：如果模数不一定是质数，那就把模数拆成若干个质数的积，然后分别求解，最后用中国剩余定理求解。

1. 拆 \\(p\\)

我们可以把 \\(p\\) 拆成 \\(a_1^b_1\\times a_2^b_2\\times…\\times a_k^b_k\\)，其中所有 \\(a_i\\) 都是质数。
然后对于每一个 \\(a_i^b_i\\)，分别求出 \\(C_n^m \\bmod a_i^b_i\\)，最后再用中国剩余定理合并。

2. 求 \\(C_n^m \\bmod a^b\\)

这等同于求 \\(\\fracn!m!\\times(n-m)!\\bmod a^b\\)。
但是我们不一定能求出 \\(m!\\times (n-m)!\\) 的逆元，因为 \\(m!\\times (n-m)!\\) 不一定与 \\(a\\) 互质。
所以我们要手动让它们互质。
原式可以表示为 \\(\\frac\\fracn!a^x\\fracm!a^y\\times\\frac(n-m)!a^z\\times a^x-y-z\\bmod a^b,x=v_a(n!),y=v_a(m!),z=v_a((n-m)!)\\)。
然后我们可以用 \\(\\operatornameF(n,a)\\) 表示 \\(\\fracn!a^v_a(n!)\\)，用 \\(\\operatornameG(n,a)\\) 表示 \\(v_a(n!)\\)，那么原式就是 \\(\\frac\\operatornameF(n,a)\\operatornameF(m,a)\\times \\operatornameF(n-m,a)\\times a^\\operatornameG(n,a)-\\operatornameG(m,a)-\\operatornameG(n-m,a)\\bmod a^b\\)。

3. 求 \\(\\operatornameG(n,a)\\)

这可以用到 Legendre 公式，也就是 \\(a\\) 是质数，\\(v_a(n!)=\\sum\\limits_i=1\\limits^\\infty\\lfloor\\fracn!a^i\\rfloor\\)。
所以直接按照这个求即可，时间复杂度 \\(O(\\log n)\\)。

4. 求 \\(\\operatornameF(n,a)\\)

由于这个数可能很大，我们求的其实是 \\(\\operatornameF(n,a)\\bmod A\\)，其中 \\(A\\) 就是 \\(a^b\\)。
好像就是 \\(n!\\times \\operatornameinv(a^\\operatornameG(n,a))\\bmod A\\)。但是这样我们要预处理出 \\(1!\\to n!\\)，时间复杂度 \\(O(n)\\)，超时。
考虑用另一种方法计算。我们注意到 \\(n!=1\\times 2\\times …\\times (A-1)\\times A\\times (A+1)\\times …\\times (2A-1)\\times 2A\\times (2A+1)\\times…\\times (\\lfloor\\fracnA\\rfloor\\times A-1)\\times(\\lfloor\\fracnA\\rfloor\\times A)\\times (\\lfloor\\fracnA\\rfloor\\times A+1)\\times…\\times n\\)。
由于 \\(\\operatornameF(n,a)\\) 中 \\(a\\) 的倍数全部被去掉（注意 \\(a\\) 是个质数），所以我们假设 \\(\\operatornameH(s,a)=\\fracsa^v_a(s)\\)。可以发现如果 \\(\\gcd(a,t)=1,t\\mid s\\)，那么 \\(\\operatornameH(s,a)=t\\times \\operatornameH(\\fracst,a)\\)。
所以
\\(\\beginalign\\operatornameF(n,a)&\\equiv\\operatornameH(1\\times 2\\times …\\times (A-1)\\times 1\\times (A+1)\\times …\\times (2A-1)\\times 2\\times (2A+1)\\times…\\times (\\lfloor\\fracnA\\rfloor\\times A-1)\\times\\lfloor\\fracnA\\rfloor\\times (\\lfloor\\fracnA\\rfloor\\times A+1))\\times …\\times n,a)\\pmod A\\\\ &\\equiv\\operatornameH((1\\times 2\\times …\\times (A-1))\\times 1\\times (1\\times 2\\times …\\times (A-1))\\times 2\\times (1\\times…\\times (A-1))\\times\\lfloor\\fracnA\\rfloor\\times (1\\times …\\times (n\\bmod A)),a)\\pmod A\\\\ &\\equiv\\operatornameH(((A-1)!)^\\lfloor\\fracnA\\rfloor\\times (n\\bmod A)!\\times …\\times \\lfloor\\fracnA\\rfloor!,a)\\pmod A\\\\ &\\equiv\\operatornameH(((A-1)!)^\\lfloor\\fracnA\\rfloor\\times (n\\bmod A)!\\times …\\times F(\\lfloor\\fracnA\\rfloor,a),a)\\pmod A\\\\ &\\equiv((A-1)!)^\\lfloor\\fracnA\\rfloor\\times (n\\bmod A)!\\times …\\times F(\\lfloor\\fracnA\\rfloor,a)\\pmod A\\\\ \\endalign\\)
递归调用 \\(O(\\log n)\\) 层，提前预处理出 \\([1\\to (A-1)]!\\bmod A\\) 每一层 \\(O(\\log n)\\)（瓶颈在于快速幂），所以总时间复杂度 \\(O(\\log^2 n)\\)。

5. 总时间复杂度

假设 \\(c\\) 是 \\(p\\) 的不同质因子个数，那么由于 \\(2\\times 3\\times 5\\times 7\\times 11\\times 13\\times 17=1021020>10^6\\)，所以 \\(c\\le 6\\)。

1. 质因数分解，时间复杂度是 \\(O(\\sqrt p)\\)。
2. 对于 \\(i\\in[1,c]\\) 预处理 \\([1\\to (a_i^b_i-1)]!\\bmod a_i^b_i\\)，总时间复杂度是 \\(O(\\sum\\limits_i=1\\limits^ca_i^b_i)\\)。
3. 求 \\(\\operatornameF,\\operatornameG\\)，因为至多调用计算函数 \\(c\\) 次，所以这个部分的上界是 \\(O(c\\log^2 n)\\)。
4. 用中国剩余定理合并，一共合并 \\(c\\) 次，一次时间复杂度 \\(O(\\log p)\\)，总时间复杂度 \\(O(c\\log p)\\)。

所以总时间复杂度是 \\(O(\\sqrt p+c\\log^2 n+c\\log p+\\sum\\limits_i=1\\limits^ca_i^b_i)\\le O(p)\\)，并且常数很小。

三、代码

变量名与上文不完全对应。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
ll t,n,m,p,c=0,s[1000010],a[1010],b[1010];
inline void exgcd(ll &x,ll &y,ll a,ll b)
	if(!b)x=1;y=0;return;
	exgcd(y,x,b,a%b);y-=a/b*x;

inline ll ksm(ll a,ll b,ll p)
	ll r=1;
	while(b)if(b&1)r=r*a%p;a=a*a%p;b>>=1;
	return r;

inline ll inv(ll n,ll p)
	ll x,y;
	exgcd(x,y,n,p);
	return (x%p+p)%p;

inline ll G(ll n,ll p)
	ll r=0;
	while(n)n/=p,r+=n;
	return r;

inline ll F(ll n,ll p,ll pk)
	if(!n)return 1;
	return F(n/p,p,pk)*ksm(s[pk-1],n/pk,pk)%pk*s[n%pk]%pk;

inline ll qwq(ll n,ll m,ll p,ll pk)
	s[0]=1;
	for(ll i=1;i<pk;i++)
		if(i%p)s[i]=s[i-1]*i%pk;
		else s[i]=s[i-1];
	ll nowf=F(n,p,pk)*inv(F(m,p,pk)*F(n-m,p,pk)%pk,pk)%pk;
	ll nowg=ksm(p,G(n,p)-G(m,p)-G(n-m,p),pk);
	return nowf*nowg%pk;

inline ll exlucas(ll x,ll y,ll p)
	ll r=0,P=p;
	for(ll i=2;i*i<=p;i++)
		if(p%i==0)
			a[++c]=1;
			while(p%i==0)p/=i,a[c]*=i;
			b[c]=qwq(x,y,i,a[c]);
		
	if(p>1)a[++c]=p,b[c]=qwq(x,y,p,p);
	for(ll i=1,x,y;i<=c;i++)r=(r+inv(P/a[i],a[i])*P/a[i]%P*b[i]%P)%P;
	return r;

int main()
	cin>>n>>m>>p;
	cout<<exlucas(n,m,p)<<\'\\n\';
	return 0;

Lucas&&Exlucas

Lucas和Exlucas可以求模p意义下大数的组合数。

先考虑p为质数的情况，那么直接上Lucas定理即可。

Lucas 定理基本内容：
[ C_n^m=C_{n mod p}^{m mod p} (mod p) p是质数 ]
对于Lucas的实现直接递归处理即可，注意特判n＜m的情况。

如果p不为质数，那么用Exlucas求解。

实际上，Exlucas和Lucas定理没有关系……

由于p不为质数，所以可以考虑把p质因数分解：
[ p=p_1^{c_1}*p_2^{c_2}...p_t^{c_t} ]

那么可以问题可以转化成求解下列每一个式子的值，然后用CRT合并答案即可：

[ C_n^m mod p_1^{c_1}C_n^m mod p_2^{c_2}...C_n^m mod p_t^{c_t}\]

继续转化，把组合数写成阶乘形式,相当于要求：

[ frac{n!}{m!·(n-m)!} mod p^k ]

发现直接做不好做，逆元问题都不好解决，考虑把左边每一项中的因子p给提出来:

[ frac{frac{n!}{p^{a_1}}}{frac{m!}{p^{a_2}}*frac{(n-m)!}{p^{a3}}}*p^{a^1-a^2-a^3} mod p^k ]

这样逆元就可以直接用Exgcd球了。

那么现在重点解决的问题又变成了：

[ n! mod p^k 其中n还要除去所有的因子p ]

举个例子：n=22，p=3，k=2

把其中所有p（也就是3）的倍数提取出来，得到：
[ 22!=3^7×(1×2×3×4×5×6×7)×(1×2×4×5×7×8)×(10×11×13×14×16×17)×(19×20×22) ]

可以发现：

[ (1×2×4×5×7×8)equiv(10×11×13×14×16×17) (mod 3) ]

所以对于这种一段一段的可以直接处理，最后求一下它（n/p^k）次方即可。

对于3的次方不需要考虑，因为在外面会乘上来。

那么为什么不把3,6彻底分解呢？

这是为了递归的方便，可以发现存在一项7！，也就是（n/p）！，可以递归处理。

所以在一层层递归中，每次每个含有因数p的数提且仅提出一个p，那么可以实现递归，且最后可以把所有p都提出!

然后就可以开心的求出组合数了O(∩_∩)O！