Highly Abundant Numbers

December 20, 2016

We compute the sum of divisors of n by factoring using a 2,3,5-wheel and computing the sum of divisors using the standard formula:

(define (factors n) ; factors of n in increasing order by 2,3,5-wheel
  (define (last-pair xs) (if (null? (cdr xs)) xs (last-pair (cdr xs))))
  (define (cycle . xs) (set-cdr! (last-pair xs) xs) xs)
  (define wheel (cons 1 (cons 2 (cons 2 (cycle 4 2 4 2 4 6 2 6)))))
  (let loop ((n n) (f 2) (wheel wheel) (fs (list)))
    (if (< n (* f f)) (reverse (cons n fs))
      (if (zero? (modulo n f)) (loop (/ n f) f wheel (cons f fs))
        (loop n (+ f (car wheel)) (cdr wheel) fs)))))
(define (sigma n) ; sum of divisors of n
  (if (= n 1) 1
    (let ((fs (factors n)))
      (let loop ((fs (cdr fs)) (f (car fs)) (prod (car fs)) (sum 1))
        (if (null? fs) (* sum (/ (- (* prod f) 1) (- f 1)))
          (if (= (car fs) f) (loop (cdr fs) f (* prod f) sum)
            (loop (cdr fs) (car fs) (car fs)
                  (* sum (/ (- (* prod f) 1) (- f 1))))))))))

You might want to look at https://programmingpraxis.com/2012/02/14/divisors for a different algorithm for computing the sequence of sums of divisors.

Then a highly abundant number is found wherever a new maximum is found in the sequence of sigma. We create a function records that returns the index and value found at each new maximum:

(define (records lt? xs) ; index and value at each new maximum
  (if (null? xs) (error 'records "no data")
    (let loop ((xs (cdr xs)) (k 1) (recs (list (cons 0 (car xs)))))
      (if (null? xs) (reverse recs)
        (if (lt? (cdar recs) (car xs))
            (loop (cdr xs) (+ k 1) (cons (cons k (car xs)) recs))
            (loop (cdr xs) (+ k 1) recs))))))

That makes it easy to compute the sequence of highly abundant numbers:

(define (hans n)
  (map add1 (map car (records < (map sigma (range 1 (+ n 1)))))))
> (hans 2100) ; A002093
(1 2 3 4 6 8 10 12 16 18 20 24 30 36 42 48 60 72 84 90 96 108
120 144 168 180 210 216 240 288 300 336 360 420 480 504 540
600 630 660 720 840 960 1008 1080 1200 1260 1440 1560 1620
1680 1800 1920 1980 2100)

We used range from the Standard Prelude. You can run the program at http://ideone.com/nJWdDu.

Pages: 1 2

3 Responses to “Highly Abundant Numbers”

  1. matthew said

    We can express the divisor sum, sigma(n), as a product of sums of powers of prime divisors, so sigma(60) = sigma(2*2*3*5) = (1+2+2*2)(1+3)(1+5) = 7*4*6 = 168. In python, with a very simple factorizing function:

    def factor(n):
        s = []; p = 2
        while (p*p <= n):
            e = 0
            while (n%p == 0):
                e += 1; n //= p
            if e > 0: s.append((p,e))
            p += 1
        if n > 1: s.append((n,1))
        return s
    
    def sigma(n):
        total = 1
        for (p,e) in factor(n):
            total *= sum([p**i for i in range(0,e+1)])
        return total
    
    max = 0; n = 1;
    while(n < 1000):
        s = sigma(n)
        if s > max:
            print(n,s)
            max = s
        n += 1
    

    However, since we want to find sigma(n) for all n (up to some point) we can use a variation of the Sieve of Eratosthenes to save a lot of factoring. Here, we build up the divisors sums, using the formula above, in array a, using array b to store the sums of powers for each prime p:

    #include <stdio.h>
    #include <stdlib.h>
    
    int main(int argc, char *argv[])
    {
      int N = 1000, max = 0;
      if (argc > 1) N = strtoul(argv[1],NULL,0);
      // a contains products of prime power sums
      // b contains prime power sum for current prime
      int *a = new int[N](), *b = new int[N]();
      // initialize divisor sums to 1
      for (int p = 1; p < N; p++) a[p] = 1;
      for (int p = 2; p < N; p++) {
        if (a[p] == 1) {
          // for each prime < N
          // for each prime power: beware overflow here
          for (long pp = p; pp < N; pp *= p) {
            // add to sum for all multiples
            for (int i = pp; i < N; i += pp) {
              b[i] += pp;
            }
          }
          // for each multiple of p
          for (int i = p; i < N; i += p) {
            a[i] *= 1 + b[i]; // add in the prime power sum
            b[i] = 0;         // reset for future use
          }
        }
        if (a[p] > max) {
          max = a[p];
          printf("%d %d\n", p, a[p]);
        }
      }
    }
    
  2. AN said

    The problem can divide into the following part
    a) find all divisors and it’s sum
    b) keep tracking the current abundant sum until we find the next one

    [/sourcecode lang="css"]
    (define (find-divisor p)
    (define (find-divisor-iter p current divisor-list)
    (cond ((eq? current 0) divisor-list)
    ((eq? (remainder p current) 0) (find-divisor-iter p (- current 1) (cons current divisor-list)))
    (else (find-divisor-iter p (- current 1) divisor-list))))
    (find-divisor-iter p p ‘()))

    (define (print-abundant current current-max)
    (let ((d-list (find-divisor current)))
    (let ((s (apply + d-list)))
    (cond ((> s current-max) (begin (display current) (display “—–“)(display current-max) (newline) (print-abundant (+ 1 current) s)))
    (else (print-abundant (+ 1 current) current-max))))))
    [/sourcecode]

  3. AN said

    Repost due to mis-formated post, why i can’t just edit the old post?

    The problem can divide into the following part
    a) find all divisors and it’s sum
    b) keep tracking the current abundant sum until we find the next one

    (define (find-divisor p)
    (define (find-divisor-iter p current divisor-list)
    (cond ((eq? current 0) divisor-list)
    ((eq? (remainder p current) 0) (find-divisor-iter p (- current 1) (cons current divisor-list)))
    (else (find-divisor-iter p (- current 1) divisor-list))))
    (find-divisor-iter p p ‘()))
    
    (define (print-abundant current current-max)
    (let ((d-list (find-divisor current)))
    (let ((s (apply + d-list)))
    (cond ((> s current-max) (begin (display current) (display “—–“)(display current-max) (newline) (print-abundant (+ 1 current) s)))
    (else (print-abundant (+ 1 current) current-max))))))
    

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 )

Twitter picture

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

Facebook photo

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

Google+ photo

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

Connecting to %s

%d bloggers like this: