## 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 '())) (if (le? (cadr args) x) (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)))

Advertisements

In Python, bruteforce solution:

In Python, td_factors is from Programming Praxis primes essay and divisor_gen generates all divisors.

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.

@Zack: you can just leave out the language specification, which will preserve indentation at least, eg “

” at the end:

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.