project-euler/110_diophantine_reciprocals...

55 lines
1.2 KiB
Python

#!/usr/bin/env python3
from lib import primegen
from math import log
from operator import mul
from functools import reduce
def calc(pows):
return (reduce(mul, (2*p + 1 for p in pows)) + 1) // 2
def increment(i, pows):
n = 0
if i < len(pows):
n = pows[i]
return [n + 1] * (i+1) + pows[i+1:]
def val(pows):
return sum(p*v for p, v in zip(pows, vals))
def fpwof2(pows, n):
if len(pows) == 1:
return n-1
r = (2*n - 1) // reduce(mul, (2*p + 1 for p in pows[1:]))
r = (r - 1) // 2
while calc([r] + pows[1:]) < n:
r += 1
return r
N = 4 * 10**6
primes = list(primegen(100))
vals = [log(p) / log(2) for p in primes]
i = 0
pows = [0]
best = fpwof2(pows, N)
best_pows = None
while not (val(pows) >= best and sum(pows) == len(pows)):
if not i:
pows[0] = fpwof2(pows, N)
if len(pows) > 1 and pows[0] < pows[1]:
pows[0] = pows[1]
if val(pows) < best:
best = val(pows)
best_pows = pows
i = 1
pows = increment(i, pows)
if val(pows) >= best:
i += 1
else:
i = 0
print(reduce(mul, (n**pw for n, pw in zip(primes, best_pows))))