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.