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

Posted by programmingpraxis
Filed in Exercises
9 Comments »

9 Responses to “GCD Sum”

  1. Rutger said
    April 22, 2016 at 1:52 PM

    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
    April 22, 2016 at 3:00 PM

    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
    April 22, 2016 at 6:21 PM

    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[0], len(list(s[1]))) 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
    April 22, 2016 at 7:07 PM

    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

  6. Zack said
    April 22, 2016 at 10:02 PM

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

  7. Daniel said
    April 22, 2016 at 10:29 PM

    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;
}
  8. matthew said
    April 23, 2016 at 7:45 AM

    @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
  9. matthew said
    April 23, 2016 at 8:02 AM

    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.

Leave a comment