Fast algorithm for sum of steps taken by the Euclidean algorithm over pairs of numbers under an upper bound

Question

Note: This may involve a good deal of number theory, but the formula I found online is only an approximation, so I believe an exact solution requires some sort of iterative calculation by a computer.

My goal is to find an efficient algorithm (in terms of time complexity) to solve the following problem for large values of n:

Let R(a,b) be the amount of steps that the Euclidean algorithm takes to find the GCD of nonnegative integers a and b. That is, R(a,b) = 1 + R(b,a%b), and R(a,0) = 0. Given a natural number n, find the sum of R(a,b) for all 1 <= a,b <= n.

For example, if n = 2, then the solution is R(1,1) + R(1,2) + R(2,1) + R(2,2) = 1 + 2 + 1 + 1 = 5.

Since there are n^2 pairs corresponding to the numbers to be added together, simply computing R(a,b) for every pair can do no better than O(n^2), regardless of the efficiency of R. Thus, to improve the efficiency of the algorithm, a faster method must somehow calculate the sum of R(a,b) over many values at once. There are a few properties that I suspect might be useful:

1. If a = b, then R(a,b) = 1
2. If a < b, then R(a,b) = 1 + R(b,a)
3. R(a,b) = R(ka,kb) where k is some natural number
4. If b <= a, then R(a,b) = R(a+b,b)
5. If b <= a < 2b, then R(a,b) = R(2a-b,a)

Because of the first two properties, it is only necessary to find the sum of R(a,b) over pairs where a > b. I tried using this in addition to the third property in a method that computes R(a,b) only for pairs where a and b are also coprime in addition to a being greater than b. The total sum is then n plus the sum of (n / a) * ((2 * R(a,b)) + 1) over all such pairs (using integer division for n / a). This algorithm still had time complexity O(n^2), I discovered, due to Euler's totient function being roughly linear.

I don't need any specific code solution, I just need to figure out the procedure for a more efficient algorithm. But if the programming language matters, my attempts to solve this problem have used C++.

Side note: I have found that a formula has been discovered that nearly solves this problem, but it is only an approximation. Note that the formula calculates the average rather than the sum, so it would just need to be multiplied by n^2. If the formula could be expanded to reduce the error, it might work, but from what I can tell, I'm not sure if this is possible.

Answer 1

Using Stern-Brocot, due to symmetry, we can look at just one of the four subtrees rooted at 1/3, 2/3, 3/2 or 3/1. The time complexity is still O(n^2) but obviously performs less calculations. The version below uses the subtree rooted at 2/3 (or at least that's the one I looked at to think through :). Also note, we only care about the denominators there since the numerators are lower. Also note the code relies on rules 2 and 3 as well.

C++ code (takes about a tenth of a second for n = 10,000):

#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;
}