This one has been sticking out like a sore thumb for a while. I couldn’t figure it out for a long while until I found an algorithm for solving Pell numbers. I also used my old favourite prime number sieve. So I can’t really explain why this works. It just does. Better not question it and go on the next one.

#Euler66 from math import sqrt class Algorithms: def pell(self, d): p, k, x1, y, sd = 1, 1, 1, 0, sqrt(d) while k != 1 or y == 0: p = k * (p/k+1) - p p = p - int((p - sd)/k) * k x = (p*x1 + d*y) / abs(k) y = (p*y + x1) / abs(k) k = (p*p - d) / k x1 = x return x def primeGen (self, lowerLimit, upperLimit): isPrime = [False, True] i = 2 while i < upperLimit: isPrime.append(True) i += 1 p = 2 while p * p < upperLimit: if isPrime[p] == True: m = p * p while m < upperLimit: if m % p == 0: isPrime[m] = False m += p p += 1 P = [] count = 0 for p in isPrime: if p and count >= lowerLimit : P.append(count) count += 1 return P class Problem66: def __init__(self, limit): self.limit = limit self.Euler = Algorithms() def Solution(self): P = [] for p in self.Euler.primeGen(2,self.limit): Pell = self.Euler.pell(p) print Pell, p P.append((Pell, p)) return "Answer = "+str(max(P)[1]) if __name__ == '__main__': P = Problem66(1000) print P.Solution()