Euclid 与爆搜魔法

· · 算法·理论

View it on my blog :)

设 E(u,v) 为对 u,v（u>v）进行辗转相除求 \gcd 的步数。我们希望研究使得 E(u,v) 较大的 u,v。

众所周知辗转相除在对相邻 Fibonacci 数运行时最劣，因此若设 f(n)=\max_{1\le a< n} E(n,a)，L(N)=\max_{1\le n\le N}f(n)，易知对 L 有 L(N)=\max k\ s.t.\,F_{k+2}\le N。

上面的结论是平凡的。我们尝试考察一点有意思的东西，即使得 E(u,v) 接近最大值的 u,v。具体而言，考察以下几种集合的大小：

\begin{aligned} A_j(N)&=\{n\le N:f(n)=L(N)-j\}\\ B_j(N)&=\{(a,n):1\le a<n\le N,\, E(n,a)=f(n)=L(N)-j\}\\ C_j(N)&=\{(u,v):1\le v<u\le N,\, E(u,v)=L(N)-j\} \end{aligned}

其中 j 较小，设 j\le d（一般当 N 大时 d 不会超过 4）。

连分数与爆搜

首先，辗转相除法的一堆迭代太难考察了，写成规范连分数的形式：u/v=[q_1,\ldots,q_k]，其中 q_k\ge 2；那么有 E(u,v)=k。用连项式描述分子有 u=K(q_1,\ldots,q_k)。记 K_i=K(q_1,\ldots,q_i)，则还有递推 K_i=q_iK_{i-1}+K_{i-2}（K_0=1,K_{-1}=0）。

连分数描述了所有互质对 (u,v)。由于 E(u,v)=E(tu,tv)，可以设 \lambda(n)=\max_{1\le a<n,\, \gcd(a,n)=1}E(n,a)，那么有 f(n)=\max_{u\mid n}\lambda(u)。从互质对出发，计算所有整数对的贡献是容易的；因此以下我们只考虑互质对。

由连分数视角，得到最坏情况是 Fibonacci 是容易的：显然 u 关于任一个 q_i 都是递增的，q_i 全部取最小值得到 u_{\min}=K(1,\ldots,1,2)=F_{k+2}。这也给我们一些启示：固定前一部分系数 q_1,\ldots,q_i（由递推关系我们只需要固定最后两个连项式，A=K_i 和 B=K_{i-1}），那么能够得到的最小分子就是连分数系数全取 1（末项取 2）。设还需要 r=k-i 项，那么归纳递推得到 u_{\min}=K(q_1,\ldots,q_i,1,\ldots,1,2)=F_{r+2}A+F_{r+1}B。

这个东西看起来很像一个剪枝，我们由此出发设计一个爆搜算法：爆搜中维护最后两个连项式，并枚举下一个连分数系数进行转移；如果步数等于 L(N)-j，给相应的集合计算贡献；如果 F_{r+2}A+F_{r+1}B>N，就必定无法产生贡献，直接剪枝。其实这个爆搜约等于反向做辗转相除（区别是辗转相除是连分数系数从后往前递推），只是加上了一个剪枝。下面是一段伪代码（求 |C_j|）：

void dfs(int A, int B, int r) {
    if (r == 0) {
        ans += N / A;
        return;
    }
    for (int q = (r == 1 ? 2 : 1); ; ++q) {
        int C = q * A + B;
        int lower = (r == 1 ? C : F[r + 1] * C + F[r] * A);
        if (lower > N) break;
        dfs(C, A, r - 1);
    }
}

看上去这个爆搜比 O(N^2) 的暴力快不了多少，但是实际上跑得很快，N=10^7,j\le 4 能在几秒钟内出来。以下我们来分析为什么它跑得这么快。

j=0

首先考察 j=0 时的情况。感性理解：这种时候，连分数系数一定不会偏离 Fibonacci 时的太多。因此设 b=(1,\ldots,1,2)。并令偏移量 e_i=q_i-b_i。下面给出一个基于偏移量的估算。

引理（连项式偏移估算）。 令 b=(1,\ldots,1,2)，设 q_i=b_i+e_i，其中 e_i\ge 0，则

K(q_1,\ldots,q_k)\ge F_{k+2}\prod_{i=1}^k(1+e_i/3)

证明。 记 K_{l,r}=K(q_l,\ldots,q_r)，设 A=K_{1,i-1}，B=K_{i+1,k}，A^-=K_{1,i-2}，B^-=K_{i+2,k}。连项式的拼接性质给出：

K(q_1,\ldots,q_k)=Aq_iB+A^-B+AB^-

设 x=A^-/A，y=B^-/B，则由连项式递推有 0\le x,y\le 1。故有 K(q_1,\ldots,q_k)=AB(q_i+x+y)。

现在从 q=b 出发，依次修改每一位 q_i\gets b_i+e_i。每次修改前后比值为

\frac{b_i+e_i+x+y}{b_i+x+y}=1+\frac{e_i}{b_i+x+y}

显然有 b_i+x+y\le 3，因此比值不小于 1+e_i/3。上述比值的下界在每一步都成立，因此得证。

根据以上估计，j=0 时的情况呼之欲出。

引理（j=0 情况）。 f(n)=L(N)（n\le N）仅有两种情况：n 为斐波那契数，或 n=K(q_1,\ldots,q_k) 对应的偏离量 e 中仅有一项为 1，其余全为 0。注意后一种情况仅为必要条件。

证明。 令 L=L(N)，那么 f(n)=L(N) 必须有 n\le N<F_{L+3}，进而有 n/F_{L+2}<F_{L+3}/F_{L+2}\le 5/3。在此基础上，只要规范连分数长度 k=L 即可。

由单调性只需要考虑两种例外：e 中一项为 2，其余为 0；或 e 中两项为 1。由之前的引理得，两者的 n/F_{L+2} 分别至少为 5/3 和 16/9，必然不满足 n/F_{L+2}<5/3。故得证。

当 e_i=1，其余项为 0 时，有 K(q_1,\ldots,q_k)=F_{L+2}+F_iF_{L+2-i}。对该类情况可以直接枚举 i 进行计数，判断分子是否不超过 N 即可，复杂度约为 O(L)=O(\log N)。求 |A_0| 需要把所有得到的 n 排序之后去重，复杂度为 O(L\log L)。

一般情况

下面考虑偏一般的情况。

引理（一般步数下的数量上界）。 设 W_k(N) 为分子不超过 N，且连分数长度为 k 的 u/v 个数。再令 t=\lfloor \log_{4/3}(N/F_{k+2})\rfloor，则有：

W_k(N)\le 2^t\binom{k+t}{k}

证明。 令 R=N/F_{k+2}，由之前的引理，我们需要满足 \prod(1+e_i/3)\le R。

进行一些放缩魔术！令 y_i=\lfloor\log_2(e_i+1)\rfloor，则固定 y_i 时 e_i 共有 2^{y_i} 种取值。放缩一下有：

1+\frac{e_i}3\ge\frac{2^{y_i}+2}3\ge \left(\frac43\right)^{y_i}

那么必定有 (4/3)^{\sum y_i}\le R，即 \sum y_i\le t。枚举 s=\sum y_i\le t，那么固定 s 时的 e 序列数量有上界为

\prod_{i=1}^k2^{y_i}\sum_{\sum y_i=s}1=2^s\binom{s+k-1}{s}

直接求和得到 W_k(N) 的一个上界，再将 2^s 全部放缩成 2^t，得到

W_k(N)\le\sum_{s=0}^t2^s\binom{s+k-1}s\le 2^t\binom{k+t}k

令 L=L(N)。现在对 L-d\le k\le L 的 k 求和，就得到我们希望考察的几个集合大小的一个估计：

\mathcal W(N,d)=\sum_{k=L-d}^{L}W_k(N)

对每个 W_k(N) 的 t 做一个估计：N/F_{k+2}<F_{L+3}/F_{k+2}<2^{L-k+1}\le 2^{d+1}，故 t<(d+1)\log_{4/3}2<3(d+1)。对所有 k 我们统一取 t 的上界 T=3(d+1)，则有

\mathcal W(N,d)\le (d+1)2^T\binom{L+T}{T}\le (d+1)2^T\left(\frac{e(L+T)}{T}\right)^T

取对数有

\log \mathcal W(N,d)\le \log(d+1)+T\log 2+T\log\frac{e(L+T)}T

由 T=O(d)，d\le T-1，右边大概是 O(d\log(1+L/d))，故有

\mathcal W(N,d)=\exp\left(O\left(d\log\left(1+\frac Ld\right)\right)\right)

之前的 dfs 复杂度也就可以得到解释：搜索树叶子节点有 \mathcal W(N,d) 个，搜索深度至多 L，那么复杂度就是 O(L\cdot \mathcal W(N,d)) 的。求 |A_j|，|B_j| 时需要把所有得到的初始对按照分子分类统计，任意初始对一共有 Q\le \lfloor N/F_{L-d+2}\rfloor \mathcal W(N,d)\le 2^d\mathcal W(N,d) 个，故复杂度为 O(L\mathcal W(N,d)+Q\log Q)。

当然我们同样可以估计 \mathcal W(N,d) 的下界。考虑构造一批满足条件的连分数系数列。从长度 k=L-d 的最小序列 (1,\ldots,1,2) 出发，考虑任意选 t\le d 个位置，将其中的 1 改成 2。一个观察是一个 2 对最终分子的贡献不如两个 1（由连项式递推容易验证），因此这样构造出的分子一定不超过 F_{k+t+2}\le F_{L+2}\le N，故是合法的。枚举位置集合得到一个下界：

\mathcal W(N,d)\ge \sum_{t=0}^{\min(d,L-d-1)}\binom{L-d-1}{t}

按 d 分情况缩放。当 d\le L/4 时，只取 t=d 的部分有

\mathcal W(N,d)\ge \binom{L-d-1}{d}\ge\left(\frac L{2d}\right)^d=\exp\left(\Omega\left(d\log\left(1+\frac Ld\right)\right)\right)

当 d>L/4 时，只取 t=\lfloor L/4\rfloor，并取连分数长度 k=L-t。此时有

\mathcal W(N,d)\ge \binom{L-t-1}{t}\ge 2^t=\exp(\Omega(L))

而该情况下又有 d\log(1+L/d)=\Theta(L)。综上可以得到 \mathcal W(N,d) 的比较精确的估计

\mathcal W(N,d)=\exp\left(\Theta\left(d\log\left(1+\frac Ld\right)\right)\right)

稀疏性

以上 \mathcal W(N,d) 估计可以自然得到一个关于辗转相除步数接近上限的密度的推论。

推论（步数近上限密度）。 设 \Delta(n)=L(n)-f(n)。当 N 足够大时，除了一个大小为 o(N) 的集合，对剩下的 n\le N 都有 \Delta(n)=\Omega(\log n)。

证明。 令 L=L(N)。设 H_d(N)=|\{n\le N:\Delta(n)\le d\}|。设 X_k=\min(N,F_{k+d+3}-1)，并考察每个长为 k 的连分数可能做的贡献，有：

\{n\le N:\Delta(n)\le d\}=\bigcup_{k=1}^{L}\bigcup_{\substack{|q|=k\\K(q)\le X_k}}\left\{tK(q):1\le t\le \left\lfloor \frac{X_k}{K(q)}\right\rfloor \right\}

这里的上界估计可以复用之前 W_k(N) 的结果，并集直接用求和估计上界，得到：

H_d(N)\le \sum_{k=1}^L\sum_{\substack{|q|=k\\K(q)\le X_k}}\left\lfloor \frac{X_k}{K(q)}\right\rfloor\le \sum_{k=1}^L\left\lfloor\frac{F_{k+d+3}}{F_{k+2}}\right\rfloor W_k(X_k)\le 2^{d+1}\sum_{k=1}^LW_k(X_k)

接下来的估计过程与 \mathcal W(N,d) 几乎相同。带入之前 W_k(N ) 的估计，并统一取 T=3(d+1)，得到

H_d(N)\le L2^{d+T+1}\binom{L+T}{T}

这就给出

H_d(N)=\exp\left(O\left((d+1)\log\left(1+\frac{L}{d+1}\right)\right)\right)

现在考虑固定 d=\varepsilon L，其中 \varepsilon\to 0，那么得到 \log H_d(N)=O(\varepsilon\log(1+1/\varepsilon)L)。由于 \log N\sim L\log\phi，而 \varepsilon\log(1+1/\varepsilon)\to 0，可以得到 H_d(N)=O(N^{1-\eta})=o(N)。也就是说，当 N 足够大时，除了 H_d(N) 个 n，其余绝大多数 n 都满足 \Delta(n)\ge \varepsilon L=\Omega(\log N)。

以上推论也说明，对随机输入运行欧几里得算法，步数很难达到理论上界。