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.
Programming Praxis 
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 += 1However, 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]); } } }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]
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