# Euler Problem 66

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 = 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.

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 routinejust 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 Explicit
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)
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

Posted in Uncategorized

## 6 thoughts on “Euler Problem 66”

1. “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().

2. Michael says:

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

3. David says:

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?

4. Michael says:

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

5. Fernando says:

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.

6. Michael says:

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

Posting code? Use <pre> tags for VBA and <code> tags for inline.