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. *Un*fortunately, 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