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

问题描述

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

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

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

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

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

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

由于前两个性质，只需要在 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。如果可以扩展公式以减少错误，它可能会起作用，但据我所知，我不确定这是否可行。

解决方法

使用 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,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,2);
    
  return total;
}


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

  return 0;
}

,

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

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


def R_sum(n):
    return sum(f(n,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)