From Pell to Precision: Classical Math Powering Modern Computation

                                                                                                                      By Srimahavishnu Vinjamuri

1.      Introduction:

The Pell equation x²-Dy²=1 has infinitely many non-trivial integer solutions, and the ratios x_k/y_k ​ of higher-order solutions converge to the square root of​ D. We refer to x_k/y_k ​ as the Pell ratio. The companion “negative” equation x²-Dy²=-1 has integer solutions exactly when the continued-fraction expansion of the square root of D​ has odd period; in particular, this holds for primes D ≡ 1(mod4).

In this article we turn these classical results into a computational tool. We compute square roots of selected primes to ~15 correct decimal digits (IEEE-754 float64) using Pell ratios, and when fundamental Pell solutions become astronomically large, we use continued-fraction convergents—best rational approximations that remain numerically tractable. We benchmark both approaches against Python’s math.sqrt(D) (float64) and mpmath.sqrt(D) (high precision), and also report, across a range of primes, the smallest solutions/convergents whose ratio x/y reaches 32-digit accuracy in selected cases.

Speed at a glance: In our benchmarks, evaluating a precomputed ratio p/q is often faster than taking a square root: about 1.5–2.5× faster than math.sqrt(D) in float64, and ~25–50% faster than mpmath.sqrt(D) at higher precision. These wins are most pronounced when the same D is evaluated many times (so (p,q) are reused/cached) or when high precision is required; for single, one-off float64 evaluations, math.sqrt(D) can still be the quicker choice.

2.      Solution to the Pell equation:

The Pell equation is:

                       

                                       - Equation 1

If the initial non-trivial solution is x₁,y₁, then the solution of the order k+1 is

   

                        -Equation 2

                 

                           -Equation 3

The negative Pell equation  is:

                  
      - Equation 4

If x_pos and y_pos are the initial non-trivial solutions to Equation 1, then the solution of the order m+1 for Equation 4 is:

        

            -Equation 5

                 

                -Equation 6

By applying the recurrence relations from Equations 2, 3, 5, and 6, we generated solutions to the Pell equations 1 and 4 for the first ten primes, as well as for the 100th and 1000th primes. For each case, the search continued until the smallest integer pair (x,y) was identified such that the ratio x/y approximates the square root of D​ with at least 15 correct decimal digits. This 15-digit threshold corresponds to the effective precision of IEEE-754 binary64 (float64) arithmetic, since Python’s math.sqrt(D) relies on hardware-accelerated 64-bit doubles with a 53-bit significand, yielding approximately 15–16 accurate decimal digits. The resulting solutions are presented in the following sections.

2.1  Square Root of 2:

For D = 2 in Equation 1, the Pell equation becomes:

The solutions that give us a minimum of 15-digit accuracy for the square root of 2 are:

131836323 and 93222358.  The solution obtained using Pell ratio is nearly twice as fast as math.sqrt(2).

2.2 Square Root of 3:

For D = 3 in Equation 1, the Pell equation becomes:

The solutions that give us a minimum of 15-digit accuracy for the square root of 3 are:

189750626 and 109552575.  The solution obtained using Pell ratio is nearly twice as fast as math.sqrt(3).

2.3 Square Root of 5:

For D=5 in Equation 1, the Pell equation becomes:

Since D is a prime of the form 4n+1, its negative Pell equation will also have integer solutions.

For D=5 in Equation 4, the negative Pell equation becomes:

The solutions that give us a minimum of 15-digit accuracy for the square root of 5 are:

299537289 and 133957148, which occur for the positive Pell equation.

 The solution obtained using Pell ratio is around 1.75 times faster than math.sqrt(5).

2.4 Square Roots of Prime Numbers up to the 1000^th Prime:

In addition to √2, √3, and √5, the square roots of the remaining primes among the first ten primes, as well as the 100th prime (541) and the 1000th prime (7919), are computed using ratios derived from the corresponding Pell equation solutions. These values are compared with results obtained from Python’s math.sqrt(D) function in terms of both accuracy and speed, and further benchmarked against the more precise but slower mpmath.sqrt(D) function. The results are tabulated below, with comparisons against mpmath.sqrt(D) highlighted for those Pell-based ratios that yield highly accurate square root values. 

Table 1: Pell’s Solutions for the First Ten Primes, 100^th Prime and 1000^th Prime using float64

 

Note: “Digits of accuracy” mentioned here and elsewhere in this article are counted from the rounded printouts shown (mpmath values printed to ~40 significant digits; float64 to 17). Because this is a visual count, values can differ by about ±1 from the true high-precision agreement.

2.5 Square Roots of Larger Prime Numbers:

From Table 1, the Pell-based approach scales for finding out the square roots of non-square numbers and we have seen it successfully producing accurate results even for the 1000^th prime number.

We now examine it for larger primes.

For the 10,000th prime, 104729, the initial non-trivial Pell solutions involve extremely large numbers, making their ratio difficult to compute and handle. In such cases, approximate solutions can instead be obtained using the Continued-Fraction method, which yields convergents that approach √D with high accuracy while avoiding the need to calculate the enormous fundamental Pell solutions.

