## GCD Sum

### April 22, 2016

`; gcd sum`
```(define (range . args)
(case (length args)
((1) (range 0 (car args) (if (negative? (car args)) -1 1)))
((2) (range (car args) (cadr args) (if (< (car args) (cadr args)) 1 -1)))     ((3) (let ((le? (if (negative? (caddr args)) >= <=)))
(let loop ((x(car args)) (xs '()))
(reverse xs)
(loop (+ x (caddr args)) (cons x xs))))))
(else (error 'range "unrecognized arguments"))))```
```(define (fold-right op base xs)
(if (null? xs)
base
(op (car xs) (fold-right op base (cdr xs)))))```
`(define (sum xs) (fold-right + 0 xs))`
```(define (factors n)
(let loop ((n n) (f 2) (fs (list)))
(cond ((< n (* f f))
(reverse (cons n fs)))
((zero? (modulo n f))
(loop (/ n f) f (cons f fs)))
(else (loop n (+ f 1) fs)))))```
```(define (divisors n)
(let ((divs (list 1)))
(do ((fs (factors n) (cdr fs)))
((null? fs) (sort < divs))
(let ((temp (list)))
(do ((ds divs (cdr ds)))
((null? ds)
(set! divs (append divs temp)))
(let ((d (* (car fs) (car ds))))
(when (not (member d divs))
(set! temp (cons d temp)))))))))```
```(define (sigma x n . args) ; sum of x'th powers of divisors of n
(define (add1 n) (+ n 1))
(define (uniq-c eql? xs)
(if (null? xs) xs
(let loop ((xs (cdr xs)) (prev (car xs)) (k 1) (result '()))
(cond ((null? xs) (reverse (cons (cons prev k) result)))
((eql? (car xs) prev) (loop (cdr xs) prev (+ k 1) result))
(else (loop (cdr xs) (car xs) 1 (cons (cons prev k) result)))))))
(define (prod xs) (apply * xs))
(if (= n 1) 1
(let ((fs (uniq-c = (if (pair? args) (car args) (factors n)))))
(if (zero? x)
(prod (map add1 (map cdr fs)))
(prod (map (lambda (p a)
(/ (- (expt p (* (+ a 1) x)) 1) (- (expt p x) 1)))
(map car fs) (map cdr fs)))))))```
`(define (tau n) (sigma 0 n))`
```(define (totient n) ; count of positive integers less than n coprime to it
(if (= n 1) 1
(let loop ((t n) (p 0) (fs (factors n)))
(if (null? fs) t
(let ((f (car fs)))
(loop (if (= f p) t (* t (/ (- f 1) f))) f (cdr fs)))))))```
```(define (moebius n) ; (-1)^k if n has k factors, or 0 if any factors duplicated
(if (= n 1) 1
(let loop ((m 1) (f 0) (fs (factors n)))
(if (null? fs) m
(if (= f (car fs)) 0
(loop (- m) (car fs) (cdr fs)))))))```
```(define (gcd x y)
(if (zero? y) x
(gcd y (modulo x y))))```
```(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))))```
```(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)))```

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, len(list(s))) 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 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.