## 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*:

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.

In Python, bruteforce solution:

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

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?):

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

@Zack: https://en.support.wordpress.com/code/posting-source-code/

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

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

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

” at the end:

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.