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.
Programming Praxis 
Link to read does not work, and there is no run link.
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)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).
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!
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 }Here’s a solution in Python.
def prime_factorize(n): assert n >= 2 factors = [] x = 2 while x * x <= n: if n % x == 0: factors.append(x) n //= x continue x += 1 + (x & 1) factors.append(n) 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))