GCD Sum

April 22, 2016

Today’s exercise is inspired by A018804: Find the sum of the greatest common divisors gcd(k, n) for 1 ≤ k ≤ n.

Your task is to write a function that finds the sum of the greatest common divisors. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

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