【笔记】数论：筛法 2023.8.3
原标题：8.1 数论 + 8.3 类欧
以下记 \(F(n)=\sum_{i=1}^n f(i)\) 为数论函数 \(f\) 的前缀和，或者 \(S(f)\) 也表达同样意思。\(n/m\) 默认表示 \(\left\lfloor\frac{n}{m}\right\rfloor\)。
参考资料
数论柿子题的发现 - 博客 - masterhuang的博客
狄利克雷卷积
积性函数
- 数论函数
f :: Int -> a可以认为是定义域为正整数的函数。 - 积性函数：满足 \(\forall(a,b)=1\Rightarrow f(ab)=f(a)f(b)\) 的数论函数。每个积性函数都应该有 \(f(1)=1\)。
- 完全积性函数：满足 \(\forall a,b,f(ab)=f(a)f(b)\) 的数论函数。
线性筛积性函数
积性函数可以在线性时间内进行线性筛得到每个点的点值，从而计算前缀和。设线性筛积性函数 \(f\)。
- 默认 \(f(p^c)\) 可以关于 \(c\) 线性求。
- 记录 \(low_x\) 表示 \(x\) 的最小质因子 \(p\) 的出现的幂次，\(f(x)=f(x/p^{low_u})f(p^{low_u})\)。
例子：线性筛 \(d,g=d*d\)。
点击查看代码
int pri[1 << 26], cnt, low[1 << 26];
LL D[1 << 26], G[1 << 26];
voId sieve(int N) {
memset(pri, 1, sizeof pri), pri[0] = pri[1] = 0;
D[1] = G[1] = low[1] = 1;
for (int i = 1; i <= N; i++) {
if (pri[i]) {
pri[++cnt] = low[i] = i;
for (LL x = i, j = 1; x <= N; x *= i, j++) {
D[x] = j + 1;
for (int k = 0; k <= j; k++) G[x] += (k + 1) * (j - k + 1);
}
}
for (LL j = 1, t; j <= cnt && (t = i * pri[j]) <= N; j++) {
pri[t] = 0;
if (i % pri[j] == 0) {
low[t] = low[i] * pri[j];
D[t] = D[i / low[i]] * D[low[t]], G[t] = G[i / low[i]] * G[low[t]];
break;
}
low[t] = pri[j];
D[t] = D[i] * D[pri[j]], G[t] = G[i] * G[pri[j]];
}
}
for (int i = 1; i <= N; i++) D[i] += D[i - 1], G[i] += G[i - 1];
}
狄利克雷卷积
对于两个数论函数 \(f,g\)，我们定义它们的狄利克雷卷积 (*) :: (Int -> a) -> (Int -> a) -> (Int -> a) 为 \((f*g)(n)=\sum_{d|n}f(d)g(n/d)\)。（注意这里没有需要积性）
- 交换律 \(f*g=g*f\)。
- 结合律 \(f*(g*h)=(f*g)*h\)。
- 对加法分配律 \((f+g)*h=f*h+g*h\)。
- 对积性函数封闭（积性函数怎么卷都是积性函数）。
- 对点乘分配律 \((f*g)\cdot (h*g)=(f*h)\cdot g\)，注意 \(g\) 必须是完全积性函数。
- 一个由 狄利克雷卷积 构造出的数论函数 \(h=f*g\)，如果 \(f,g\) 可以求得前 \(n\) 项的值，则 \(h\) 的前 \(n\) 项可以在 \(O(n\log n)\) 的时间内计算，这是显然的。
常见积性函数
以下是一些非常重要的完全积性函数：
- 元函数 \(\epsilon(n)=[n=1]\) 是单位元，有 \(f*\epsilon=f\)。
- 恒等函数 \(I(n)=1\)。
- 单位函数 \(Id(n)=n\)。
- 幂函数 \(Id_k(n)=n^k\)。
积性函数：
- 莫比乌斯函数 \(\mu\) 的定义是：有平方因子的数的 \(\mu\) 为 \(0\)，否则为 \((-1)^c\) 其中 \(c\) 是质因子个数。\(\mu*I=\epsilon\)，意思是 \(\mu\) 是 \(I\) 的逆元。
- 欧拉函数 \(\varphi=\sum_{1\leq i\leq n}[(i,n)=1]\)。\(\varphi*I=Id\)，推论是 \(Id*\mu=\varphi\)，\(\varphi*\mu=Id*\mu^2\)（迫真）。
- 因数个数函数 \(d=\sigma_0=I*I=\sum_{d|n}1\)。
- 除数函数 \(\sigma_k=Id_k*I=\sum_{d|n}d^k\)。
补充：常见积性函数关系链
将以上看到的几个函数放在一起得到两条链，链上的一条边如 \(I\to Id\) 表示 \(I*I=Id, Id*\mu=I\)，往右走 \(*I\)，往左走 \(*\mu\)。其实已经有一点接近贝尔级数了。
- \(\mu\to\epsilon\to I\to d\)
- \(\varphi\to Id\)
筛法前置理论
整除分块
- \(n/i\) 取值只有 \(O(\sqrt{n})\) 种，形如：\(1,2,3,\cdots,\sqrt n,n/\sqrt n,n/(\sqrt n-1),\cdots,n\)。
- 封闭性：令 \(S=\{n/x|x\in\mathbb Z\}\)，那么 \(S\) 中的任意一个数 \(x\) 除 \(n\) 的结果 \(n/x\) 仍在 \(S\) 中，且：
- \(x\leq \sqrt n\Rightarrow n/x>\sqrt n\)。
- \(x>\sqrt n\Rightarrow n/x\leq \sqrt n\)。
- \((a/b)/c=a/(bc)\)。
- 满足 \(n/i=n/j\) 的 \(j\) 的最大值为 \(n/(n/i)\)。
于是在计算形如 \(\sum\limits_{i=1}^n f(i)g(n/i)\) 的式子，可以对每一个相同的 \(n/i\) 合并计算，将问题转化为对 \(F,g\) 块筛。
Dirichlet 双曲线法
是整除分块的一种替代品，若 \(f*g=h\)，则 \(H(n)=\sum_{i=1}^nf(i)G(n/i)\)，这是我们的引理，我们可以将引理写成：
尝试以另一种视角计算该式。选定正实数 \(x, y\) 满足 \(xy=n\)。先计算 \(i\leq x\) 的部分，再计算 \(j\leq y\) 的部分，重叠的部分满足 \(i\leq x, j\leq y\implies ij\leq n\)，重叠的都是合法部分，可以直接减去。于是我们有：
复杂度 \(O(x+y)\)，取 \(x=y=\left\lfloor\sqrt n\right\rfloor\) 即可得到最优复杂度。
块筛
对 \(f\) 块筛定义为对固定的 \(n\) 的求出数论函数 \(f\) 的前缀和 \(F\) 的所有 \(n/x(x\in\mathbb Z)\) 处的取值，以下的所有筛法都是在想方设法的在亚线性时间复杂度内解决这个问题，因为线性做法是非常明显的。
- 普通块筛卷积：数论函数 \(f,g\) 可块筛，\(f*g\) 可线性筛，则 \(f*g\) 可在 \(O(\sqrt n)\) 内单点求前缀和，\(O(n^\frac 2 3)\) 块筛。
- 杜教筛：数论函数 \(f,g\) 可块筛，\(f/g\) 可线性筛，则 \(f/g\) 可在 \(O(n^\frac 2 3)\) 内块筛。
- PN 筛：数论函数 \(f,g\)，\(g\) 可块筛，且 \(g(p)=f(p)\)，则 \(f\) 的前缀和可在 \(O(\sqrt n)\) 内单点求值。
普通块筛卷积
若 \(f*g=h\)，则 \(H(n)=\sum_{i=1}^nf(i)G(n/i)\)，即若 \(f,g\) 已完成块筛，则 \(H\) 可以在 \(O(\sqrt n)\) 的时间内完成单点求值。以下证明：
杜教筛
杜教筛
杜教筛可以在低于线性的时间（\(O(n^{\frac 2 3})\)）内块筛数论函数 \(f\)，这需要我们找到另外两个数论函数 \(g,h\) 满足 \(f*g=h\)，并且 \(g,h\) 可以块筛。以下是做法：
时间复杂度
通过一些微积分技巧得到复杂度为（这里 \(T(n/i)\) 被放缩为 \(\sqrt{n/i}\)）。
如果预处理 \(O(n^{2/3})\) 的函数 \(F\) 点值（注意一定要预处理）那么复杂度为 \(O(n^{2/3})\)。
优化：整数除法
可以预处理 \(1\sim \sqrt n\) 每个数的倒数，注意应该写成 (1 + 1e-15) / i，然后再进行除法。经过测试确实快了一点。
代码实现
递归写法是一个容易实现的写法。
点击查看代码
#include <bits/stdc++.h>
using namespace std;
#ifdef LOCAL
#define debug(...) fprintf(stderr, ##__VA_ARGS__)
#else
#define endl "\n"
#define debug(...) voId(0)
#endif
typedef long long LL;
template <int N>
struct siever {
int pri[N + 10], cnt, mu[N + 10], phi[N + 10];
siever() : cnt(0) {
memset(pri, 0x3f, sizeof pri);
pri[0] = pri[1] = false;
mu[1] = phi[1] = 1;
for (int i = 1; i <= N; i++) {
if (pri[i]) {
pri[++cnt] = i;
mu[i] = -1;
phi[i] = i - 1;
}
for (int j = 1; j <= cnt && i * pri[j] <= N; j++) {
pri[i * pri[j]] = false;
if (i % pri[j] == 0) {
mu[i * pri[j]] = 0;
phi[i * pri[j]] = phi[i] * pri[j];
break;
} else {
mu[i * pri[j]] = -mu[i];
phi[i * pri[j]] = phi[i] * phi[pri[j]];
}
}
}
}
};
siever<4000010> sie;
bool vis[1010];
LL Sphi2[1010], Sphi1[4000010];
LL Smu2[1010], Smu1[4000010];
int n;
LL getSphi(unsigned x) {
if (x <= 4e6) return Sphi1[x];
if (vis[n / x]) return Sphi2[n / x];
vis[n / x] = true;
LL &ans = Sphi2[n / x] = 1ll * x * (x + 1) / 2;
for (unsigned l = 2, r; l <= x; l = r + 1) {
r = x / (x / l);
ans -= getSphi(x / l) * (r - l + 1);
}
return ans;
}
LL getSmu(unsigned x) {
if (x <= 4e6) return Smu1[x];
if (vis[n / x]) return Smu2[n / x];
vis[n / x] = true;
LL &ans = Smu2[n / x] = x >= 1;
for (unsigned l = 2, r; l <= x; l = r + 1) {
r = x / (x / l);
ans -= getSmu(x / l) * (r - l + 1);
}
return ans;
}
int main() {
#ifndef LOCAL
cin.tie(nullptr)->sync_with_stdio(false);
#endif
for (int i = 1; i <= 4e6; i++) {
Sphi1[i] = Sphi1[i - 1] + sie.phi[i];
Smu1[i] = Smu1[i - 1] + sie.mu[i];
}
int t;
cin >> t;
while (t--) {
cin >> n;
memset(vis, false, sizeof vis);
cout << getSphi(n) << " ";
memset(vis, false, sizeof vis);
cout << getSmu(n) << endl;
}
return 0;
}
整除分块写法抛弃了递归函数，做到与递归函数相同的时间消耗。
点击查看代码
struct master_dyh {
LL n;
mint* sumg;
vector<mint> sum;
master_dyh() = default;
template <class Func1, class Func2>
master_dyh(LL _n, Func1&& sumf, Func2&& sumh, mint* _sumg)
: n(_n), sumg(_sumg), sum(_n / B + 1) {
for (int i = (int)(n / B); i >= 1; i--) {
LL m = n / i;
sum[i] = sumh(m);
for (LL l = 2, r; l <= m; l = r + 1) {
LL v = m / l;
r = m / v;
sum[i] -= (sumf(r) - sumf(l - 1)) * operator()(v);
}
}
}
mint operator()(LL m) const { return m <= B ? sumg[m] : sum[n / m]; }
};
Dirichlet 双曲线写法抛弃了递归函数，做到 与递归函数相同 不低于递归函数写法（迫真）的时间消耗。
点击查看代码
constexpr int B = 1e7;
double inv[B + 10];
LL dv(LL x, int y) { return x * inv[y]; }
struct master_du {
const LL n;
const mint* sumf;
vector<mint> dp;
master_du() = default;
template <class Func1, class Func2>
master_du(LL _n, mint* _sumf, Func1&& sumg, Func2&& sumh)
: n(_n), sumf(_sumf) {
int upb = dv(n, B);
dp.resize(upb + 1);
for (int i = upb; i >= 1; i--) {
LL m = dv(n, i);
int lim = sqrtl(m);
dp[i] = sumh(m) + sumf[lim] * sumg(lim);
for (int j = 1; j <= lim; j++)
dp[i] -= (sumf[j] - sumf[j - 1]) * sumg(dv(m, j));
for (int j = 2; j <= lim; j++)
dp[i] -= (sumg(j) - sumg(j - 1)) *
(i * j <= upb ? dp[i * j] : sumf[dv(m, j)]);
}
}
mint operator()(LL m) const { return m <= B ? sumf[m] : dp[n / m]; }
};
for (int i = 1; i <= B; i++) inv[i] = (1 + 1e-15) / i;
由于洛谷的【模板】杜教筛太过脏乱，我推荐 P6860 象棋与马 - 洛谷 作为杜教筛模板题。本题答案 \(ans_n=ans_{n/2}+\sum_{i=1}^n\varphi(i), ans_1=0\)，其具体推导过程可观题解。这个形式刚好很符合块筛的胃口，我们只需要块筛 \(\varphi\) 的前缀和即可。
点击查看代码
#include <bits/stdc++.h>
using namespace std;
#ifdef LOCAL
#define debug(...) fprintf(stderr, ##__VA_ARGS__)
#else
#define endl "\n"
#define debug(...) voId(0)
#endif
using LL = long long;
using mint = unsigned long long;
template <int N>
struct siever { /*{{{*/
int pri[N + 10], phi[N + 10], cnt;
siever() : cnt(0) {
memset(pri, true, sizeof pri);
pri[0] = pri[1] = false;
phi[1] = 1;
for (int i = 1; i <= N; i++) {
if (pri[i]) pri[++cnt] = i, phi[i] = i - 1;
for (int j = 1; j <= cnt && i * pri[j] <= N; j++) {
pri[i * pri[j]] = false;
if (i % pri[j] == 0) {
phi[i * pri[j]] = phi[i] * pri[j];
break;
}
phi[i * pri[j]] = phi[i] * (pri[j] - 1);
}
}
}
}; /*}}}*/
constexpr int B = 1e7;
mint sumphi[B + 10];
siever<B> sie;
double inv[B + 10];
LL dv(LL x, int y) { return x * inv[y]; }
struct master_du {
const LL n;
const mint* sumf;
vector<mint> dp;
master_du() = default;
template <class Func1, class Func2>
master_du(LL _n, mint* _sumf, Func1&& sumg, Func2&& sumh)
: n(_n), sumf(_sumf) {
int upb = dv(n, B);
dp.resize(upb + 1);
for (int i = upb; i >= 1; i--) {
LL m = dv(n, i);
int lim = sqrtl(m);
dp[i] = sumh(m) + sumf[lim] * sumg(lim);
for (int j = 1; j <= lim; j++)
dp[i] -= (sumf[j] - sumf[j - 1]) * sumg(dv(m, j));
for (int j = 2; j <= lim; j++)
dp[i] -= (sumg(j) - sumg(j - 1)) *
(i * j <= upb ? dp[i * j] : sumf[dv(m, j)]);
}
}
mint operator()(LL m) const { return m <= B ? sumf[m] : dp[n / m]; }
};
int main() {
#ifndef LOCAL
cin.tie(nullptr)->sync_with_stdio(false);
#endif
for (int i = 1; i <= B; i++) sumphi[i] = sumphi[i - 1] + sie.phi[i];
for (int i = 1; i <= B; i++) inv[i] = (1 + 1e-15) / i;
int t;
cin >> t;
while (t--) {
LL n;
cin >> n;
master_du dyh(
n, sumphi, [&](LL n) { return n; },
[&](mint n) { return n & 1 ? (n + 1) / 2 * n : n / 2 * (n + 1); });
mint ans = 0;
while (n) ans += dyh(n), debug("dyh(%lld) = %llu\n", n, dyh(n)), n >>= 1;
cout << ans - 1 << endl;
}
return 0;
}
Powerful Number 筛
Powerful Number
定义 Powerful Number 为所有满足所有质因子指数都 \(> 1\) 的数，可以证明这样的数一定能被写成 \(a^2b^3\) 的形式，那么这样的数一共不超过 \(O(\int_0^\sqrt{n}\sqrt[3]{n/x^2}\textrm dx)=O(\sqrt{n})\) 个。
想要求出所有 Powerful number，只需要线性筛 \(O(\sqrt n)\) 以内的质数然后 DFS 枚举质数指数。
Powerful Number 筛
Powerful number 筛可以求出积性函数 \(f\) 的前缀和（单点），需要构造积性函数 \(g\) 使得 \(g(p)=f(p)\)（以下 \(p\) 为质数）且 \(g\) 可以块筛。令 \(h*g=f\)，那么有：
- \(h\) 是积性函数。\(h(1)=1\)。
- \(h(p)=0\)，这是因为 \(f(p)=g(1)h(p)+g(p)h(1)=h(p)+g(p)\)，又有 \(g(p)=f(p)\)，所以 \(h(p)=0\)。
- \(h\) 只可能在 Powerful number 处有值。因为如果 \(n\) 不是 Powerful number，那么 \(h(n=pc)=h(p)h(c)=0\)。
从而有：
所以预处理每个 \(h(p^c)\)，爆搜所有 Powerful number 求值即可。
\(h\) 的求法
因为 \(h\) 是积性函数，在 DFS 的过程中我们可以将若干个 \(h(p^c)\) 乘起来就能计算 Powerful number 的点值，所以现在要算 \(h\) 的质数幂点值。两种做法：
- 大力解析 \(h(p^c)\) 使其只和 \(p,c\) 有关。
- \(f=g*h\to f(p^c)=\sum_{k=0}^cg(p^k)h(p^{c-k})\to g(1)h(p^c)=f(p^c)-\sum_{k=1}^cg(p^k)h(p^{c-k})\)。
时间复杂度
对于 \(h\) 第二种求法，复杂度为 \(O(\sqrt{n}\log n)\)。
对于爆搜部分，复杂度为 \(O(\sqrt{n})\)。
例题：SPOJ divcnt3
这题可以使用素数拟合（PN 筛），观察到 \(f(p)=4\)，那么构造 \(g=d*d\)，那么 \(g(p)=4\)。考虑块筛 \(g\)？
\(G(n)=\sum_{i=1}^nd(i)D(n/i)\)，所以要块筛 \(d\)；
\(D(n)=\sum_{i=1}^n d(i)=\sum_{i=1}^n\sum_{d|i}1=\sum_{d=1}^n(n/d)\)。
两个都可以整除分块，都可以线性筛，那么做完。
代码实现：
点击查看代码
#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstring>
#include <vector>
using namespace std;
#ifdef LOCAL
#define debug(...) fprintf(stderr, ##__VA_ARGS__)
#else
#define debug(...) voId(0)
#endif
typedef long long LL;
int pri[1 << 26], cnt, low[1 << 26];
LL D[1 << 26], G[1 << 26];
voId sieve(int N) {
memset(pri, 1, sizeof pri), pri[0] = pri[1] = 0;
D[1] = G[1] = low[1] = 1;
for (int i = 1; i <= N; i++) {
if (pri[i]) {
pri[++cnt] = low[i] = i;
for (LL x = i, j = 1; x <= N; x *= i, j++) {
D[x] = j + 1;
for (int k = 0; k <= j; k++) G[x] += (k + 1) * (j - k + 1);
}
}
for (LL j = 1, t; j <= cnt && (t = i * pri[j]) <= N; j++) {
pri[t] = 0;
if (i % pri[j] == 0) {
low[t] = low[i] * pri[j];
D[t] = D[i / low[i]] * D[low[t]], G[t] = G[i / low[i]] * G[low[t]];
break;
}
low[t] = pri[j];
D[t] = D[i] * D[pri[j]], G[t] = G[i] * G[pri[j]];
}
}
for (int i = 1; i <= N; i++) D[i] += D[i - 1], G[i] += G[i - 1];
}
LL n;
LL rD[1 << 13], rG[1 << 13];
/*
auto remember=[](LL*F,LL*rF){return [&](LL x)->LL&{
debug("call %lld\n",x);
if(x<=6e7) return debug("return F[%lld]=%lld\n",x,F[x]),F[x];
else return debug("return rF[n/%lld]=%lld\n",x,rF[n/x]),rF[n/x];
};};
auto gD=remember(D,rD),gG=remember(G,rG);
*/
LL& gD(LL x) {
if (x <= 6e7)
return D[x];
else
return rD[n / x];
}
LL& gG(LL x) {
if (x <= 6e7)
return G[x];
else
return rG[n / x];
}
LL h[50];
LL dfs(LL x, int now, LL hv) {
LL res = hv * gG(n / x);
for (int i = now; i <= cnt; i++) {
LL p = pri[i];
if (x > n || x * p > n || x * p * p > n) break;
for (LL j = 2, t = x * p * p; t <= n; j++, t *= p)
res += dfs(t, i + 1, hv * h[j]);
}
return res;
}
int main() {
vector<LL> g(60);
for (int i = 0; i <= 50; i++) {
h[i] = 3 * i + 1;
for (int j = 0; j <= i; j++) g[i] += (j + 1) * (i - j + 1);
for (int j = 1; j <= i; j++) h[i] -= g[j] * h[i - j];
}
int T;
scanf("%d", &T);
if (T <= 300)
sieve(6e7);
else
sieve(1e4);
while (T--) {
scanf("%lld", &n);
for (LL i = 1; i * i <= (n); i++) {
LL x = n / i;
if (x <= 6e7) break;
gD(x) = 0;
for (LL l = 1, r; l <= x; l = r + 1)
r = x / (x / l), gD(x) += (r - l + 1) * (x / l);
}
for (LL i = 1; i * i <= (n); i++) {
LL x = n / i;
if (x <= 6e7) break;
gG(x) = 0;
for (LL l = 1, r; l <= x; l = r + 1)
r = x / (x / l), gG(x) += (gD(r) - gD(l - 1)) * gD(x / l);
}
printf("%lld\n", dfs(1, 1, 1));
}
return 0;
}
/*
LL D(LL n) {
LL res = 0;
for (LL l = 1, r; l <= n; l = r + 1) {
r = n / (n / l);
res += (r - l + 1) * (n / l);
}
return res;
}
LL G(LL n) {
LL res = 0;
for (LL l = 1, r; l <= n; l = r + 1) {
r = n / (n / l);
res += (D(r) - D(l - 1)) * D(n / l);
}
return res;
}
*/
后面的内容，有 min_25 筛、类欧几里得算法，都学不动了。标记了。
万能欧几里得
see -> 【模板】万能欧几里得算法 - caijianhong - 博客园
min25 筛
see -> 【模板】min25 筛 - caijianhong - 博客园
类欧几里得算法
以下纯口胡。
直线下整点计数
即计算 \(\sum\limits_{i=0}^n\left\lfloor\frac{ai+b}{c}\right\rfloor\)。
无理数有理逼近
\(\sqrt n\) 逼近为 \(\frac p q\)？
\(a_0=\left\lfloor\sqrt n\right\rfloor\)。
\(a_i=\left\lfloor\frac{1}{\sqrt n-a_{i-1}}\right\rfloor\)。
就是写成了连分数的形式 \(\sqrt n=a_0+\dfrac{1}{a_1+\frac{1}{a_2+\cdots}}\)。
作用是可以计算 \(\left\lfloor i\sqrt r\right\rfloor\) 之类的东西了，有理逼近 \(\sqrt r\) 后用类欧几里得那一套。
模意义有理逼近？
抄了 ix35 课件。
给定素数 \(p\leq 10^{18}\) 和 \(a, b, l, r\)，求 \(x \in [l, r]\) 使得 \((ax + b) \bmod p\) 最小。以下全部在模 \(p\) 意义下讨论。
- 如果 \(a(r − l) < p\)，那么 \((ax+b)\bmod p\) 分为至多两个单调段，暴力即可。
- 否则，一定存在一个 \(x\) 使得 \(ax + b < a\)，于是答案一定不超过 \(a\)，接着让我们关注那些满足 \(ax + b \in[0, a)\)，将它们按照 \(x\) 从小到大排序，则下一项在 \(\pmod a\) 意义下等于上一项减去 \(p \bmod a\)。于是我们把原来的问题 \((a, p)\) 变为了问题 \((−p\bmod a, a)\)，这里 \(b, l, r\) 需要重新计算一下，剩下的无非就是欧几里得过程递归。
- 就是一直跳到 \(al+b\) 前面，然后做分讨？不会。
本文来自博客园，作者：caijianhong，转载请注明原文链接：https://www.cnblogs.com/caijianhong/p/17598863.html
