A natural number, N, that can be written as the sum and product of a given set of at least two natural numbers, {a1, a2, ... , ak} is called a product-sum number: N = a1 + a2 + ... + ak = a1 × a2 × ... × ak.
For example, 6 = 1 + 2 + 3 = 1 × 2 × 3.
For a given set of size, k, we shall call the smallest N with this property a minimal product-sum number. The minimal product-sum numbers for sets of size, k = 2, 3, 4, 5, and 6 are as follows.
k=2: 4 = 2 × 2 = 2 + 2
k=3: 6 = 1 × 2 × 3 = 1 + 2 + 3
k=4: 8 = 1 × 1 × 2 × 4 = 1 + 1 + 2 + 4
k=5: 8 = 1 × 1 × 2 × 2 × 2 = 1 + 1 + 2 + 2 + 2
k=6: 12 = 1 × 1 × 1 × 1 × 2 × 6 = 1 + 1 + 1 + 1 + 2 + 6
Hence for 2≤k≤6, the sum of all the minimal product-sum numbers is 4+6+8+12 = 30; note that 8 is only counted once in the sum.
In fact, as the complete set of minimal product-sum numbers for 2≤k≤12 is {4, 6, 8, 12, 15, 16}, the sum is 61.
What is the sum of all the minimal product-sum numbers for 2≤k≤12000?
In [1]:
import itertools
import numpy as np
def prime_factors(n): # from eratosthenes
factors = []
d = 2
while n > 1:
while n % d == 0:
factors.append(d)
n /= d
d = d + 1
if d * d > n:
if n > 1:
factors.append(n)
break
return factors
def powerset(iterable): # from itertools
"powerset([1,2,3]) --> () (1,) (2,) (3,) (1,2) (1,3) (2,3) (1,2,3)"
s = list(iterable)
return itertools.chain.from_iterable(itertools.combinations(s, r) for r in range(len(s)+1))
def proper_divisors(n):
pf = prime_factors(n)
factors = []
for p in set(powerset(pf)):
factors.append(int(np.prod(np.array(p))))
return sorted(factors)[:-1]
print proper_divisors(n=220)
In [141]:
import math
import time
import calendar
def factorizations(n):
divisors = proper_divisors(n)
facts = set()
for i,d in enumerate(divisors):
if d == 1:
continue
b = n/d
facts.add(tuple(sorted([b,d])))
for f in factorizations(b):
f = list(f)
f.append(d)
facts.add(tuple(sorted(f)))
return facts
def is_prod_sum(n, k):
for f in factorizations(n):
if len(f) == k and sum(f) == n:
#print f
return True
elif len(f) < k and sum(f) == n - (k - len(f)):
#print f, '+ 1 *', k - len(f)
return True
return False
time_start = calendar.timegm(time.gmtime())
min_set = set()
max_k = 12000
for k in range(2,max_k):
for n in range(k,int(max_k*1.5)):
if is_prod_sum(n, k):
#print "k =",k,"n =",n
min_set.add(n)
break
print sum(min_set), len(min_set)
time_stop = calendar.timegm(time.gmtime())
print time_stop - time_start, 'seconds'
In [119]:
max(min_set) # 12200. 7587457. 1409
Out[119]:
In [ ]: