Problem #66 concerns the minimal solutions to Pell Equations. The question reads:
Project Euler Problem 66: Diophantine equation
Consider quadratic Diophantine equations of the form:
x2 – Dy2 = 1
For example, when D=13, the minimal solution in x is 6492 – 13×1802 = 1.
It can be assumed that there are no solutions in positive integers when D is square.
By finding minimal solutions in x for D = {2, 3, 5, 6, 7}, we obtain the following:
32 – 2×22 = 1
22 – 3×12 = 1
92 – 5×42 = 1
52 – 6×22 = 1
82 – 7×32 = 1
Hence, by considering minimal solutions in x for D ≤ 7, the largest x is obtained when D=5.
Find the value of D ≤ 1000 in minimal solutions of x for which the largest value of x is obtained.
My solution for this problem is a bit risky because it relies on arbitrary precision. It uses the useful fact that the solutions to the Pell Equation are related to the continued fraction expansions for the square root of D. Here is my solution:
Solution #1: Arbitrary Precision Approach
We simply use the Python decimal library to obtain very precise continued fractions of the irrational square roots less than 1000. Then we keep calculating continued fractions until one provides a solution to the Pell Equation. Here is an implementation of this approach in Python 2.7:
'''
Author: Walker Kroubalkian
Brute Force Approach to Project Euler Problem #66
'''
import time
from math import sqrt
from decimal import *
getcontext().prec = 74
def isSquare(n):
a = int(sqrt(n))
return a*a==n
def gcd(a,b):
if(a<0 or b<0):
return gcd(abs(a),abs(b))
if(min(a,b)==0):
return max(a,b)
if(a>b):
return gcd(b,a%b)
return gcd(a,b%a)
def continuedFractionDenominators(n,x):
a = Decimal(n).sqrt()
denominators = []
while(len(denominators)<x):
v = int(a)
denominators.append(v)
a = Decimal(1.)/(a-Decimal(v))
return denominators
def evaluateContinuedFraction(n,x):
a = continuedFractionDenominators(n, x)
b = a[::-1]
curNumerator = 1
curDenominator = 0
for x in b:
n = curDenominator + curNumerator*x
d = curNumerator
g = gcd(n,d)
curNumerator = n/g
curDenominator = d/g
return [curNumerator,curDenominator]
def projectEulerProblemSixtySix(n):
maxFound = 3
maxIndex = 2
for c in range(3,n+1):
if(not isSquare(c)):
d = 2
while(True):
e = evaluateContinuedFraction(c, d)
f = e[0]
g = e[1]
if(f*f-c*g*g==1):
if(f>maxFound):
maxFound = f
maxIndex = c
break
d+=1
return maxIndex
start = time.time()
print projectEulerProblemSixtySix(1000)
print ("--- %s seconds ---" % (time.time()-start))
'''
Prints
661
--- 4.85595107079 seconds ---
for input of n = 1000
'''
And with that, we’re done. I needed 74 digits after the decimal place to get an accurate answer, so the program took over 4 seconds to run. I’d like to come back to this problem eventually to look for a more efficient and rigorous method.
Thanks for reading! See you tomorrow.
The Inner Machinations of Walker 