Highly Abundant Numbers
December 20, 2016
We studied highly composite numbers in a series of several previous exercises, and had a lot of fun. In today’s exercise we look at a related concept: highly abundant numbers (A002093).
Highly abundant numbers are those numbers n such that sigma(m) < sigma(n) for all m < n, where m and n are positive integers and sigma(n) is the sum of the divisors of n. For instance, 12 is a highly abundant number since the sum of its divisors is 1 + 2 + 3 + 4 + 6 + 12 = 28, and no number less than 12 has a greater sum of divisors.
Your task is to compute the sequence of highly abundant numbers. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.
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