扩展卢卡斯定理 笔记

发表于 | 阅读次数：
本文字数： 1.9k 字 | 阅读时长 ≈ 9 分钟

先问一个简单的问题，求：

不保证 $q$ 是质数。$1\le q\le 10^6,1\le m\le n\le 10^{18}$。

（原题：洛谷上的板子）

这很毒瘤了！我们平常熟知的几个方法，预处理阶乘和阶乘逆元，或者是 Lucas 定理，都需要 $q$ 为质数才能这么干。现在 $q$ 连质数都不是，可咋整呢…

正片开始

主题思路

  1. 我们用公式 $C_n^m=\dfrac{n!}{m!(n-m)!}$
  2. 把 $q$ 分解质因数，分解成若干的 $p^k$ 相乘的形式
  3. 现在假设模数为 $p^k$，求 $C_n^m \pmod {p^k}$
  4. 用中国剩余定理合并答案

那么，卢卡斯定理呢？引用某神仙写题解时，说过的一句话：

首先，得确定先会扩展欧几里得算法和扩展中国剩余定理。至于卢卡斯定理，那并不重要。

—— Great_Influence 神仙

（扩展卢卡斯定理和卢卡斯定理并没有任何关系2333）

对每个质因数求答案

我们发现，$n$ 的范围比 $q$ 要大很多。也就是说，分子分母两个阶乘很可能都是 $0$。

那么就出现了一个想法：我们把分子分母中的含 $p$ 的因子都单独拿出来。

设 $n!$ 中最多包含 $a$ 个 $p$ 因子，$m!$ 和 $(n-m)!$ 中分别包含 $b,c$ 个。设 $f(x)$ 表示把 $x!$ 中的 $p$ 因子都去掉后的值，于是答案等于

然后现在我们要求两个东西：

  1. $f(x)$ （定义见上）
  2. $x!$ 中有多少 $p$ 因子，设为 $g(x)$（这个小学就会了吧，$g(x)=g(x/p)+(x/p)$）

那怎么求 $f(x)$ 呢？

第一部分：拿掉所有倍数

先把所有 $p$ 的倍数都拿掉，大概等于：

$[1\times 2\times 3\times …\times (p-1)]\times [(p+1)\times (p+2)\times (p+3)…\times (2p-1)]\times …$

（一直乘到 $x$ 往下最后一个不是 $p$ 的倍数的数）

然后你们就会发现我偷偷的把每 $p-1$ 个数分在了同一个中括号里面。我们称这一个中括号为“一组”。组从 $1$ 开始编号（也就是 $1\times 2\times …\times (p-1)$ 是第 $1$ 组，然后剩下的依次编号）

找一下规律： 第 $x$ 组的积就等于 $[xp-p+1,xp-1]$ 之间的数的积。

然后考虑 $p^k-1$ 这个数，它显然是第 $p^{k-1}$ 组的最后一个。那它下一组，就是以下这些数的乘积：

$p^k+1,p^k+2,… p^k+p-1$。

别忘了我们的 $f(x)$ 模数是 $p^k$ （见上面的式子）。所以，很显然，这些数同余于：

$1,2,…p-1$

这和第一组是一样的了。那就出现了 循环节！

我们设“一节”的积为 $S$，那它就等于第 $1$ 组，第 $2$ 组,…第 $p^{k-1}$ 组中的所有数的乘积。

显然，$x!$ 一共能分出 $\dfrac{x}{p^k}$ 个 $S$。

当然，还有 $x \bmod {p^k}$ 个不完整的“一节”，设为 $R$。暴力计算即可。

这一部分的答案（设为 $A$）就等于 $(S)^{\frac{x}{p^k}} \times R$

那么求 $S$ 和求最后几个不完整的块，就都是 $O(p^k)$ 的了。因为 $q\le 10^6$，所以不会有问题。

第二部分：只考虑 p 的倍数

那这一部分就等于

$p\times 2p \times 3p \times …\times (n/p)p$

然后我们要去掉所有 $p$ 的因子。于是剩下了一个 $(n/p)!$

但是我们发现这个玩意要递归计算…于是这一部分答案为 $f(n/p)$

综上

时间复杂度 $O(p^k\log x)$

中国剩余定理合并答案

什么？你中国剩余定理不会？

额…好的，我会写（gu）博（gu）客（gu）的！

