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

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: