Calculating the Convergents of x² −Dy² = ±1 (Part XIII)

A Computer Programming Method for Pellian Equations

The Pellian equation is the Diophantine equation x² − Dy² = z² where z equals 1. The least solutions of the Pell equation are posted in Wikipedia and also listed in Table 91, page 254 of Recreations in the Theory of Numbers by Albert H. Beiler (1966), where the values for D on page 252-253 have been computed using the following two expressions:

The quadratic surd √D, however, can be converted to a continued fraction, then into a series of converging fractions yielding the Pell equation. The convergents p_n∕q_n are calculated from the a_ns, i.e., the coefficients or terms of a continued fraction (see Wikipedia for a good explanative article), and their evaluation requires a series of formulas along with P and Q which are shown on pages 261-262 of the above Theory of Numbers by Beiler. The formulas are defined in the book and are as follows:

In addition, for finding the nth convergent p_n ∕q_n from the previous two convergents and the as, Beiler gives the formulas:

This last equation is important since it determines whether n is the negative Pell solution (n is odd) or the regular Pell solution (n is even). In addition, when Q_n+1 = 1 the Pell x,y values are located at p_n,q_n, respectively.

An example which he gives in the book is that of conversion of √23 into a continued fraction. The result using the formulas from above is shown in Table I:

Thus when Q₅ and Q₉ both equal one then p₄ = 24 and q₄ = 5, and p₈ = 1151 and q₈ = 240, respectively, fulfilling the last formula above. As can be seen only the positive Pell solutions are present since n is even.

Calculation of all five values using Computer Code

The calculations for finding the convergents and generating a table of results using JS computer code is the purpose of this article. Note that the file is a txt file and must first be converted to an html file if one wishes to input their own Ds. A caveat, however: Errors sometime crop up when the numbers are large. Also because JS uses 64 bits for a number, very large number results are often depicted in exponential form.

We can generate a computer program that will give the above results with an even greater number of columns. The way the program is set up is as follows:

• Open the file and input the values for colmn (the column number n) and D.
• The program generates five arrays Pn(), Qn(), an(), pn() and qn() where the appropriate values are to be stored.
• All the equations above are used except for the ones in firebrick red since these formulas are redundant and can be produced simply by:
• Inputting the initial values into the JS arrays. Since JS arrays begin at index 0, initialize Pn(0)=0, Qn(0)=0, an(0)=0, pn(0)=1 and qn(0)=0. These last two values are the trivial values for the equation x² −Dy² = 1 where the curve crosses the x axis at y=0 and are required when calculating for pn(1) and qn(1).
• Likewise for P₁ and Q₁ since these values can be computed from P_n and Q_n when n is set to 1 in the program.
• To obtain the integral part of a₁ and a_n Math.trunc() is performed which removes the decimal point and anything after the decimal point.
• The program calculates the values and stores them into the five arrays using for loops.
• The program tabulates all the data row by row, first Pn() and finally qn() using a series of for loops.

Computer Programming Examples

• Thus, one can run the program and generate the same result as Table I D=23 but with n=22.
• A second uses D=26,n=15. Gives both + and − Pell solutions.
• A third uses D=41,n=18. Gives both + and − Pell solutions.
• Those numbers which produce negative Pell equations are shown in A031396 in the OEIS database.
• A fourth calculates the negative Pell at n=17 and the positive one at n=34 where the p₃₄ is off by 1. The last two digits should be 49 using the values D=313,n=35.

The JS computer code simplififies the calculations. For example, D=211 was given as an example on page 259 in the Bailer Book for us to try to solve. It's easy peasy with the program. It's at n=26 one less than Q₂₇ = 1.
As a final humoristic aside the results for D=1021 by this method gives values at n=49 (the − Pell Solution) for p₄₉ and q₄₉ of 3.15217×10²³ and 9.86500×10²¹, respectively. That's very close to the Avogadro number of 6.02×10²³.

The values calculated by this program can be checked by searching Google books:Canon Pellianus which contains about 1000 entries. It's an old book which probably required lots of manual, tedious calculations in monklike surroundings using the very top two expressions.

This concludes Part XIII. Go back to Part XIIB.

Go back to homepage.