洛谷上的板子代码

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
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
// 以下注释中的 ^ 均表示幂
// （这个代码里面没有异或.......
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cstdarg>
#include <cmath>
using namespace std;
namespace Flandre_Scarlet
{
	#define N 155555
	#define int long long 
	#define F(i,l,r) for(int i=l;i<=r;++i)
	#define D(i,r,l) for(int i=r;i>=l;--i)
	#define Fs(i,l,r,c) for(int i=l;i<=r;c)
	#define Ds(i,r,l,c) for(int i=r;i>=l;c)
	#define MEM(x,a) memset(x,a,sizeof(x))
	#define FK(x) MEM(x,0)
	#define Tra(i,u) for(int i=G.Start(u),v=G.To(i);~i;i=G.Next(i),v=G.To(i))
	#define p_b push_back
	#define sz(a) ((int)a.size())
	#define iter(a,p) (a.begin()+p)
	void R1(int &x)
	{
	    x=0;char c=getchar();int f=1;
	    while(c<'0' or c>'9') f=(c=='-')?-1:1,c=getchar();
	    while(c>='0' and c<='9') x=(x<<1)+(x<<3)+(c^48),c=getchar();
	    x=(f==1)?x:-x;
	}
	void Rd(int cnt,...)
	{
	    va_list args;
	    va_start(args,cnt);
	    F(i,1,cnt) 
	    {
	        int* x=va_arg(args,int*);R1(*x);
	    }
	    va_end(args);
	}

	int n,k,p;
	void Input()
	{
		Rd(3,&n,&k,&p);
	}

	int qpow(int a,int b,int m) // 快速幂，不用细看，求 a^b%m
	{
		int r=1;
		while(b)
		{
			if (b&1) r=r*a%m;
			a=a*a%m,b>>=1;
		}
		return r;
	}
	void exgcd(int a,int b,int &x,int &y) // exgcd，用来求逆元，以及中国剩余定理
	{
		if (b==0) {x=1;y=0;return;}
		int x2,y2;
		exgcd(b,a%b,x2,y2);
		x=y2; y=x2-(a/b)*y2;
	}
	int inv(int a,int p) // 求 a 模 p 的逆元
	{
		int x,y; exgcd(a,p,x,y);
		return (x%p+p)%p;
	}
	int lcm(int a,int b){return a/__gcd(a,b)*b;}
    // 以下函数中的 pk 参数均表示 p^k
    // 分解质因数的时候顺便求的，这样就不用再求一遍快速幂了
    // 常数小~
	int fac(int n,int p,int pk) // 这个就是 f 函数
	{
		if (n==0 or n==1) return 1;

		int r1=1; // r1 就是上面说的 A，表示完整部分
		F(i,1,pk) if (i%p) r1=(r1*i)%pk;
		r1=qpow(r1,n/pk,pk); // 乘一个 n/pk
		int r2=1; // r2 就是上面说的 R，表示剩余部分
        // 写解析的时间和上代码的时间并不同步，所以...名字不同，稍微转化下
		F(i,(n/pk)*pk,n) if (i%p) r2=(r2*(i%pk))%pk;
		return fac(n/p,p,pk)*r1%pk*r2%pk;
	}
	int g(int n,int p) // 求 n! 中有多少 p 因子
	{
		if (n<p) return 0;
		return n/p+g(n/p,p);
	}
	int C(int n,int m,int p,int pk) // 求 C(n,m)%pk
	{
		if (m>n) return 0;
		int f1=fac(n,p,pk);
		int i1=inv(fac(m,p,pk),pk),i2=inv(fac(n-m,p,pk),pk);
		int gg=g(n,p)-g(m,p)-g(n-m,p);
		return f1%pk*i1%pk*i2%pk*qpow(p,gg,pk)%pk;
	}
	int m[N],a[N];
	int cnt=0;
	void solve(int a,int b,int c,int &x,int &y) // 求方程 ax+by=c 的其中一组特殊解，直接修改 x,y 的值并返回
	{
		int g=__gcd(a,b); 
		if (c%g) {x=-1e9; y=-1e9; return;}
		a/=g; b/=g; c/=g;
		exgcd(a,b,x,y);
		x*=c; y*=c;
	}
	int China() 
    // 中国剩余定理（看这伟大的函数名
	{
		int M=m[1],A=a[1];
		F(i,2,cnt)
		{
			int x,y;
			solve(M,m[i],a[i]-A,x,y);
			if (x==-1e9 and y==-1e9) {return -1;}
			x%=m[i];
			A=(A+M*x);
			M=lcm(M,m[i]);
			A=(A%M+M)%M;
		}
		return A;
	}
	int exLucas(int n,int k,int mod) // 正片：求 C(n,k)%mod
	{
		cnt=0;
		for(int i=2;i*i<=mod;++i) if (mod%i==0) // 分解质因数
		{
			int p=i,pk=1;
			while(mod%i==0) pk*=i,mod/=i;
            // 找到一个质因数，就顺便把 p^k 给求了
			++cnt;
			m[cnt]=pk;
			a[cnt]=C(n,k,p,pk)%pk;
		}
		if (mod!=1)
		{
			++cnt;
			m[cnt]=mod;
			a[cnt]=C(n,k,mod,mod)%mod;
		}
		return China();
	}
	void Soviet()
	{
		printf("%lld\n",exLucas(n,k,p));
	}

	#define Flan void
	Flan IsMyWife()
	{
		Input();
		Soviet();
	}
	#undef int //long long
}
int main()
{
	Flandre_Scarlet::IsMyWife();
	getchar();getchar();
	return 0;
}
w