Non-Abundant Sums

January 28, 2020

Our solution comes in three parts: first find abundant numbers, second find sums, and third compute the solution. Here we compute the list of abundant numbers less than or equal to 28123, using factors and sigma functions we have written previously, as well as uniq-c from the Standard Prelude:

(define (factors n)
  (let ((wheel (vector 1 2 2 4 2 4 2 4 6 2 6)))
    (let loop ((n n) (f 2) (w 0) (fs (list)))
      (if (< n (* f f)) (reverse (cons n fs))
        (if (zero? (modulo n f))
            (loop (/ n f) f w (cons f fs))
            (loop n (+ f (vector-ref wheel w))
                  (if (= w 10) 3 (+ w 1)) fs))))))
(define (sigma x n) ; sum of x'th powers of divisors of n
  (if (= n 1) 1
    (let ((fs (uniq-c = (factors n))))
      (if (zero? x)
          (apply * (map add1 (map cdr fs)))
          (apply * (map (lambda (p a)
                          (/ (- (expt p (* (+ a 1) x)) 1) (- (expt p x) 1)))
                        (map car fs) (map cdr fs)))))))
(define abuns (filter (lambda (n) (< (+ n n) (sigma 1 n))) (range 1 28124)))

Function 2abun returns two abundant numbers that sum to n, if they exist, or #f if there is no solution:

(define (2abun n)
  (let loop ((abuns abuns))
    (if (null? abuns) #f
      (if (< 28123 (+ (car abuns) (car abuns))) #f
        (if (member (- n (car abuns)) abuns)
            (list (car abuns) (- n (car abuns)))
            (loop (cdr abuns)))))))

We compute the requested sum with the statement (filter (complement 2 abun) (range 1 28124)) in 95 seconds; it’s slow because member repeatedly searches the 6965 elements of the abuns list in linear time. To speed up the program, we convert the abuns list to a hash set:

(define abunset
  (let ((abunset (make-hashtable (lambda (n) n) = 7000)))
    (do ((abuns abuns (cdr abuns))) ((null? abuns) abunset)
      (hashtable-set! abunset (car abuns) (car abuns)))))
(define (2abun n)
  (let loop ((abuns abuns))
    (if (null? abuns) #f
      (if (< 28123 (+ (car abuns) (car abuns))) #f
        (if (hashtable-contains? abunset (- n (car abuns)))
            (list (car abuns) (- n (car abuns)))
            (loop (cdr abuns)))))))

Now we compute the result in about a seventh of a second. There are 1456 positive integers that cannot be written as the sum of two abundant numbers, of which the largest is 20161.

Ideone seems to be having problems; I’ll post the exercise there when it returns.


Pages: 1 2

6 Responses to “Non-Abundant Sums”

  1. Steve said

    Link to read does not work, and there is no run link.

  2. matthew said

    Here’s a Python solution that makes use of the multiplicativity of divisor sum:

    import sympy
    def dsum(N):
        primes = sympy.primerange(2,N)
        a = [0]+[1]*(N-1)
        b = [1]*N
        for p in primes:
            pp = p
            while pp < N:
                for i in range(pp,N,pp):
                    b[i] += pp
                pp *= p
            for i in range(p,N,p):
                a[i] *= b[i]
                b[i] = 1
        return a
    def abundant(N):
        s = dsum(N)
        return [n for n in range(N) if s[n] > 2*n]
    N = 28124
    a = abundant(N)
    b = [1]*N
    for p in a:
        lim = min(p,N-p-1)
        for q in a:
            if q > lim: break
            b[p+q] = 0
    print(sum(b)) # 1457 (counting from zero)
  3. matthew said

    That prints the total number of non-abundant sum numbers (like the Praxis solution at, except mine includes 0).

  4. Zack said

    Here is my take on this in Julia:

    Total run time (including the time-consuming printing of the numbers): 44.5 sec.
    Memory usage: 59.04 MB

    Naturally, this code could be improved, something I may venture another time.

    Cheers for the intriguing problem! Definitely saved me a lot of time going through the PE archive and skimming through all the overly mathematical problems there!

  5. matthew said

    Reusing some code from, here’s a C++ version of that Python solution. Once again, just the count of non-sum numbers is printed:

    #include <stdio.h>
    #include <stdlib.h>
    #include <numeric>
    #include <vector>
    int main(int argc, char *argv[]) {
      int N = (argc == 1) ? 1000 : strtoul(argv[1],NULL,0);
      // a contains products of prime power sums
      // b contains prime power sum for current prime
      // initialize divisor sums to 1
      std::vector<int> a(N,1), b(N,1);
      a[0] = 0; // dsum(0) = 0
      for (int p = 2; p < N; p++) {
        // for each prime < N
        if (a[p] == 1) {
          // 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] *= b[i]; // add in the prime power sum
            b[i] = 1;     // reset for future use
      // Cross off sums in b, which is initially all ones.
      for (int i = 0; i < N; i++) {
        if (a[i] > 2*i) {
          for (int j = 0; j <= i && i+j < N; j++) {
            if (a[j] > 2*j) {
              b[i+j] = 0;
      printf("%d\n", std::accumulate(b.begin(),b.end(),0)); // 1457 for N = 30000
  6. Daniel said

    Here’s a solution in Python.

    def prime_factorize(n):
        assert n >= 2
        factors = []
        x = 2
        while x * x <= n:
            if n % x == 0:
                n //= x
            x += 1 + (x & 1)
        return tuple(factors)
    def calc_divisors(n):
        divisors = {1}
        for factor in prime_factorize(n):
            divisors.update([factor * d for d in divisors])
        return tuple(sorted(divisors))
    abundants = [x for x in range(2, 28124) if sum(calc_divisors(x)[:-1]) > x]
    abundant_sums = {x + y for x in abundants for y in abundants}
    print(sum(x for x in range(1, 28123) if x not in abundant_sums))

Leave a Reply

Fill in your details below or click an icon to log in: Logo

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

Facebook photo

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

Connecting to %s

%d bloggers like this: