from lib_prime import prime_nth def varray(i): return i * 2 + 1 def rem(i): n = varray(i) p = prime_nth(n) return (pow(p - 1, n) + pow(p + 1, n)) % (p * p) def bin_search(lo, hi, target): if hi - lo == 1: return varray(hi) i = lo + ((hi - lo) // 2) r = rem(i) # print(f"{i=:<6} {r=}") if r < target: return bin_search(i, hi, target) elif r > target: return bin_search(lo, i, target) assert False def euler_123(): """ I found out that when n is even, the remainder is always 2. When the remainder is odd, the series grows monotonically. That means we can use a binary search on odd values to find the solution. """ return bin_search(1000, 30000, 10**10) if __name__ == "__main__": solution = euler_123() print("e123.py: " + str(solution)) assert(solution == 21035)