Euler Problem 66 asks:
Consider quadratic Diophantine equations of the form:
x^2 Dy^2 = 1
For example, when D=13, the minimal solution in x is 649^2 13*180^2 = 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:
3^2 2*2^2 = 1
2^2 3*1^2 = 1
9^2 5*4^2 = 1
5^2 6*2^2 = 1
8^2 7*3^2 = 1Hence, 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.
If you rearrange x^2 – D*y^2 = 1 to D = (x^2 – 1)/y^2, the square root of D is nearly proportional to x/y. This square root ratio can be determined by continued fractions, and as long as D is not a perfect square, has an infinte number of ratios that are solutions. Fortunately, Euler only asks for the first one. The algorithm to do this is described at http://mathworld.wolfram.com/PellEquation.html. Unfortunately, the author wasn’t wasn’t wearing a programmer’s cap when writing, and the equations had to be rearranged. Here is my code that implements Pell’s equation. There’s not much to see in the main routinejust looking for the maximum x and picking off its index. The heavy lifting is in the PellSoln() function. PellSoln() returns an {x,y} pair (Euler only asked for the D associated with the maximum x.) and optionally will return the n-th higher order solution. It’s annotated to the appropriate MathWorld equation. You can see there was some re-ordering. The code runs in a blink.
Option Base 1
Sub Problem_066()
Dim D As Long
Dim x(1000) As Double
Dim Max_X As Double
Dim Answer As Long, T As Single
T = Timer
For D = 1 To 1000
x(D) = PellSoln(D)(1) ‘Option Based 1
If x(D) GT Max_X Then
Max_X = x(D)
Answer = D
End If
Next D
Debug.Print Answer; ” Time:”; Timer – T
End Sub
Function PellSoln(D As Long, Optional nth) As Variant
‘Returns trivial {1,0} if D is a perfect square
‘D = d^2 means x^2 – D*y^2 = x^2-(d^2*y^2) = x^2 -(d*y)^2
‘No other solutions when difference between squares = 1
‘From http://mathworld.wolfram.com/PellEquation.htm
‘Int() function same as Mathworld.Floor() function
Dim x1 As Double, x_nth As Double
Dim y1 As Double, y_nth As Double
Dim n As Long, r As Long
Dim SqrRoot_D As Double
Dim Pn(0 To 100) As Double ‘100 is generous
Dim Qn(0 To 100) As Double
Dim a(0 To 100) As Double
Dim p(0 To 100) As Double
Dim q(0 To 100) As Double
SqrRoot_D = Sqr(D)
If SqrRoot_D = Int(SqrRoot_D) Then ‘D is a perfect square
x1 = 1 ‘trivial solution: 1^2 + D*zero^2 = 1
y1 = 0
PellSoln = Array(x1, y1)
Exit Function
End If
If IsMissing(nth) Then nth = 1
‘n = 0
Pn(0) = 0 ‘Mathworld (15)
Qn(0) = 1 ‘Mathworld (18)
a(0) = Int(SqrRoot_D) ‘Mathworld (8)
p(0) = a(0) ‘Mathworld (9)
q(0) = 1 ‘Mathworld (12)
n = 1
Pn(1) = a(0) ‘Mathworld (16)
Qn(1) = D – a(0) ^ 2 ‘Mathworld (19)
a(1) = Int((a(0) + Pn(1)) / Qn(1)) ‘Mathworld (21)
p(1) = a(0) * a(1) + 1 ‘Mathworld (10)
q(1) = a(1) ‘Mathworld (13)
While a(n) != 2 * a(0)
n = n + 1 ‘ n = 2…
Pn(n) = a(n – 1) * Qn(n – 1) – Pn(n – 1) ‘Mathworld (17)
Qn(n) = (D – Pn(n) * Pn(n)) / Qn(n – 1) ‘Mathworld (20)
a(n) = Int((a(0) + Pn(n)) / Qn(n)) ‘Mathworld (21)
p(n) = a(n) * p(n – 1) + p(n – 2) ‘Mathworld (11)
q(n) = a(n) * q(n – 1) + q(n – 2) ‘Mathworld (14)
Wend
r = n – 1 ‘ a(r+1) = a(n) = 2*a(0)
If r Mod 2 = 0 Then ‘r is even
For n = r + 2 To 2 * r + 1 ‘ r+2 = old n+1
Pn(n) = a(n – 1) * Qn(n – 1) – Pn(n – 1) ‘Mathworld (17)
Qn(n) = (D – Pn(n) * Pn(n)) / Qn(n – 1) ‘Mathworld (20)
a(n) = Int((a(0) + Pn(n)) / Qn(n)) ‘Mathworld (21)
p(n) = a(n) * p(n – 1) + p(n – 2) ‘Mathworld (11)
q(n) = a(n) * q(n – 1) + q(n – 2) ‘Mathworld (14)
Next n
x1 = p(2 * r + 1) ‘Mathworld (28)
y1 = q(2 * r + 1)
Else ‘ r is odd
x1 = p(r) ‘Mathworld (28)
y1 = q(r)
End If
If nth = 1 Then
PellSoln = Array(x1, y1)
Exit Function
End If
x_nth = ((x1 + y1 * SqrRoot_D) ^ nth + (x1 – y1 * SqrRoot_D) ^ nth) / 2 ‘Mathworld (35)
y_nth = ((x1 + y1 * SqrRoot_D) ^ nth – (x1 – y1 * SqrRoot_D) ^ nth) / (2 * SqrRoot_D) ‘Mathworld (36)
PellSoln = Array(x_nth, y_nth)
End Function
Every Diophantine equation has the the “trivial solution” of {1,0}, and so PellSoln() turns trivial if D is a perfect square. The usual accommodations are made for angle brackets. Pell’s equation is big in number theory. I don’t expect to have a use for it again. ;-)
Is there a better way to return a 1×2 array?
…mrt
“Is there a better way to return a 1×2 array?”
Is there a problem with the way you did it?
Actually it creates a 1D array with 2 elements, rather than a 2D array with 1 row and 2 columns.
If I want to return an array from a function I usually decalare and dimension (base 1) an array variable, then assign it to the function:
Function PassArray() as variant
Dim VBAArray(1 to 1, 1 to 2) as double
VBAArray(1,1) = 123
VBAArray(1,2) = 456
PassArray = VBAArray
End Function
It has the advantage that it has the same structure as the array that ends up in the spreadsheet, but I don’t know that it’s really any better than using Array().
Hi Doug –
Thank you. No, no great problem, just a sense that that wasn’t how the pro’s do it. ;-) And that sense turned out to be right.
Your “Actually…” clears up something I’d never understood, being of too literal a mindset.
And contrary to what I thought, I’ll probably have use for Pell again in #64.
Thanks.
…mrt
It is of some interest to write out all of the results and it turns out to be at least as fast to do so. The conditional formatting feature of Excel highlights the maximum value of x (1.64E+37).
Sub Problem_066()
Dim D As Long
Dim T As Single
Dim SolnArray As Variant
Dim tempArray As Variant
ReDim SolnArray(1 To 1000, 1 To 3)
ReDim tempArray(1 To 2)
T = Timer
For D = 1 To 1000
SolnArray(D, 1) = D
tempArray = PellSoln(D)
SolnArray(D, 2) = tempArray(1)
SolnArray(D, 3) = tempArray(2)
Next D
Sheets(1).Range(“A1:C1000?).Value = SolArray
Debug.Print ” Time:”; Timer – T
End Sub
It ran for me in anywhere from 0.012 to 0.03 secs
I have always wanted to know if these three lines can be reduced to one:
tempArray = PellSoln(D)
SolnArray(D, 2) = tempArray(1)
SolnArray(D, 3) = tempArray(2)
If they cannot, is there a way to write the array return from PellSoln to a single row of a 1000 x 2 array in one step?
Hi David –
I played around with it, using Doug’s insights in the first comment, and while they were a different three lines, there were still three.
Don’t take that as the conclusive answer. Certainly PellSoln can be re-written to do the looping there, but you’d want to pass it limits as arguments. Problem 64 goes to 10,000. I’ve discovered that in these tough economic times, 100 just isn’t generous anymore. ;-) For Problem, you’d need to raise the five arrays to (0 to 500), and in order change the returns to PellSoln = Array(x1,y1,n) ; = Array(x1,y1,r +1) ; = Array(x_nth,y_nth, r +1). Then test PellSoln(D)(3) for oddness. R+1 is the n at which the periodicity starts.
…mrt
Hi there!
Thanks for the algorithm it helped me quite a lot.
I translated it into C# and it worked straight away.
If anyone out there is still going through these problems let me know. It will be nice to exchange ideas. I have done about 90 so far so I still have a long way to go.
Hello Fernando –
You’re welcome. Glad to help out. I’ve ripped C# off enough. Nice to know it goes the other way.
I’m at about 140 or so. I’ve slowed down as real life intruded into my hobby. You’ll find the ones that interested me, or taught me something, here.
…mrt