from math import sqrt from lib_misc import gcd def proper_divisors(n): """ Returns the list of divisors for n excluding n. """ if n < 2: return [] divisors = [1, ] d = 2 while d * d <= n: if n % d == 0: divisors.append(d) d += 1 return divisors def triples_brute_force(b_max): squares = {n * n: n for n in range(1, b_max)} triples = [] for b in range(3, b_max): for a in range(2, b): c_squared = a * a + b * b if c_squared in squares: triples.append((a, b, squares[c_squared])) return triples def dicksons_method(r): triples = [] assert(r % 2 == 0) st = r * r // 2 for d in proper_divisors(st): s = d t = st // s a = r + s b = r + t c = r + s + t L = a + b + c triples.append(L) return sorted(triples) def dickson_method_min_max_perimeter(r): """ I used this method to check that the min and max perimeter for increasing r are also increasing strictly monotonic. """ assert(r % 2 == 0) st = r * r // 2 L_min = 10**8 L_max = 0 for d in proper_divisors(st): s = d t = st // s L = r + s + r + t + r + s + t if L > L_max: L_max = L if L < L_min: L_min = L return (L_min, L_max) def dicksons_method_efficient(r): L_max = 1500000 perimeters = [] assert(r % 2 == 0) st = r * r // 2 s_max = int(sqrt(st)) for s in range(s_max, 0, -1): if st % s == 0: t = st // s L = r + s + r + t + r + s + t if L > L_max: break perimeters.append(L) else: pass return perimeters def euler_075_dicksons(): """ We are using Dickson's method to get all Pythagorean triples. The challenge is that this method requires finding the divisors for a given r. To make it more efficient we stop iterating after the perimeter exceeds the threshold L_max for any given r. In addition we found a r_max of 258000 empirically. This runs for 15 minutes with cpython and in under 5 minutes in pypy. It would have been better to use a generator function that finds native solutions and then extrapolate with factors k in range(2, 1500000) while L < L_max. https://en.wikipedia.org/wiki/Formulas_for_generating_Pythagorean_triples#Dickson's_method XXX: Can be strongly optimized. """ r_max = 258000 perimeter_counts = {} for r in range(2, r_max, 2): if r % 10000 == 0: print(r) for perimeter in dicksons_method_efficient(r): try: perimeter_counts[perimeter] += 1 except KeyError: perimeter_counts[perimeter] = 1 one_integer_sided_right_triangle_count = 0 for perimeter, count in perimeter_counts.items(): if count == 1: one_integer_sided_right_triangle_count += 1 return one_integer_sided_right_triangle_count def triples_euclids_formula(m, n): """ Generates only primitive triples. https://en.wikipedia.org/wiki/Pythagorean_triple#Generating_a_triple """ assert(m > n) assert(n > 0) assert(gcd(m, n) == 1) assert(m % 2 == 0 or n % 2 == 0) a = m * m - n * n b = 2 * m * n c = m * m + n * n if a < b: return (a, b, c) else: return (b, a, c) def euler_075(): """ The range for m was chosen empirically (aka by trial and error). I tried terminating when L > L_max the first time, but there might be combinations for m = m + 1 that yield lower L. I guess the right criteria to return would be if L > L_max and n == 1. This is good enough for now. I still like my other solution, but it is so much slower because finding the divisors is so expensive. """ L = 0 L_max = 1500000 m = 1 perimeter_counts = {} for m in range(1, 1000): for n in range(1, m): try: L = sum(triples_euclids_formula(m, n)) L_k = L while L_k <= L_max: try: perimeter_counts[L_k] += 1 except KeyError: perimeter_counts[L_k] = 1 L_k += L except AssertionError: pass one_integer_sided_right_triangle_count = 0 for perimeter, count in perimeter_counts.items(): if count == 1: one_integer_sided_right_triangle_count += 1 return one_integer_sided_right_triangle_count if __name__ == "__main__": print("e075.py: " + str(euler_075())) assert(euler_075() == 161667)