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()