Below is the approach to be followed for finding out the convergents.

2.5.1 Continued-Fraction Expansion of √D:

For a non-square positive integer D:

 √D =[a₀​;a₁​,a₂​,…,a_i,…,a_L​​],

where a₀=integer portion of √D and the rest of the a_i repeat with period L.

2.5.2 Recurrence Relations for Convergents

Define numerators p_k ​ and denominators q_k recursively such that: 

p_-2=0 and p_-1=1

q_-2=1 and q_-1=0

 For k≥0, the recurrence equations are:

       

                    - Equation 7

     

                    - Equation 8

 Every convergent satisfies the following:

 

 

And

 

Each convergent p_k/q_k ​ is a best rational approximation to √D. No fraction with denominator<q_k is closer. Consequently, q≥10^k/2 is a practical threshold for k decimal digits.

 In Equations 7 and 8 above,

         
                     -Equation 9

Where

 
              -Equation 10

                                                                        

              -Equation 11

The initial values are:

 

Using Equations 7 to 11, we generate continued-fraction convergents and stop at the first one whose denominator meets the threshold q≥10^k/2 . We run this for k=16 and k=32. Because IEEE-754 binary64 (float64) provides only ~15–16 correct digits, the k=32 case is evaluated in high-precision arithmetic (mpmath), with p/q and mpmath.sqrt(D) computed at matching precision. 

Let us discuss the results in the subsequent sections.

2.5.3 Square Root of 10000^th Prime 104729

The mpmath based square root of 104729 is: 323.6186026791414226842975608589216364423

The convergents to the equation for k=32 are,

p = 9771686109987793705 and q = 30195069223743422.

The mpmath based ratio of convergents is:     323.6186026791414226842975608589216360847 

The float64-based square root of 104729 is:   323.61860267914142

The convergents to the equation for k=16 are,

 p = 133316627027, q = 411956006

The float64-based ratio of convergents is:       323.61860267914142 

The ratio agrees to 33 digits. 

The computational speeds are as given below.

Table 2 Computation Time for Calculating Square Root of 104729

Time Taken (microseconds)				Efficiency of ratio-based computation
mpmath.sqrt(D)	mpmath convergent ratio	float64-based math.sqrt(D)	float64 convergent ratio	mpmath vs mpmath	p/q vs math.sqrt(D) (float64)
2.67	2.04	0.08	0.05	1.31	1.70

 2.5.4 Square Root of 100000^th Prime 1299709

The mpmath based square root of 1299709 is: 1140.047806015168821788763613850391783587

The convergents to the equation for k=32 are,

p = 60561156670261640493 and q = 53121593981170162.

The mpmath based ratio of convergents is:        1140.047806015168821788763613850391783596 

The float64-based square root of 1299709 is:   1140.0478060151688

The convergents to the equation for k=16 are,

 p = 502728132791, q = 440971098

The float64-based ratio of convergents is:       1140.0478060151688

 The ratio agrees to  34 digits. 

The computational speeds are as given below. 

Table 3 Computation Time for Calculating Square Root of 1299709 

Time Taken (microseconds)				Efficiency of ratio-based computation
mpmath.sqrt(D)	mpmath convergent ratio	float64-based math.sqrt(D)	float64 convergent ratio	mpmath vs mpmath	p/q vs math.sqrt(D) (float64)
3.19	2.01	0.13	0.05	1.58	2.47

 2.5.5 Square Root of 1000000^th Prime 15485863

The mpmath based square root of 15485863 is: 3935.208126643367258027724271285199066394

The convergents to the equation for k=32 are,

p = 60119015904572548365 and q = 15277213801612202.

The mpmath based ratio of convergents is:        3935.208126643367258027724271285199063875 

The float64-based square root of 15485863 is:   3935.2081266433675

The convergents to the equation for k=16 are,

 p = 562762556485, q = 143007063

The float64-based ratio of convergents is:       3935.2081266433675 

The ratio agrees to  32 digits. 

The float64 ratio matches math.sqrt(D) and is accurate to ~12–13 digits (float64 limit). This is due to the limitation in the float64-based computation 

The computational speeds are as given below.

 Table 4 Computation Time for Calculating Square Root of 15485863 

Time Taken (microseconds)				Efficiency of ratio-based computation
mpmath.sqrt(D)	mpmath convergent ratio	float64-based math.sqrt(D)	float64 convergent ratio	mpmath vs mpmath	p/q vs math.sqrt(D) (float64)
5.71	4.31	0.17	0.10	1.32	1.58

 2.5.6 Square Root of the Largest 32-digit Prime 

The largest 32-digit prime is 99999999999999999999999999999979.

 The mpmath based square root of this number is: 9999999999999999.99999999999999895

The convergents to the equation for k=32 are,

p = 199999999999999999999999999999979 and q = 20000000000000000.

The mpmath based ratio of convergents is:              9999999999999999.99999999999999895

The float64-based square root  is:                              10000000000000000

The convergents to the equation for k=16 are,

 p = 9523809523809519999999999999999, q = 952380952380952

The float64-based ratio of convergents is:                10000000000000000

 From the above results, the mpmath-precision convergent ratio p/q agrees with the mpmath reference square root of D​ to all reported digits. 

Float64 rounding near 10¹⁶: float64 both p/q and math.sqrt(D) print as 10,000,000,000,000,000. Around 1016, the binary64 spacing (ULP) is 2, so the nearest representable numbers are 9,999,999,999,999,998 and 10,000,000,000,000,000. Since the square root of D=9,999,999,999,999,999.99999999999999895 lies within 1 of 10¹⁶, both values round to 10¹⁶ exactly. The agreement in float64 is therefore a rounding artifact; the true square root is slightly smaller than the square root of 10³². 

The computational speeds are as given below 

Table 5 Computation Time for Calculating Square Root of the Largest 32-digit Prime 

Time Taken (microseconds)				Efficiency of ratio-based computation
mpmath.sqrt(D)	mpmath convergent ratio	float64-based math.sqrt(D)	float64 convergent ratio	mpmath vs mpmath	p/q vs math.sqrt(D) (float64)
2.84	2.00	0.09	0.05	1.42	1.79

2.6 Square Roots of Primes up to 1000^th Prime with 32-digit Accuracy 

In Section 2.5, we examined convergents of Pell equation whose ratios approximate square roots of prime numbers with a minimum accuracy of 32 digits for selected cases. In this section, we extend the analysis to identify the smallest convergents (or solutions) of Pell equation corresponding to the first ten primes, as well as the 100^th and 1000^th  primes. For reference, accuracies of 32 digits and beyond have been obtained directly using Python’s mpmath.sqrt(D) function.

The table below lists the solutions/convergents up to the desired level of accuracy for the first 10 primes, and the 100^th, 1000^th, 10,000^th, 100,000^th and 1,000,000^th  primes. 

Table 6: Pell’s Solutions for the First Ten Primes, and larger primes using mpmath  

Notes:

• The numbers 541, 7,919, 104,729, 1,299,709, and 15,485,863 correspond to the 100th, 1,000th, 10,000th, 100,000th, and 1,000,000th prime numbers, respectively.
• Section 2.5.6 presents the largest 32-digit prime together with its smallest convergents, which approximate its square root to 32-digit accuracy.
•  For each D, we report the first convergent whose denominator satisfies q≥10^k/2.

3.      Disclaimers:

•  The reported computation times for the division of Pell’s solutions/convergents do not include the one-time cost of generating these solutions. They represent only the time required to evaluate the ratio.
• The computational speeds presented in this article may vary depending on the system and runtime conditions. The overall trend, however, is expected to remain consistent. 

4.      Conclusions and inferences:

•  We present the smallest solutions/convergents to the Pell equation whose ratio x/y approximates the square root of D​ for the first ten primes and the 100th, 1,000th, 10,000th, 100,000th, and 1,000,000th primes, as well as the largest 32-digit prime.

• Accuracy: Solutions derived from Pell equations and continued-fraction convergents yield highly accurate square root approximations, achieving up to 15–16 correct digits in float64 arithmetic and extending to 33–36 digits under arbitrary-precision computation (mpmath).

•  Computational Efficiency: Ratio-based methods demonstrate clear performance advantages, running 1.5–2.5 times faster than Python’s math.sqrt(D) in float64 environments, and 25–50% faster than mpmath.sqrt(D) when higher precision is required. The division p/q is competitive when many evaluations reuse the same (p,q). A library of convergents to the desired level of accuracy can be prepared for the repeatedly used square roots of numbers like 2,3,5 etc.

•  Limitations: For a single use (or a few uses), generating a convergent to hit 15–32 digits can be costly (big-integer CF steps, possibly huge q). In that regime, math.sqrt(D) (float64) will usually win.

•  Scalability: The approach remains stable and effective across a wide range of cases, from small primes to extremely large numbers, including the largest 32-digit prime tested.

•  Precision Boundaries: While float64 arithmetic is inherently limited to approximately 15–16 digits of precision, the use of continued-fraction convergents with high-precision libraries enables far greater accuracy.

•  Significance: The study establishes Pell equation–based methods as a stable, scalable, and often faster evaluation path - especially when ratios are reused or when high precision is required.

•  Applications in Cryptography and Data Security: High-precision, efficient square root computation has potential relevance in public-key cryptosystems and primality testing,  where large prime computations are central.

•   Applications in Computational Mathematics and Simulation: Ratio-based square root evaluations may also be useful in scientific computing, high-precision simulations, and numerical algorithms where stability, scalability, and efficiency are critical.

•  Applications in Integer-only Environments: In integer-only settings, a precomputed library of Pell’s solutions/convergents for each prime or non-square number enables rapid division-based calculations, achieving the desired level of accuracy for square root estimation.