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.

Winniewinnifred

Dynamic programming seems to me to be the most promising (and simplest) approach to this problem, though that would still be greater than O(n^2). – Ayotte

I think some people assume everything is a homework problem and downvote for that reason. I've noticed that on the math stackexchange as well. (This question might be better suited for the math stackexchange actually) – Truancy

@Ayotte could you please elaborate a little on your idea for a dynamic program? – Preternatural

@גלעדברקן Basically Just start at (n,n) and do Euclid's algorithm for it, memoizing each intermediate value, then repeat it for all (a={n to 1), b={n to 1}) that does not already have a result. The problem is I cannot think of a way to skip over already completed values without at least checking them, so it's still at least O(n^2) (and probably slightly worse). Doing it bottom-up might be faster, but you'd still need to either check or sum every value in the [n*n] array, so still O(n^2) at best. In it's best form, it'd likely just be the Stern-Brocot tree similar to our answer. – Ayotte

I added C++ code, which seems to take about a tenth of a second for n = 10,000. – Preternatural

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

Regenerate

Can this work (time-wise, practically) for N=10000? If so, could you include it in your snippet? – Killough

Oh, I thought Javascript at least uses 32-bit integer, so we can reach at least 2 billion? Seems just at the right range for the result for 10000. – Killough

Whoops, I tried in Python (after increasing recursion depth), and the result is 782828208, bigger than 32-int. Also it took a few seconds. – Killough

Yea, Javascript is more optimized for recursion (function calls) than Python I think. Also, for n=1000: f(n)=5893024. Oh, you mean you got the same answer in Javascript? That's cool then. (Oh, I misread my number, haha. I thought it's 7 billion. It's only 700 million) – Killough

An iterative version of this might be cool to explore too. – Killough

Answer 2

I think this is a hard problem. We can avoid division and reduce the space usage to linear at least via the Stern--Brocot tree.

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)

Janniejanos

Added an answer that uses symmetry in the tree to at least reduce a little. – Preternatural