Euler Problem 66

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

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

Leave a Reply

Your email address will not be published. Required fields are marked *