数论入门:从 gcd 到 exLucas

· · 算法·理论

Part.0 前言

你是否曾经对数学感到无比恐惧?那太好了,我现在也是。

虽说近几年 NOIP 以下的比赛数学题在大幅减少,毕竟信息学不是数学竞赛,但是掌握一些基础的数学算法对竞赛生涯还是很重要的。

本文介绍信息学竞赛中最为基础的数论算法,包括四个板块:

希望对大家有帮助。

$\mathrm{2026.6.30}$:修改格式,简要讲解时间复杂度。 --- ### 记号约定 通常不区分时间复杂度渐进上界和渐进紧确界。时间复杂度中为显示代码流程有时不省略低阶项。 通常不严格要求各类求和等记号格式。 通常所有记号都是整数,如果没有取整符号,默认取下整。$x\mid y$ 表示 $x$ 整除 $y$。 --- ### 代码 数学题重在思路,代码通常都比较模板,就不展示自己奇丑无比的代码了。仅展示模板代码,尽量保持和正文一样的变量名。 ### 参考文献 * [oi wiki](https://oi.wiki/math/number-theory/basic/) * luogu 题解 # Part.1 欧几里得算法 *本版块从零开始介绍最大公因数,欧几里得算法和拓展欧几里得算法。* --- ### [B4025 最大公约数](https://www.luogu.com.cn/problem/B4025) 如果一组数都能被 $x$ 整除,那么 $x$ 是它们的公约数。最大的 $x$ 称作**最大公约数**,简写作 $\gcd$。 给定 $a,b$,求 $\gcd(a,b)$。 $a,b\leq10^9$。 --- 首先,有欧几里得定理: $$\gcd(a,b)=\gcd(b,a\bmod b)$$ 证明如下: >$a\leq b$ 显然成立。 > >$a>b$ 时,设 $d\mid a,d\mid b$。 > >$$a=\lfloor\frac{a}{b}\rfloor\times b+a\bmod b$$ >$$a\bmod b=a-\lfloor\frac{a}{b}\rfloor\times b$$ >$$\frac{a\bmod b}{d}=\frac{a}{d}-\lfloor\frac{a}{b}\rfloor\times\frac{b}{d}$$ > >右边都是整数,故 $d\mid (a\bmod b)$,即 $a,b$ 的公约数是 $b,a\bmod b$ 的公约数。 > >同理设 $d'\mid b,d'\mid (a\bmod b)$。 >$$\frac{a}{d'}=\frac{a\bmod b}{d'}-\lfloor\frac{a}{b}\rfloor\times\frac{b}{d'}$$ > >故 $d'\mid a$,即 $b,b\bmod a$ 的公约数是 $a,b$ 的公约数。 > >公约数相同,最大公约数也相同,得证。 于是我们有**欧几里得算法**,又称辗转相除法。 递归函数 $\gcd(a,b)$ 流程如下: * 如果 $b$ 非 $0$,则返回 $\gcd(b,a\bmod b)$。 * 如果 $b$ 是 $0$,说明上一步 $b\mid a$,最大公因数自然是上一步 $b$,即这一步的 $a$。 假设 $a,b$ 同阶。在递归过程中,若 $a<b$,则下一层交换;若 $a\geq b$,则 $a\bmod b<\frac{a}{2}$。因此每次递归较大值至少减半,递归深度为 $O(\log a)$,故时间复杂度为 $O(\log a)$。 ```cpp int gcd(int a,int b){ return b ? gcd(b,a % b) : a; } ``` --- ### [P4549 【模板】裴蜀定理](https://www.luogu.com.cn/problem/solution/P4549) 给定长为 $n$ 的序列 $a_i$,求 $x_i$ 使得 $S=\sum\limits_{i=1}^na_ix_i>0$ 且 $S$ 最小。 $n\leq20,\lvert a_i\rvert\leq10^5$。 --- 裴蜀定理,也叫 Bézout 定理,内容为: >对于不全为 $0$ 的 $a,b$,存在 $x,y$ 使得 $ax+by=d$ 的充要条件为 $\gcd(a,b)\mid d$ 。 数学归纳证明: >必要显然。 > >当 $x=1,y=0$ 时显然成立。 > >假设存在 $x',y'$,使得 $x'b+y'(a\bmod b)=d$。 >$$a\bmod b=a-\lfloor\frac{a}{b}\rfloor\times b$$ >$$x'b+y'(a-\lfloor\frac{a}{b}\rfloor\times b)=d$$ >$$ay'+(x'-\lfloor\frac{a}{b}\rfloor)b=d$$ >故得证。 此外多个数也成立: >对于不全为 $0$ 的 $n$ 个数 $a_i$,存在 $x_i$ 使得 $\sum\limits_{i=1}^na_ix_i=d$ 的充要条件为 $\gcd\limits_{i=1}^na_i\mid d$ 。 归纳法可证($\gcd\limits_{i=1}^na_i=\gcd(\gcd\limits_{i=1}^{n-1}a_i,a_n)$)。 对于本题,发现取 $\gcd\limits_{i=1}^n\lvert a_i\rvert$ 是最优的。 --- ### 拓展欧几里得算法 假设有两方程: $$\begin{cases}ax_1+by_1=\gcd(a,b)\\bx_2+(a\bmod b)y_2=\gcd(b,a\bmod b)\end{cases}$$ 根据裴蜀定理可知该方程有解。 根据欧几里得定理知: $$\gcd(a,b)=\gcd(b,a\bmod b)$$ $$ax_1+by_1=bx_2+(a\bmod b)y_2$$ $$a\bmod b=a-\lfloor\frac{a}{b}\rfloor\times b$$ $$ax_1+by_1=bx_2+(a-\lfloor\frac{a}{b}\rfloor\times b)y_2=ay_2+b(x_2-\lfloor\frac{a}{b}\rfloor y_2)$$ $$\begin{cases}x_1=y_2\\y_1=x_2-\lfloor\frac{a}{b}\rfloor y_2\end{cases}$$ 递归求解至 $b=0$ 说明上一步 $b\mid a$ 返回 $x=1,y=0$ 表示 $a\times0+b\times1=\gcd(a,b)=b$。 值得一提的是若 $b\neq 0$,拓展欧几里得给出的解一定 $\lvert x\rvert\leq b,\lvert y\rvert \leq a$。 注意到复杂度不变,为 $O(\log a)$。 ```cpp int exgcd(int a,int b){ int d = a; if(b){ d = exgcd(b,a % b); int t = x; x = y; y = t - a / b * y; } else x = 1,y = 0; return d; } ``` # Part.2 同余方程 *本版块介绍如何利用拓展欧几里得算法处理线性不定方程和同余方程。* --- ### [P1082 [NOIP 2012 提高组] 同余方程](https://www.luogu.com.cn/problem/P1082) 给定 $a,b$,求线性同余方程 $ax\equiv 1\pmod b$ 的最小正整数解。 $a,b\leq2\times10^9$。 --- 转换为二元一次不定方程: $$ax+by=1$$ 而拓展欧几里得算法求的是: $$ax+by=\gcd(a,b)$$ 根据裴蜀定理可知该方程有解当且仅当 $\gcd(a,b)\mid 1$。 可知 $\gcd(a,b)=1$,故有解。除此以为,原方程 $ax\equiv1\pmod b$ 只有恰好有一个解,证明: >假如有两个解 $x_1,x_2$,则: >$$a(x_1-x_2)\equiv0\pmod b$$ >说明 $b\mid a(x_1-x_2)$,而 $a,b$ 互质,故 $b\mid (x_1-x_2)$,则 $x_1\equiv x_2\pmod b$,故只有一个解。 这说明 $i$ 存在乘法逆元 $i^{-1}$ 当且仅当 $i$ 和模数 $b$ 互质,并且只存在 $1$ 个乘法逆元。 --- ### [P5656 【模板】二元一次不定方程 (exgcd)](https://www.luogu.com.cn/problem/P5656) 给定 $a,b,c$,对于二元一次不定方程 $ax+by=c$: - 若该方程无整数解,输出 $-1$。 - 若该方程有整数解,且有正整数解,则输出其**正整数**解的数量、所有**正整数**解中 $x$ 的最小值、所有**正整数**解中 $y$ 的最小值、所有**正整数**解中 $x$ 的最大值、以及所有**正整数**解中 $y$ 的最大值。 - 若方程有整数解,但没有正整数解,你需要输出所有**整数解**中 $x$ 的最小正整数值, $y$ 的最小正整数值。 $a,b,c\leq10^9$。 --- 由裴蜀定理知无解当且仅当 $c$ 不是 $\gcd(a,b)$ 的倍数。 先求: $$ax+by=\gcd(a,b)$$ 拓展欧几里得算法给出的特解记作 $x_0,y_0$。 则原方程特解为: $$\begin{cases}x_1=\frac{cx_0}{\gcd(a,b)}\\y_1=\frac{cy_0}{\gcd(a,b)}\end{cases}$$ 其通解为: $$\begin{cases}x=x_1+t\frac{b}{\gcd(a,b)}\\y=y_1-t\frac{a}{\gcd(a,b)}\end{cases}$$ 其中 $t$ 是整数。 可以将 $x_1$ 对 $\frac{b}{\gcd(a,b)}$ 取模,$y_1$ 对 $\frac{a}{\gcd(a,b)}$ 取模。但如果是 $0$,我们要求正整数解,应当设为模数。 现在要求最小正整数解,以 $x$ 为例: $$x>0$$ $$x_1+t\frac{b}{\gcd(a,b)}>0$$ $$t>-x_1\times\frac{\gcd(a,b)}{b}$$ $$t\geq\lfloor{\frac{-x_1\gcd(a,b)}{b}}\rfloor+1$$ $y$ 同理得: $$t\leq\lceil\frac{y_1\gcd(a,b)}{a}\rceil-1$$ 注意到 $x$ 变大,$y$ 变小。故所求解一定是 $t$ 取边界值时。 # Part.3 同余方程组 *本版块介绍中国剩余定理,即如何处理线性同余方程组。* --- ### [P1495 【模板】中国剩余定理(CRT)/ 曹冲养猪](https://www.luogu.com.cn/problem/P1495) 给定 $n$ 个线性同余方程: $$x\equiv a_i\pmod{p_i}$$ 求该线性同余方程组最小正整数解。 $n\leq10,a_i,p_i\leq10^5,\prod\limits_{i=1}^kp_i\leq10^{18},p_i$ 互质。 --- 先来看中国剩余定理的过程: 计算 $P=\prod\limits_{i=1}^kp_i$。 对于第 $i$ 个方程: * 计算 $m_i=\frac{P}{p_i}$,即 $\prod\limits_{j\neq i}p_j$。 * 计算 $m_i$ 模 $p_i$ 意义下的逆元 $m_i^{-1}$。 * 计算 $c_i=m_im_i^{-1}$,不用取模。 该方程组在模 $P$ 意义下有且只有唯一解 $x=\sum\limits_{i=1}^na_ic_i$。 现在来证明正确性: >因为设计出的 $m_j$ 对于 $j\neq i$,一定是 $p_i$ 倍数,故 $c_j\equiv 0\pmod{p_i}$,于是有: > >$$x\equiv\sum_{i=1}^na_ic_i\equiv a_ic_i\equiv a_i\pmod{p_i}$$ >故得证。 可以知道只有当 $p_i$ 互质时才能使用中国剩余定理,因为否则会不存在逆元。 每个方程都要进行一次拓展欧几里得算法计算逆元,故时间复杂度 $O(n\log P)$。 ```cpp for(int i = 1;i <= n;++i)P *= p[i]; for(int i = 1;i <= n;++i){ m[i] = P / p[i]; exgcd(m[i],p[i]); c[i] = m[i] * (x % p[i] + p[i]) % p[i]; ans = (ans + c[i] % P * a[i] % P) % P; } printf("%lld",ans); ``` 有了中国剩余定理,我们就可以将对合数取模变成对其质因数取模的同余方程组,比如后面要提的拓展 Lucas 定理。 --- ### [P4777 【模板】扩展中国剩余定理(EXCRT)](https://www.luogu.com.cn/problem/P4777) 同上题。 $n\leq10^5,a_i,p_i\leq10^{12},\mathrm{lcm}_{i=1}^np_i\leq10^{18}$。 --- 对于两个方程考虑进行合并: $$\begin{cases}x\equiv a_i\pmod{p_i}\\x\equiv a_j\pmod{p_j}\end{cases}$$ $$x=p_iu+a_i=p_jv+a_j$$ $$p_iu-p_jv=a_j-a_i$$ 根据裴蜀定理知方程有解当且仅当 $$\gcd(p_i,p_j)\mid a_j-a_i$$。 用拓展欧几里得算法得出一组特解 $(u_0,v_0)$。 于是我们将上面两个方程合并成新的方程: $$x\equiv p_iu_0+a_i\pmod{\operatorname{lcm}(p_i,p_j)}$$ 原方程组的解在模 $\operatorname{lcm}_{i=1}^np_i$ 意义下为最后方程的 $a$。 每进行一次拓展欧几里得算法可以合并掉一个方程,时间复杂度 $O(n\log\operatorname{lcm}_{i=1}^np_i)$。 此题建议开`__int128`。 ```cpp #define ii __int128 ii merge(ii ai,ii pi,ii aj,ii pj){ ii d = exgcd(pi,pj); if((aj - ai) % d)return -1; ii u = (x * (aj - ai) / d % (pj / d) + pj / d) % (pj / d); return u * pi + ai; } ii excrt(){ ii A = a[1],P = p[1]; for(int i = 2;i <= n;++i){ A = merge(A,P,a[i],p[i]); if(A == -1)return -1; P = P / gcd(P,p[i]) * p[i]; } return A; } ``` 中国剩余定理其实不常用,我们用的更多的是拓展中国剩余定理,因为它适用范围广,复杂度一致,代码也不是很复杂。 简单总结一下,解决同余方程的主要方法是转化为不定方程,利用拓展欧几里得算法得到特解。解决同余方程组的方法是合并为方程。 --- ### 习题 [P2480 [SDOI2010] 古代猪文](https://www.luogu.com.cn/problem/P2480) [P3868 [TJOI2009] 猜数字](https://www.luogu.com.cn/problem/P3868) [P4774 [NOI2018] 屠龙勇士](https://www.luogu.com.cn/problem/P4774) # Part.4 组合数取模 *前面我们已经知道了逆元、同余方程与中国剩余定理, 现在可以用这些工具来解决竞赛中最常见的计数问题之一:模意义下的组合数。* --- ### 乘法逆元求组合数 $$C_n^m = \frac{n!}{m!(n-m)!}$$ 求 $C_n^m\bmod p$,$q$ 组问询。 $n,m\leq10^7$ 且 $p$ 为质数。 --- 我们可以利用乘法逆元解决,不过要做到线性。 $$\lfloor\frac{p}{i}\rfloor\times i+p\bmod i\equiv0\pmod p$$ 乘以 $i^{-1}$: $$\lfloor\frac{p}{i}\rfloor+p\bmod i\times i^{-1}\equiv0\pmod p$$ 变形得: $$i^{-1}\equiv-\lfloor\frac{p}{i}\times(p\bmod i)^{-1}\rfloor\pmod p$$ 时间复杂度 $O(n+q)$。 ```cpp fac[1] = inv[1] = invf[1] = 1; for(int i = 2;i <= n;++i){ fac[i] = fac[i - 1] * i % p; inv[i] = (-p * i % p * inv[p % i] % p + p) % p; invf[i] = invf[i - 1] * inv[i] % p; } long long C(long long x,long long y){ if(x < y || y < 0)return 0; return fac[x] * invf[y] % p * invf[x - y] % p; } ``` 但是有些毒瘤的题目的数据范围很大。 --- ### [P3807 【模板】卢卡斯定理 / Lucas 定理](https://www.luogu.com.cn/problem/P3807) 同上题。 $n\leq10^{18}$,$p\leq10^7$ 且 $p$ 是质数。 --- 先给出 **Lucas 定理**: $$C_n^m\equiv C_{\lfloor\frac{n}{p}\rfloor}^{\lfloor\frac{m}{p}\rfloor}C_{n\bmod p}^{m\bmod p}\pmod p$$ 证明: >引理:当 $p$ 为质数时,$(x+1)^p\equiv x^p+1\pmod p$。 > >证明: >>根据**二项式定理**展开,系数为 $C_p^k=\frac{p!}{k!(p-k)!}$。当 $k\neq0,p$ 时,$p$ 整除 $C_p^k$,故得证。 > >$$(x+1)^n\equiv((x+1)^p)^{\lfloor\frac{n}{p}\rfloor}(x+1)^{n\bmod p}\equiv(x^p+1)^{\lfloor\frac{n}{p}\rfloor}(x+1)^{n\bmod p}\pmod p$$ > >注意 $x^m$ 的系数,左边为 $C_n^m$,右边两项**只能**取 $(x^p)^{\lfloor\frac{m}{p}\rfloor}$ 和 $x^{m\bmod p}$ 两项,系数为 $C_{\lfloor\frac{n}{p}\rfloor}^{\lfloor\frac{m}{p}\rfloor}\times C_{n\bmod p}^{m\bmod p}$,故得证。 > >请仔细理解上一步的“只能”,~~其实挺好理解的~~。 所以我们可以写一个递归,将 $C_n^m$ 转化为一些规模小于 $p$ 的组合数乘积,这可以用上一题解决。 会有多少个组合数呢? 每递归一层 $n,m$ 会变成 $\lfloor\frac{n}{p}\rfloor,\lfloor\frac{m}{p}\rfloor$,故至多递归 $O(\log_pn)$ 层。时间复杂度 $O(p+q\log_pn)$。 ```cpp long long lucas(long long x,long long y){ if(x < y || y < 0)return 0; return y ? C(x % p,y % p) * lucas(x / p,y / p) % p : 1; } ``` --- ### [P4720 【模板】扩展卢卡斯定理 / exLucas](https://www.luogu.com.cn/problem/P4720) 同上题。 $n\leq10^{18}$,$p\leq10^7$ 且任意。 --- 主要思路:我们可以先分解 $p=\prod\limits_{i=1}^kp_i^{\alpha_i}$,再求出 $V_i\equiv C_n^m\pmod{p_i^{\alpha_i}}$,这是个同余方程组,最后用中国剩余定理合并即可。 怎么求 $V_i$?硬来! $C_n^m=\frac{n!}{m!(n-m)!}$,分母除 $p_i$ 倍数其他都有逆元,故先将 $p_i$ 的倍数分离出来约掉。记 $p_i^{\alpha_i}=P_i$。 定义 $\nu_{p_i}(n)$ 为 $n$ 中质因子 $p_i$ 的个数。对于 $\nu_{p_i}(n!)$,就等于 $\sum\limits_{j=1}^{\lfloor\log_{p_i}n\rfloor}\lfloor\frac{n}{p_i^j}\rfloor$,$O(\log_{p_i}n)$ 可以计算。 定义 $f_{p_i}(n)=\frac{n!}{p_i^{\nu_{p_i}(n!)}}$,即 $n!$ 去除所有质因数 $p_i$。 那么转化为: $$C_n^m\equiv\frac{f_{p_i}(n)}{f_{p_i}(m)f_{p_i}(n-m)}\times p_i^{\nu_{p_i}(n!)-\nu_{p_i}(m!)-\nu_{p_i}((n-m)!)}\pmod{P_i}$$ 现在求 $f_{p_i}(n)$,分离 $p_i$ 倍数(可以去除这个倍数转化为 $f_{p_i}(\lfloor\frac{n}{p_i}\rfloor)$)与非倍数得到: $$f_{p_i}(n)=f_{p_i}(\lfloor\frac{n}{p_i}\rfloor)\prod\limits_{j=1...n,(j,p_i)=1}j$$ 前面那项很好算,后面这项参考之前的方法转化为: $$\prod\limits_{j=1...n,(j,p_i)=1}j\equiv(\prod\limits_{j=1...P_i,(j,p_i)=1}j)^{\lfloor\frac{n}{P_i}\rfloor}(\prod\limits_{j=1,...,n\bmod P_i,(j,p_i)=1}j)\pmod{P_i}$$ 直接暴力算就行,才 $O(P_i)$。 卢卡斯定理部分代码没变,只是多加一个中国剩余定理,总时间复杂度 $O(p+q\sum\limits_i^kP_i\log_{p_i}n)$。 --- ### 习题 [P1641 [SCOI2010] 生成字符串](https://www.luogu.com.cn/problem/P1641) [P7386 「EZEC-6」0-1 Trie](https://www.luogu.com.cn/problem/P7386)