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.

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: