1 (1998) · pp.12–14


An Infinite Product for Square-Rooting
with Cubic Convergence


F. L. BAUER
Nordliche Villenstrasse 19
Munich, Germany


In 1913, Friedrich Engel (1861–1941) gave [1] an infinite product for the positive square root √N ,

N  = ( 1 +  1 

 q1

)( 1 +  1 

 q2

)( 1 +  1 

 q3

)( 1 +  1 

 q4

) ...     (N > 0,   N ≠ 1),

where
q1 N – 1

N + 1

 ,     qn+1 = 2qn2– 1.

Oskar Perron (1880–1975), in his book Irrationalzahlen, p. 127 [2], refers to this result and gives Engel's elementary derivation: With the help of the identity

( 1 +  1 

 qn

)( 1 +  1 

 qn+1

)  =  qn+1 – 1

 qn+1

( qn – 1

 qn

) –1 ,

he obtains
( 1 +  1 

 qn

)  · ( 1 +  1 

 qn+1

)  =  ( q1 – 1

 q1

) –1
n=1 n=1

and, thus,
é ù 2
ê ( 1 +  1 

 qn

) ú  =  ( 1 +  1 

 q1

)( q1 – 1

 q1

) –1 = N.
ë n=1 û

Perron, although a specialist on continued fractions, does not mention here that Engel's infinite product is related to the well-known periodic continued fraction, convergent for positive real  z

 z = 1 +  (z – 1)/2
  1 +  (z – 1)/4
  1 +  (z – 1)/4

1 + ...


(see his book [3], p. 320), whose approximants kn are kn = an/bn, where

 an (1 + √ z)n + (1 – √ z)n

2(√ z – 1)n

,      bn (1 + √ z)n – (1 – √ z)n

2√ z(√ z – 1)n

.

There are the recursion equations

 an+2 2

  z – 1

an+1 + an,      bn+2 2

  z – 1

bn+1 + bn; 

the formulas for index-doubling

a2n = an2 + zbn2,
b2n = 2anbn;

and the formulas for index-trebling

a3n = an(an2 + 3zbn2),
b3n = bn(3an2 + zbn2).

From index-doubling immediately follows

 k2n 1

2

( kn  z 

kn

) .

Thus, Newton iteration for the polynomial P(x) ≡ x2 – z,

 ζn+1 1

2

( ζn  z 

ζn

) ,

starting with ζ0 = 1, coincides with index-doubling of the continued fraction approximants, starting with k = 1. The correspondence is

k2n = ζn

[note that a1 = b1 1

  z – 1

 and a2 z + 1

z – 1

, b2 2

z – 1

; thus, ζ1 = k2 z + 1

2

].
 

Moreover, the following identity holds:

an2zbn2= (–1)n.

Thus, for even n,

an2zbn2= 1,       a2n = 2an2 – 1.

There results

z 

k2n

 =  zb2n

a2n

 =  zbn

an

 ·  2an

a2n

2 z 

kn

 ·  a2n + 1

a2n

 =  z 

kn

  ( 1 +  1 

 a2n

)  .

Since the z/k2n = zn converge (from below) toward √ z, the infinite product of Engel with q1 = a2, q2 = a4, q3 = a8, ... is obtained. Explicitly,

 √ z 2z 

 z + 1

 ·  2(z + 1)2

z2 + 6z + 1

 ·  2(z2 + 6z + 1)2

z4 + 28z3 + 70z2 + 28z + 1

 · ... .

L'appétit vient en mangeant.” Is there an even faster way of calculating square roots? Index-trebling should give cubic convergence, corresponding to some cubically converging higher-order Newton-Schröder iteration. In fact, index-trebling amounts to

 k3n  a3n

 b3n

 =  an(an2+ 3zbn2)

bn(3an2+ zbn2)

 = kn  ( 1 +  2(zbn2an2)

3an2+ zbn2

)  .

For odd n,

an2zbn2= –1,       a2n = 2an2 + 1.

Thus,

 k3n = kn  ( 1 +  2 

 3an2+ zbn2

)  = kn  ( 1 +  2 

 4zbn2– 3

)  .

Of course, there is another infinite product for the analytic function √ z; it reads

 √ z 3z + 1 

 z + 3

 ·  3z3 + 27z2 + 33z + 1

z3 + 33z2 + 27z + 3

 ·  3z9 + 297z8 + 6156z7 + 37044z6 + 87642z5 + 87390z4 + 37212z3 + 6084z2 + 315z + 1

z9 + 315z8 + 6084z7 + 37212z6 + 87390z5 + 87642z4 + 37044z3 + 6156z2 + 297z + 3

 · ...

[note that a3 3z + 1

(√ z – 1)3

 and b3 z + 3

(√ z – 1)3

; thus, k3 3z + 1

z + 3

].
 

This corresponds to the cubic Newton-Schröder iteration for ηn = a3n/b3n,

 ηn+1 = ηn  ηn2 + 3z

 3ηn2 + z

 ,

where ηn = k3n, starting with

η0 = k1 = 1.

In order to bring the infinite product to Engel's form, we may write

 ηn+1 = ηn  ( 1 +  2 

 hn+1

)  ,

where

 hn+1 n2+ z

 z – ηn2

 =  4z 

 z – ηn2

 – 3 =  4z 

 zb3ⁿ2a3ⁿ2

 b3n2– 3 = 4zb3n2– 3

and

a3n+1 = a3n(hn+1 + 2),
b3n+1 = b3n·hn+1 = b1h1h2h3·...·hnhn+1.

This iteration finally leads to the infinite product

 z = ( 1 +  2 

h1

)( 1 +  2 

h2

)( 1 +  2 

h3

)( 1 +  2 

h4

) ... ,

where

 hn+1 4z

z – 1 

 h12h22... hn2– 3.

Explicitly, we have

 h1 z + 3

z – 1

,      h2 z3 + 33z2 + 27z + 3

(z – 1)3

,
 h3 z9 + 315z8 + 6084z7 + 37212z6 + 87390z5 + 87642z4 + 37044z3 + 6156z2 + 297z + 3

(z – 1)9

.

Calculating the next hn+1 now needs only one additional multiplication if the previous results are stored. The cubic convergence is rapid. Three examples follow:

 2 = ( 1 +  2

5

)( 1 +  2

197

)( 1 +  2

7761797

)( 1 +  2

467613464999866416197

) ... ,

the first four terms give 62 correct figures;

 3 = ( 1 +  2

3

)( 1 +  2

51

)( 1 +  2

140451

)( 1 +  2

2770663499604051

) ... ,

the first four terms give 46 correct figures;

 5 = ( 1 +  2

2

)( 1 +  2

17

)( 1 +  2

5777

)( 1 +  2

192900153617

) ... ,

the first four terms give 33 correct figures.

There is no need to restrict the iteration to integers. For roots of positive rational numbers, rational arithmetic can be used. If z is close to 1, convergence is dramatic: For example, the first two terms in
 
Ö  3 

2

 =  ( 1 +  2

9

)( 1 +  2

969

)( 1 +  2

912670089

)( 1 +  2

760223786832147978143718729

)  ...

give 8 correct figures, the first three 27 correct figures, the first four 80 correct figures. Furthermore, for 0<z<1, the hn are negative, for example in
 
Ö  1 

2

 =  ( 1 –  2

7

)( 1 –  2

199

)( 1 –  2

7761799

)( 1 –  2

467613464999866416199

)  ... ,
 
Ö  2 

3

 =  ( 1 –  2

11

)( 1 –  2

971

)( 1 –  2

912670091

)( 1 –  2

760223786832147978143718731

)  ... .

What one observes here holds quite generally: The cubic approximants of √1/z are the reciprocals of the approximants of √ z.

The cubic approximation could be used with advantage for numerical calculation of square roots √ z with high-precision arithmetic of floating-point binary numbers, where the mantissa is restricted to the interval [1, 4]. It is, therefore, sufficient to consider the case zÎ[1, 2] and to reduce the case zÎ[2, 4] to the case 4/zÎ[1, 2].

Computationally, in high-precision arithmetic with floating-point numbers, division is the most demanding operation. Now, two Newton steps give quartic convergence, but they need two divisions; the iteration above has the slight advantage that it reaches a cubic order of convergence with only one division per step.

Approximations of still higher order for square-rooting are possible, like one with quintic convergence using the formulas

a5n = an(an4+ 10zan2bn2+ 5z2bn4),
b5n = bn(5an4+ 10zan2bn2+ z2bn4).

Practically, they do not seem to be of importance; the gain in convergence can be made up by repeated use of the simpler cubic iteration.


REFERENCES

  1. Friedrich Engel, in: Verhandlungen der 52. Versammlung deutscher Philologen und Schulmänner in Marburg 1913.
  2. Oskar Perron, Irrationalzahlen, Berlin: de Gruyter (1921).
  3. Oskar Perron, Die Lehre von den Kettenbrüchen, Leipzig: Teubner (1913).
Hosted by uCoz