## 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

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 = +*(N-1)
b = *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 = *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,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; // 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
}
```