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

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) 的总和。我怀疑有一些属性可能有用:

如果 a = b，则 R(a,b) = 1

如果 a < b，则 R(a,b) = 1 + R(b,a)

R(a,b) = R(ka,kb) 其中 k 是某个自然数

如果 b <= a，则 R(a,b) = R(a+b,b)

如果 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 的 totient 函数大致是线性的。
我不需要任何特定的代码解决方案，我只需要找出更有效算法的过程。但如果编程语言很重要，我尝试解决这个问题的方法是使用 C++。
旁注:我发现一个 formula已经发现几乎解决了这个问题，但这只是一个近似值。请注意，该公式计算的是平均值而不是总和，因此只需乘以 n^2。如果可以扩展公式以减少错误，它可能会起作用，但据我所知，我不确定这是否可行。

Answer 1

使用 Stern-Brocot，由于对称性，我们可以只查看以 1/3、2/3、3/2 或 3/1 为根的四个子树之一。时间复杂度仍然是 O(n^2) 但显然执行的计算更少。下面的版本使用以 2/3 为根的子树(或者至少这是我仔细考虑的那个:)。另请注意，我们只关心那里的分母，因为分子较低。另请注意，代码也依赖于规则 2 和 3。
C++ 代码(对于 n = 10,000，需要 about a tenth of a second):

#include <iostream>
using namespace std;


long g(int n, int l, int mid, int r, int fromL, int turns){
  long right = 0;
  long left = 0;
        
  if (mid + r <= n)
    right = g(n, mid, mid + r, r, 1, turns + (1^fromL));
          
  if (mid + l <= n)
    left = g(n, l, mid + l, mid, 0, turns + fromL);
    
  // Multiples
  int k = n / mid;

  // This subtree is rooted at 2/3
  return 4 * k * turns + left + right;
}


long f(int n) {  
  // 1/1, 2/2, 3/3 etc.
  long total = n;
      
  // 1/2, 2/4, 3/6 etc.
  if (n > 1)
    total += 3 * (n >> 1);
        
  if (n > 2)
  // Technically 3 turns for 2/3 but
  // we can avoid a subtraction
  // per call by starting with 2. (I
  // guess that means it could be 
  // another subtree, but I haven't
  // thought it through.)
    total += g(n, 2, 3, 1, 1, 2);
    
  return total;
}


int main() {
  cout << f(10000);

  return 0;
}