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.
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))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.
static long gcdSum(int n) { int[] array = new int[n+1]; // array[0] 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; }@Zack: you can just leave out the language specification, which will preserve indentation at least, eg “
” 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 endGah, 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.