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

5 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 https://programmingpraxis.com/2020/01/28/non-abundant-sums/2/, except mine includes 0).

  4. Zack said

    Here is my take on this in Julia: https://pastebin.com/9XL9fWs5

    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 https://programmingpraxis.com/2016/12/20/highly-abundant-numbers/, 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
    }
    

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 )

Google photo

You are commenting using your Google 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 )

Connecting to %s

%d bloggers like this: