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.

Advertisements

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[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

    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[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;
    }
    
  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.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: