What is the most performant method to calculate the sum of lcm(i, j), i and j ranging from 1 to k?

THE PROBLEM :

MY SOLUTION :

The solution works, but unfortunately it is really too slow when k is large (k can be a number as big as 69533550916004 for example)

var k = 14;
var res = 0;
for(let i = 1; i <= k; i++){
  for(let j = 1; j <= k; j++){
    res += lcm(i, j);
  }
}
console.log(res)

function lcm(a,b){
  return (a * b) / gcd(a,b);
}

function gcd(a,b){
  if(!b){
    return a; 
  }  
  return gcd(b, a % b); 
}

Thanks

That's a large value of k -- quadratic time won't do. Here's an O(k log log k)-time algorithm, the same as the Sieve of Eratosthenes.

First, some number theory. We know that lcm(i, j) = (i*j) / gcd(i, j). If we were trying to evaluate

 k   k
sum sum (i*j),
i=1 j=1

then we could rewrite using the distributive law:

 k   k            k    2
sum sum (i*j) = (sum i)  = ((k+1)*k / 2)^2.
i=1 j=1          i=1

We can try to salvage the idea by summing over possible greatest common divisors g:

 k   k                k         [k/g]      2
sum sum lcm(i, j) =? sum (1/g)*( sum  (g*i) ),
i=1 j=1              g=1         i=1

but of course this is wrong, since we end up with (e.g.) lcm(2, 2) being represented twice, once with g=1, once with g=2. To straighten it out, we turn to Möbius inversion, which gives us the correct formula

 k   k               k                [k/g]      2
sum sum lcm(i, j) = sum f(g)*(1/g)*( sum  (g*i) ),
i=1 j=1             g=1                i=1

                     k                              2
                  = sum g*f(g)*(([k/g]+1)*[k/g] / 2) ,
                    g=1

where the definition of f(g) is

f(g) =  product  (1-p).
       prime p|g

We can evaluate f in bulk by modifying the Sieve of Eratosthenes. Python 3 below:

def fast_method(k):
    f = [1] * (k + 1)
    for i in range(2, k + 1):
        if f[i] != 1:
            continue
        for j in range(i, k + 1, i):
            f[j] *= 1 - i
    return sum(g * f[g] * ((k // g + 1) * (k // g) // 2) ** 2 for g in range(1, k + 1))


import math


def slow_method(k):
    return sum(math.lcm(i, j) for i in range(1, k + 1) for j in range(k + 1))


def test():
    for k in range(100):
        assert fast_method(k) == slow_method(k), k


test()