## GCD Sum

### April 22, 2016

We begin with a simple recursive function to find the greatest common divisor of two positive integers; this algorithm was discovered by Euclid about 300BC:

(define (gcd x y)
(if (zero? y) x
(gcd y (modulo x y))))

Then the brute force solution simply sums the gcd’s; here are four ways to write it:

(define (gcd-sum n)
(let loop ((k 1) (s 0))
(if (< n k) s
(loop (+ k 1) (+ s (gcd k n))))))
(define (gcd-sum n)
(let ((sum 0))
(do ((k 1 (+ k 1))) ((< n k) sum)
(set! sum (+ sum (gcd k n))))))
(define (gcd-sum n)
(sum (map (lambda (k) (gcd k n)) (range 1 (+ n 1)))))
(define (gcd-sum n)
(fold-right (lambda (k s) (+ s (gcd k n))) 0 (range 1 (+ n 1))))

Any of these methods are fine if n is small, but quickly grow slower as n grows larger:

> (time (gcd-sum (* 17 97)))
(time (gcd-sum (* 17 ...)))
no collections
2 ms elapsed cpu time
1 ms elapsed real time
144 bytes allocated
6369
> (time (gcd-sum (* 97 997)))
(time (gcd-sum (* 97 ...)))
no collections
118 ms elapsed cpu time
119 ms elapsed real time
144 bytes allocated
384649
> (time (gcd-sum (* 991 997)))
(time (gcd-sum (* 991 ...)))
no collections
1394 ms elapsed cpu time
1396 ms elapsed real time
144 bytes allocated
3948133

If you can factor n, looking at A018804 will show you two formulas that greatly speed the computation for large n: $\sum_{k=1}^n\gcd(k,n)=\sum_{d|n}d\,\phi(n/d)=\sum_{d|n}d\,\tau(d)\,\mu(n/d)$

where φ(n) is Euler’s totient function, τ(n) is the number of divisors and μ(n) is the möbius function.

(define (gcd-sum n)
(fold-right
(lambda (d s) (+ s (* d (totient (/ n d)))))
0
(divisors n)))
(define (gcd-sum n)
(let ((ds (divisors n)))
(fold-right
(lambda (d s) (+ s (* d (tau d) (moebius (/ n d)))))
0
ds)))

Either of these is very much faster than the brute-force method:

> (time (gcd-sum (* 991 997)))
(time (gcd-sum (* 991 ...)))
no collections
0 ms elapsed cpu time
1 ms elapsed real time
2512 bytes allocated
3948133

You can run the program at http://ideone.com/UgFgSy; the entire program, including the number-theoretic functions, is also shown on the next page.

Pages: 1 2 3

### 9 Responses to “GCD Sum”

1. Rutger said

In Python, bruteforce solution:

from fractions import gcd

for n in range(20):
print n, ":", sum([gcd(n,k) for k in range(1,n+1)])

# Result
# 0 : 0
# 1 : 1
# 2 : 3
# 3 : 5
# 4 : 8
# 5 : 9
# 6 : 15
# 7 : 13
# 8 : 20
# 9 : 21
# 10 : 27
# 11 : 21
# 12 : 40
# 13 : 25
# 14 : 39
# 15 : 45
# 16 : 48
# 17 : 33
# 18 : 63
# 19 : 37
# [Finished in 0.2s]

2. Paul said

In Python, td_factors is from Programming Praxis primes essay and divisor_gen generates all divisors.

def totient(n):
product = n
for prime in set(td_factors(n)):
product -= product // prime
return product

def sumgcd(n):
return sum(d * totient(n // d) for d in divisor_gen(n))

3. matthew said

http://oeis.org/A018804 also tells us that sumgcd (or Pillai’s arithmetic function) is multiplicative (f is multiplicative if f(nm) = f(n)f(m) for coprime n and m) and that its value for p^e is p^(e-1)*((p-1)e+p), for prime p – this suffices to calculate the function for any n, given a prime factorization. Here we use the Python primefac library, converting the output to prime-exponent pairs, and use that to define a generic multiplicative function. Fairly snappy, even for “colossally abundant” numbers with lots of divisors (has anyone seen the Ramanujan film yet?):

from primefac import primefac
from itertools import groupby
from operator import mul

def factors(n): return ((s, len(list(s))) for s in groupby(primefac(n)))
def multiplicative(f): return lambda n: reduce(mul, (f(p,e) for (p,e) in factors(n)))
gcdsum = multiplicative(lambda p,e: p**(e-1)*((p-1)*e+p))

print gcdsum(991*997)
print gcdsum(224403121196654400) # 22nd colossally abundant number

\$ ./gcdsum.py
3948133
4812563181512005920000

4. Zack said

Here is an alternative in Julia, without invoking a bunch of packages to do all the hard work:

function gcd_(x::Int64, y::Int64)
m = min(x,y)

for d in m:-1:1
if (x%d == 0) && (y%d == 0)
return d
end
end
end

function gcd_sum(n::Int64)
s = n + 1

for k = 2:(n-1)
s += gcd(k,n)
end

return s
end

5. Zack said

@matthew: This CSS add-on doesn’t support formatting for Julia…

6. Daniel said

Here’s a Java solution that is similar in spirit to Sieve of Eratosthenes.

static long gcdSum(int n) {
int[] array = new int[n+1]; // array ignored

for (int i = 1; i <= n/2; i++) {
if (n % i == 0) {
for (int j = i; j <= n; j += i) {
array[j] = i;
}
}
}

long sum = n;
for (int i = 1; i < n; i++) {
sum += array[i];
}

return sum;
}

7. matthew said

@Zack: you can just leave out the language specification, which will preserve indentation at least, eg “

" at the start, "

” at the end:

function gcd_(x::Int64, y::Int64)
m = min(x,y)

for d in m:-1:1
if (x%d == 0) && (y%d == 0)
return d
end
end
end

function gcd_sum(n::Int64)
s = n + 1

for k = 2:(n-1)
s += gcd(k,n)
end

return s
end

8. matthew said

Gah, I was attempting to show the formatting commands [code] and [/code] – I had a notion they only applied when on a new line, but evidently not. Using HTML entities for the square brackets here.