扩展Lucas

发表于 分类于

求 $\binom{n}{m}\bmod p$ 的值。

$m\le n\le 10^{18},p\le 10^7$

主要是存一个模板。

链接

luogu

引入

Lucas定理：当 $p$ 是质数时，有 $\binom{n}{m}\equiv \binom{\frac{n}{p}}{\frac{m}{p}}\binom{n\bmod p}{m\bmod p}\pmod p$

那么当 $p$ 不是质数时，是否也能计算呢？

算法流程

Part 1

如果要算 $\binom{n}{m}\bmod p$，只需要把 $p$ 质因数分解，$p=\sum_{i=1}^mp_i^{a_i}$。

那么对于任意的 $i$，算出 $\binom{n}{m}\bmod p_i^{a_i}$，然后用CRT合并。

Part 2

现在只需要算组合数对 $p^k$ 取模的结果。

将 $\binom{n}{m}$ 分解成阶乘形式：$\binom{n}{m}=\frac{n!}{m!(n-m)!}$

那么将阶乘搞成 $p^t\times y$ 的形式（其中 $\gcd(y,p)=1$），跟 $p$ 有关的幂用 $n!$ 的减去 $m!$ 和 $(n-m)!$ 的幂就可以算了。

剩下的 $y$ 是有逆元的，可以直接算。

Part 3

现在只需要分解 $n!$ 就好了。

我们将 $n!=1\times 2\times \cdots \times n$ 分成几部分。

第一部分是 $p$ 的倍数，分别是 $1p,2p\cdots \left \lfloor \frac{n}{p} \right \rfloor p$，那么可以表示为 $p^{\left \lfloor \frac{n}{p} \right \rfloor}\times (\left \lfloor \frac{n}{p} \right \rfloor)!$。

第二部分则不是 $p$ 的倍数。这一部分构成了一个模 $p^k$ 的循环节。

由于 $p^k$ 不会超过 $10^7$，因此可以预处理这个循环节。

设 $sum_j\equiv \sum_{i=1,i\bmod p\not = 0}^ji\pmod p^k$

然后就变成 $sum_{p^k}$ 的 $\left \lfloor \frac{n}{p} \right \rfloor$ 次幂乘上 $sum_{n\bmod p^k}$。

时间复杂度

预处理的时间为 $O(p^k)$，算一次组合数的复杂度为 $O(\log_pn\times log_2n)$。

优化

根据不知道是什么定理，$sum_{p^k}$ 要么是 $1$ 要么是 $-1$。

那么算组合数时就省去了快速幂，复杂度变成 $O(\log_2n)$。

程序

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
#include <bits/stdc++.h>
using namespace std;
#define FO(x) freopen(#x".in","r",stdin),freopen(#x".out","w",stdout)
#define fo(i,j,k) for(int i=(j),end_i=(k);i<=end_i;i++)
#define ff(i,j,k) for(int i=(j),end_i=(k);i< end_i;i++)
#define fd(i,j,k) for(int i=(j),end_i=(k);i>=end_i;i--)
#define DEBUG(x) cerr<<#x<<"="<<x<<endl
#define all(x) (x).begin(),(x).end()
#define cle(x) memset(x,0,sizeof(x))
#define lowbit(x) ((x)&-(x))
#define ll long long
#define ull unsigned ll
#define db double
#define lb long db
#define pb push_back
#define mp make_pair
#define fi first
#define se second
inline int read()
{
	int x=0; char ch=getchar(); bool f=0;
	for(;ch<'0'||ch>'9';ch=getchar()) if(ch=='-') f=1;
	for(;ch>='0'&&ch<='9';ch=getchar()) x=(x<<3)+(x<<1)+(ch^48);
	return f?-x:x;
}
#define CASET fo(___,1,read())

inline ll Pow(ll x,ll y,ll mod)
{
	ll ans=1;
	for(;y;y>>=1,x=x*x%mod) if(y&1) ans=ans*x%mod;
	return ans;
}
ll gcd(ll x,ll y){return y?gcd(y,x%y):x;}
ll exgcd(ll a,ll b,ll &x,ll &y)
{
	if(!b)
	{
		x=1; y=0;
		return a;
	}
	ll d=exgcd(b,a%b,y,x);
	y-=a/b*x;
	return d;
}
inline ll inv(ll a,ll b)
{
	ll x,y,d;
	d=exgcd(a,b,x,y);
	if(d!=1) return -1;
	return (x%b+b)%b;
}
inline pair<ll,ll> excrt(int n,ll *c,ll *m)
{
	ll nc=c[1],nm=m[1],d;
	fo(i,2,n)
	{
		d=gcd(m[i],nm);
		if((c[i]-nc)%d!=0) return mp(-1,-1);
		nc=inv(nm/d,m[i]/d)*((c[i]-nc)/d)*nm+nc;
		nm=nm/d*m[i];
		nc=(nc%nm+nm)%nm;
	}
	return mp(nc,nm);
}
const int M=1e6+5;
namespace ExLucas{
	ll sum[M];
	int t,p;
	void init(int _p,int k)
	{
		p=_p; t=1;
		fo(i,1,k) t=t*p;
		sum[0]=1;
		fo(i,1,t) sum[i]=(sum[i-1]*((i%p==0)?1:i))%t;
	}
	int nq,sq;
	ll fac(ll n)
	{
		if(!n) return 1;
		nq+=n/p;
		if(sum[t]==1||(n/t)%2==0) return sum[n%t]*fac(n/p)%t;
		return (t-sum[n%t]%t*fac(n/p)%t)%t;
	}
	inline ll C(ll n,ll m,int p,int k)
	{
		init(p,k); sq=0;
		ll ans,sum;
		nq=0; ans=fac(n); sq+=nq;
		nq=0; sum=fac(m); sq-=nq;
		nq=0; sum=sum*fac(n-m)%t; sq-=nq;
		return ans*inv(sum,t)%t*Pow(p,sq,t)%t;
	}
}
ll n,m,mod,q,c[100],mo[100];
int top;
int main()
{
	cin>>n>>m>>mod; q=mod;
	int t;
	for(int i=2;i*i<=q;i++)
		if(q%i==0)
		{
			for(t=0;q%i==0;q/=i) t++;
			++top;
			c[top]=ExLucas::C(n,m,i,t);
			mo[top]=ExLucas::t;
		}
	if(q!=1)
	{
		++top;
		c[top]=ExLucas::C(n,m,q,1);
		mo[top]=ExLucas::t;
	}
	pair<ll,ll> ans=excrt(top,c,mo);
	printf("%lld",ans.fi);
	return 0;
}