Non-Abundant Sums
January 28, 2020
It’s been a long time since we did an exercise from Project Euler; here is number 23:
A perfect number is a number for which the sum of its proper divisors is exactly equal to the number. For example, the sum of the proper divisors of 28 would be 1 + 2 + 4 + 7 + 14 = 28, which means that 28 is a perfect number.
A number n is called deficient if the sum of its proper divisors is less than n and it is called abundant if this sum exceeds n.
As 12 is the smallest abundant number, 1 + 2 + 3 + 4 + 6 = 16, the smallest number that can be written as the sum of two abundant numbers is 24. By mathematical analysis, it can be shown that all integers greater than 28123 can be written as the sum of two abundant numbers. However, this upper limit cannot be reduced any further by analysis even though it is known that the greatest number that cannot be expressed as the sum of two abundant numbers is less than this limit.
Find the sum of all the positive integers which cannot be written as the sum of two abundant numbers.
Your task is to solve Project Euler 23; in the spirit of Project Euler, please do not display the solution. 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.
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))