from math import gcd,sqrt from functools import reduce # My code actually returns solutions from 2 to max # such that x < n <= y <= z def ErdosStraus(max): def div(n): # efficient divisors finder - taken from stackoverflow - agf, Steinar Lima return sorted(set(x for tup in ([i, n // i] for i in range(2, int(sqrt(n)) + 1) if not n % i) for x in tup)) def extendES(P,M,D): x=M*P y=(M*P+D)//(4*M-1) z=(M*P+D)*M*P//((4*M-1)*D) ES[P]=(x,y,z) def mizony(p): k=((p-1)/4) for a in range(1,p): cc=div((k+a)**2) for c in cc: if (4*k+1)*(4*c+1)%(4*a-1)==0: x = int(k+a) y = int((k+a)*(a+4*k*c+c+k)/(c*(4*a-1))) z = int((4*k+1)*(a+4*k*c+c+k)/(4*a-1)) if sum([i

0: ES[n]=(x,y,z) return ES = {} for p in range(2,max+1): for m in range(4,11): dm=div(m**2) for d in dm: if (p+4*d)%(4*m-1)==0: # general condition if p not in ES.keys(): extendES(p,m,d) k=0 ; P=2 while P<=max: if P not in ES.keys(): extendES(P,m,d) # 1st identity k+=1 P=p+(4*m-1)*k k=0 ; P=2 a=reduce(gcd,(m,d,m**2//d)) while P<=max: if P not in ES.keys(): M=m+m*k//a D=d+d*k//a extendES(P,M,D) # 2nd identity k+=1 P=p+4*(m*p+d)*k//((4*m-1)*a) else: if p not in ES.keys(): if p%3==2: # Using modulo 3 property y=(p+1)//3 z=p*(p+1)//3 ES[p]=(p,y,z) m=2 while m*p<=max: if m*p not in ES: ES[m*p]=(m*p,m*y,m*z) m+=1 else: # Using the multiples d = div(p) if len(d): m = d[0] ES[p] = tuple(m*i for i in ES[p//m]) else: if (p-1)%4==0: # Using mizoni's (n-1)%4 mizony(p) else: basic(p) return ES max=14001 ES=ErdosStraus(max) for n in range(2,max+1): if n not in ES.keys(): print(n,'no solution') else: s=sum([j

Erdős–Straus Challenge - 2