欧几里得算法对上限下的数对所采取的步数总和的快速算法

Question

注意：这可能涉及大量的数论，但我在网上找到的公式只是一个近似值，所以我相信一个精确的解决方案需要计算机进行某种迭代计算。

我的目标是找到一种有效的算法（就时间复杂度而言）来解决大 n 值的以下问题：

令 R(a,b) 是欧几里德算法为找到非负整数 a 和 b 的 GCD 所采取的步数。即 R(a,b) = 1 + R(b,a%b)，且 R(a,0) = 0。给定一个自然数 n，求所有 1 的 R(a,b) 之和<= a,b <= n。

例如，如果 n = 2，则解为 R(1,1) + R(1,2) + R(2,1) + R(2,2) = 1 + 2 + 1 + 1 = 5。

由于有 n^2 对对应于要加在一起的数字，无论 R 的效率如何，简单地为每一对计算 R(a,b) 不会比 O(n^2) 好。因此，为了提高为了算法的效率，更快的方法必须以某种方式一次计算多个值的 R(a,b) 的总和。我怀疑有一些属性可能有用：

1. 如果 a = b，则 R(a,b) = 1
2. 如果 a < b，则 R(a,b) = 1 + R(b,a)
3. R(a,b) = R(ka,kb) 其中 k 是某个自然数
4. 如果 b <= a，则 R(a,b) = R(a+b,b)
5. 如果 b <= a < 2b，则 R(a,b) = R(2a-b,a)

由于前两个属性，只需在 a > b 的对上找到 R(a,b) 的总和。除了 a 大于 b 之外，我尝试在计算 R(a,b) 的方法中除了第三个属性外使用它，仅针对其中 a 和 b 也是互质的对。总和是 n 加上 (n / a) * ((2 * R(a,b)) + 1) 在所有这些对上的总和（对 n / a 使用整数除法）。我发现这个算法的时间复杂度仍然是 O(n^2)，因为 Euler 的整体函数大致是线性的。

我不需要任何特定的代码解决方案，我只需要找出更有效算法的过程。但如果编程语言很重要，我尝试解决这个问题的方法是使用 C++。

旁注：我发现发现了一个几乎可以解决这个问题的公式，但这只是一个近似值。请注意，该公式计算的是平均值而不是总和，因此只需乘以 n^2。如果可以扩展公式以减少错误，它可能会起作用，但据我所知，我不确定这是否可行。

Answer 1

我认为这是一个难题。我们至少可以通过 Stern--Brocot 树来避免划分并将空间使用减少到线性。

def f(n, a, b, r):
    return r if a + b > n else r + f(n, a + b, b, r) + f(n, a + b, a, r + 1)


def R_sum(n):
    return sum(f(n, d, d, 1) for d in range(1, n + 1))


def R(a, b):
    return 1 + R(b, a % b) if b else 0


def test(n):
    print(R_sum(n))
    print(sum(R(a, b) for a in range(1, n + 1) for b in range(1, n + 1)))


test(100)