Если Вы найдете статью полезной и интересной - не сочтите за труд, переведите материал или хотя бы часть его и отправьте на адрес algolist@manual.ru. Задать вопросы или просто написать письмо можно также на странице контактов.

There are certainly people who regard Ц2 as something perfectly obvious but jib at Ц-1. This is because they think they can visualise the former as something in physical space but not the latter. Actually Ц-1 is a much simpler concept.

Edward Charles Titchmarsh (1899-1963).

1  Introduction

The constant Ц2 is famous because it's probably one of the first irrational numbers discovered. According to the Greek philosopher Aristotle (384-322 BC), it was the Pythagoreans around 430 BC who first demonstrated the irrationality of the diagonal of the unit square and this discover was terrible for them because all their system was based on integers and fractions of integers. Later, about 2300 years ago, in Book X of the impressive Elements, Euclid (325-265 BC) showed the irrationality of every nonsquare integer (consult [7] for an introduction to early Greek Mathematics). This number was also studied by the ancient Babylonians. The history of the famous sign Ц goes back up to 1525 in a treatise named Coss where the German mathematician Christoff Rudolff (1499-1545) used a similar sign to represent square roots.

Definition 1 Ц 2 is the unique positive root of the polynomial equation :

x2-2=0.

Theorem 2 Ц 2 is an irrational and algebraic number.

Proof. Suppose that Ц 2=p/q where p and q are relatively primes, then p²=2q² therefore p is even and p=2p^ў which leads to q²=2p^ў2 and to the fact that q=2q^ў. This is in contradiction with p and q being relatively primes.  

We will now introduce some of the techniques available to compute this number.

2  Sequences with integers

2.1  Elementary approach

An elementary and ancient algorithm consists in using the double sequence

м н о pn+1=pn+2qn qn+1=pn+qn

(1)

which verifies

pn+1 qn+1 =  pn/qn+2 pn/qn+1

so that :

lim n®Ґ  xn= lim n®Ґ   pn qn =Ц2.

It is of interest to note that this quite simple recursion only uses additions with integers and eventually a final division to convert the fraction into the usual decimal representation. Denoting the error by e_n=|Ц2-p_n/q_n|, one easily see that e_n+1 < e_n/5, which involves a geometrical convergence (that is, one digit will be added at every N iterations).

Starting with p₀=q₀=1 we obtain :

к к к к к к к x1=1.(500ј), x3=1.41(666...), x6=1.4142(011...), x9=1.414213(624...),

x₃=17/12 was used by Mesopotamians to replace Ц2. It may be convenient to use a matrix representation of the sequence, we write the sequence like this :

ж з и pn+1 qn+1 ц ч ш = ж з и 1 2 1 1 ц ч ш ж з и pn qn ц ч ш

so that

ж з и pn+1 qn+1 ц ч ш =An ж з и p1 q1 ц ч ш = ж з и 1 2 1 1 ц ч ш n   ж з и p1 q1 ц ч ш .

2.1.1  Quadratic convergence

This matrix representation may be used with n=2^p and thanks to the relation (exercise : show this by recurrence) :

An= ж з и an 2bn bn an ц ч ш

comes the following algorithm

м н о ap+1=ap2+2bp2 bp+1=2apbp

starting with a₀=b₀=1. It's easy to see that this new process has a quadratic convergence rate and its first iterates are given here :

(3,2),(17,12),(577,408),(665857,470832),(886731088897,627013566048),...

and

к к  886731088897 627013566048 -Ц2 к к < 10-24.

Note that, in each step of this algorithm, the product of large integers is required and this may be done very efficiently using fast multiplication methods. This method is related to Newton's iteration by observing that :

xp+1=  ap+1 bp+1 =  ap2+2bp2 2apbp =  1 2  ap bp +  bp ap =  xp 2 +  1 xp ,

which will be derived by different means later in this essay (see section 6).

2.1.2  Cubic convergence and quintic convergence

The same approach leads to the Cubic sequence

м н о ap+1=ap( ap2+6bp2) bp+1=bp( 3ap2+2bp2)

starting with a₀=b₀=1 the first iterates are given by :

(7,5),(1393,985),(10812186007,7645370045),...

and also to the Quintic sequence

м н о ap+1=ap( ap4+20ap2bp2+20bp4) bp+1=bp( 5ap4+20ap2bp2+4bp4)

giving the iterates

(41,29),(1855077841,1311738121),...

We have generated algorithms of this nature at any order of convergence using in the sequence polynomials of higher degrees.

2.2  Modification of the sequence

If we set u_n=p_n+q_n and v_n=q_n then the double sequence (1) is equivalent to this new one :

м н о un+1=2un+un-1 vn+1=un

and u_n/v_n=u_n/u_n-1=(p_n+q_n)/q_n converge to Ц2+1. So starting with u₀=1,u₁=2 gives the set of u_n:

(1,2,5,12,29,70,169,408,985,2378,5741,13860,33461,80782,195025,470832,...)

those numbers are sometime called Pell numbers and the sequence u_n is a Pell sequence. The first 41 numbers may be found in [10] as well as other similar sequences like Lucas sequences used in Primality tests. For example

u15 u14 -1=  470832 195025 -1=1.4142135623(637...)

and it was found with very little efforts.

This method is probably due to ancient Babylonians.

2.3  Improvement of the sequence

It's possible to improve this method by using a more general sequence :

м н о pn+1=(A+B)pn+2Aqn qn+1=2Bpn+(A+B)qn

Easy manipulations show that

lim n®Ґ   pn qn =   ж Ц  A B

and that for the error

en= к к к   ж Ц  A B   -xn к к к ,

we have the bound

en+1 < en ж з и   ж Ц  A B   -1 ц ч ш 2   .

The more A/B is close to 1 the more the speed of convergence is increased.

Example 3 To illustrate the method let's use the relation

Ц 2=  7 5   ж Ц  50 49

which corresponds to A=50 and B=49 and the iterative system becomes :

p1=q1=1,       м н о pn+1=99pn+100qn qn+1=98pn+99qn

as for the error, we have

en+1 < en ж з и   ж Ц 1+  1 49   -1 ц ч ш 2   < en/9604.

Here are the first iterates :

к к к к к к к x1=1.4, x2=1.414213(197...), x3=1.4142135623(637...), x4=1.41421356237309(481...)

and almost 4 new digits are added with each step.

3  Continued fraction

3.1  Elementary continued fraction

The idea is to write Ц2=1+1/a₁ with a₁=2.4142135...=2+1/a₂ and so on... This gives the well-known development :

Theorem 4 (Bombelli 1572, [4]).

Ц 2=1+  1 2+1/2+...

It is more convenient to use the notation

Ц2=[1;2,2,2,2,...].

The development is periodic of period 1 with number {2}. It's a general result : square roots can always be represented with a periodic continued fraction development [5]. It's possible to show that the previous sequence also gives the continued development of Ц2. We give the first rational approximations :

ж и  pk qk ц ш k=1...  = ж и  3 2 ,  7 5 ,  17 12 ,  41 29 ,  99 70 ,  239 169 ,  577 408 ,  1393 985 ,  3363 2378 ,  8119 5741 ,  19601 13860 ,  47321 33461 ,ј ц ш .

3.2  Other continued fractions

Faster converging continued fractions may be obtained, for example :

5Ц2-7=[0;14,14,14,14,...].

The development is periodic of period 1 with number {14}. This time, the first rational approximations are :

ж и  pk qk ц ш k=1...  = ж и  1 14 ,  14 197 ,  197 2772 ,  2772 39005 ,  39005 548842 ,  548842 7722793 ,  7722793 108667944 ,ј ц ш

and of course Ц2 » 7/5+p_k/(5q_k), the error is about 1/q_k². Observe that here p_k=q_k-1 which simplifies the evaluation of the approximations.

Other fast converging continued fractions are given by

29Ц2-41=[0;82,82,82,82,...]

and

169Ц2-239=[0;478,478,478,478,...].

It's possible to find interesting continued fractions for numbers of the form mЦ2-n and it's not hard to show that if the integers m and n are such as

n2-2m2=-1

(this is a nonlinear Diophantine equation and it's called a Pell's equation) then

mЦ2-n=[0;2n,2n,2n,2n,...].

From a fundamental result on Pell's equations the solutions are contained in the convergents of the simple continued fraction of Ц2 (there are exactly the convergents of even rank [5] in the case of Ц2). It's also possible to show that there are given by the sequence

м н о mk+1=3mk+2nk nk+1=4mk+3nk

starting with (1,1).

Hence the first couples (m_k,n_k) are

(1,1),(5,7),(29,41),(169,239),(985,1393),(5741,8119),(33461,47321),...

Example 5 With the couple (33461,47321) we have the approximation

33461Ц 2-47321 »  1 2.47321

that is

Ц 2 »  47321 33461 +  1 3166815962 =1.4142135623730950488(369...)

(19 correct digits).

4  Infinite products

The two well-known Infinite products for cosx and sinx are :

Theorem 6 (Euler 1748, [2])

cosx= ж и 1-  4x2 p2 ц ш ж и 1-  4x2 9p2 ц ш ж и 1-  4x2 25p2 ц ш ј sinx=x ж и 1-  x2 p2 ц ш ж и 1-  x2 4p2 ц ш ж и 1-  x2 9p2 ц ш ј

Setting x=p/4 in the cosine product leads to

1 Ц2 = ж и 1-  1 4 ц ш ж и 1-  1 36 ц ш ж и 1-  1 100 ц ш ...

therefore

Ц2= ж и  2.2 1.3 ц ш ж и  6.6 5.7 ц ш ж и  10.10 9.11 ц ш ж и  14.14 13.15 ц ш ...= Х k і 0   (4k+2)2 (4k+1)(4k+3) ,

(2)

or, which is equivalent, to the aesthetic formula :

Ц2= Х k і 0  ж и 1+  1 4k+1 ц ш ж и 1-  1 4k+3 ц ш = ж и 1+  1 1 ц ш ж и 1-  1 3 ц ш ж и 1+  1 5 ц ш ж и 1-  1 7 ц ш ...

(3)

The convergence is extremely slow as we may observe with some iterates :

к к к к к к к x1=1.(333...), x10=1.4(054...), x100=1.41(332...), x1000=1.414(125...).

Euler gave the interesting product (2) in 1748 [2], and it's very similar to John Wallis formula for p:

p 4 = ж и  2.4 3.3 ц ш ж и  4.6 5.5 ц ш ж и  6.8 7.7 ц ш ж и  8.10 9.9 ц ш ...,

while the product (3) can be compared to the celebrated series

p 4 =1-  1 3 +  1 5 -  1 7 +  1 9 -...

With x=p/8 in the infinite product we extract the also slow converging infinite product

2+Ц2=4. ж и  3.5 4.4 ц ш 2   ж и  11.13 12.12 ц ш 2   ж и  19.21 20.20 ц ш 2   ж и  27.29 28.28 ц ш 2   ...

5  Taylor's expansions

From the Taylor expansion of (1+x)^1/2 or (1-x)^-1/2, it's possible to find another class of algorithms based on series computation.

Theorem 7 (Newton 1665). Let - 1 < x < 1, then

Ц   1+x   =1+  1 2 x-  1 2.4 x2+  1.3 2.4.6 x3-ј 1/ Ц   1-x   =1+  1 2 x+  1.3 2.4 x2+  1.3.5 2.4.6 x3+ј

(4)

This theorem was first stated by Isaac Newton (1643-1727) and a correct proof appears later and is due to Leonhard Euler (1707-1783).

Applying the first expansion in (4) which is still valid for x=1 produces the very slowly convergent and alternating series

Ц2=1+  1 2 -  1 2.4 +  1.3 2.4.6 -  1.3.5 2.4.6.8 +...,

which first terms are :

к к к к к к к x1=1.(500...), x2=1.(375...), x3=1.4(375...), x4=1.(398...).

A nice improvement is to use the previous results on continuous fractions. We can write :

Ц2=  pk qk   ж Ц  2qk2 pk2 ,

and for each value of k, this will provide a formula where the number inside the square root of the right hand side of the formula is near 1. When k increases this number tends to 1. Using this remark leads to a sequence of formulae for Ц2, which are more and more efficient :

Ц2= ж з и  3 2   ж Ц  8 9   ,  7 5   ж Ц  50 49   ,  17 12   ж Ц  288 289   ,  41 29   ж Ц  1682 1681   ,  99 70   ж Ц  9800 9801   ,  239 169   ж Ц  57122 57121   ,ј ц ч ш .

Example 8 Using this last formula jointly with Newton's series expansion gives the very fast converging and easy to implement sequence

Ц 2=  239 169 ж и 1+  1 2  1 57121 -  1 8  1 571212 +  3 48  1 571213 -ј ц ш

for which the first terms are :

к к к к к к к x1=1.4142(011...), x2=1.414213562(427...), x3=1.41421356237309(457...), x4=1.41421356237309504880(687...).

With a very little calculation we have computed 20 digits of Ц 2.

Example 9 (Euler 1755, [4], [8]). The relation Ц 2=(7/5)/Ц (49/50) produces the nice and easy to compute by hand series expansion :

Ц 2=  7 5 ж и 1+  1 100 +  1.3 100.200 +  1.3.5 100.200.300 +... ц ш .

The main advantage of using Taylor's formulae is to require only basic multi-precision operations between numbers.

6   Newton's iteration

A very efficient way to compute square roots is to use the Newton iteration [9] on the function f(x)=x²-2. It gives the following sequence :

x0=1,

xn+1=xn-  f(xn) fў(xn) =  1 2 ж и xn+  2 xn ц ш =  xn 2 +  1 xn ,

x_n+1 can also be interpreted as the simple mean between the approximation x_n and 2/x_n which is also another approximation.

The first terms of this sequence are :

к к к к к к к x1=1.(500...), x2=1.41(666...), x3=1.41421(568...), x4=1.41421356237(468...).

Number D of correct digits after n iterations with the elementary Newton iteration :

n 1 2 3 4 5 6 7 8 9 10 D 0 2 5 11 24 48 97 195 391 783

It's a well-known result that the rate of convergence of Newton's method and therefore of this sequence is quadratic (the number of digits is doubling at each iteration) for a starting point x₀ sufficiently close to the root of the equation.

Computing Ц2 with this algorithm is convenient if one can compute easily the inverse of x_n. In fact, it is simpler to use only multiplications. This can be reached by using a modified sequence which converges to 1/Ц2 : it consists in using the Newton iteration but this time with the function f(x)=1/x²-2, yielding the sequence

x0=  1 2 ,       xn+1=xn ж и  3 2 -xn2 ц ш =xn  +xn ж и  1 2 -xn2 ц ш .

Observe that in the right hand side relation the increment tends to zero. We give the first iterates :

к к к к к к к к к x1=0.(625...), x2=0.(693...), x3=0.70(670...), x4=0.707106(444...), x5=0.707106781186(307...),

again the convergence is quadratic. A final multiplication by 2 will give the value of Ц2. The advantage of this method is to avoid a multi-precision division at each step and replace it by two multi-precision multiplications. This algorithm is extremely efficient and may be used to compute Ц2 up to billion's of digits.

7  Cubic iteration

7.1  Halley's iteration

Using the second derivative of f(x)=x²-2, it's possible to find a sequence with cubic convergence (the number of digits is multiplied by 3 at each step). The general formula is given by Halley's iteration (the original form was introduced by the astronomer Edmund Halley (1656-1742) in 1694) :

xn+1=xn-  f(xn) fў(xn) ж и 1-  f(xn)fўў(xn) 2fў(xn)2 ц ш -1   .

Applying this formula with function f(x) gives :

x0=1,       xn+1=xn  xn2+6 3xn2+2 =xn ж и 1+2  (2-xn2) 3xn2+2 ц ш .

The first terms of this sequence are :

к к к к к x1=1.4(000...), x2=1.414213(197...), x3=1.414213562373095048(795...).

Number D of correct digits after n iterations with Halley's iteration :

n 1 2 3 4 5 6 7 8 9 10 D 1 6 20 61 185 557 1673 5022 15067 45204

7.2  Another cubic iteration

In [1], an interesting cubic iteration is given which converge to 1/ЦA

xn+1=  xn 8 (15-10Axn2+3A2xn4).

It may be deduced from the general Householder's iteration [6]

xn+1=xn-  f(xn) fў(xn) ж и 1+  f(xn)fўў(xn) 2fў(xn)2 ц ш ,

which has cubical convergence for a close enough starting point x₀ (observe the similarity with Halley's iteration).

By applying this general pattern for f(x)=1/x²-A we obtain :

xn+1=xn+  xn 8 ( 7-10Axn2+3A2xn4) = xn+  3A2xn 8 ж и xn2-  1 A ц ш ж и xn2-  7 3A ц ш .

If we explicit this formula with A=2, and starting with x₀=1/2,

xn+1=xn+  xn 8 (7-20xn2+12xn4)=xn+  xn 8 (2xn2-1)(6xn2-7),

and, in this iteration, no multiprecision division is required. The first iterates are :

к к к к к x1=0.(671...), x2=0.70(689...), x3=0.7071067811(398...).

Number D of correct digits after n iterations with Householder's cubic algorithm :

n 1 2 3 4 5 6 7 8 9 10 D 0 2 10 30 90 269 808 2425 7276 21829

This method may also be very efficient for high precision computation.

8  Quartic and high order iterations

The direct application of the modified iterations on the function f(x)=1/x²-1/2 produces a set of high order algorithms. It's interesting to note that relatively easy algorithms of any order may be derived.

In the following lines we will use the notation

hn=1-  xn2 2 .

8.1  Quartic iteration

The quartic modified iteration is

xn+1=xn+  xn 16 ( 8hn+6hn2+5hn3) .

For example, if we take the initial value x₀=3/2, the first two iterations are

к к к x1=1.414(123...), x2=1.41421356237309(494...).

Number D of correct digits after n iterations with the quartic algorithm :

n 1 2 3 4 5 6 7 D 3 14 63 254 1019 4078 16312

8.2  Quintic iteration

The same method also produces the quintic (fifth-order) modified iteration

xn+1=xn+  xn 128 ( 64hn+hn2( 48+40hn+35hn2) ) ,

and again with the initial value x₀=3/2, the first two iterations are :

к к к x1=1.4142(236...), x2=1.414213562373095048801688(932...).

Number D of correct digits after n iterations with the quintic algorithm :

n 1 2 3 4 5 6 D 4 24 123 615 3076 15380

Algorithms with any order of convergence may also be generated easily (for example : sextic, septic, octic, ... iterations are available !) and we end this section with the octic (eighth-order) iteration :

м п н п о dn=1024hn+hn2( 768+640hn+560hn2+hn2(504hn+hn2( 462+429hn) )) , xn+1=xn+  1 2048 xndn.

9  Double Iteration

9.1  Quadratic double iterative procedure

During the first years of computer time (EDSAC group at Cambridge University 1951), the following double sequence (easily deduced from Newton's iteration on the inverse of the square root) was proposed to compute square roots :

x0=A,h0=A-1,

and

м п п н п п о xn+1=xn ж и 1-  1 2 hn ц ш hn+1=  1 4 hn2( hn-3)

then the sequence h_n tends to 0 and the sequence x_n tends to ЦA (A < 3).

Example 10 Let A=2, therefore

к к к к к к к к к x1=1.(000...),h1=-0.5 x2=1.(250...),h2=-0.21875... x3=1.(386...),h3=-0.03850... x4=1.41(341...),h4=-0.001126... x5=1.41421(288...),h5=-0.000000951...

and then the convergence becomes quadratic. The number of non quadratic cycles is reduced when A tends to 1. For example, using relation  Ц 2=7/5Ц{50/49} so that A=50/49, will lead to a quadratic converge at once.

9.2  Other rates of convergence

Small modifications in the previous procedure may increase the order of convergence, for example the following algorithm has cubic convergence :

x0=A,h0=A-1,

and

м п п н п п о xn+1=xn ж и 1-  1 2 hn+  3 8 hn2 ц ш hn+1=  1 64 hn3( 40-15hn+9hn2) .

Algorithms at any order of convergence may also be generated.

References

[1]

J.M. Borwein and P.B. Borwein, Pi and the AGM - A study in Analytic Number Theory and Computational Complexity, A Wiley-Interscience Publication, New York, (1987)

[2]

L. Euler, Introduction à l'analyse infinitésimale (french traduction by Labey), Barrois, ainé, Librairie, (original 1748, traduction 1796), vol. 1, p. 89-90

[3]

X. Gourdon and P. Sebah, Numbers, Constants and Computation, Internet site at http://numbers.computation.free.fr/Constants/constants.html.

[4]

E. Hairer and G. Wanner, L'analyse au fil de l'histoire, Bibliothèque Scopos, Springer, (2000)

[5]

G.H. Hardy and E.M. Wright, An Introduction to the Theory of Numbers, Oxford Science Publications, (1979)

[6]

A.S. Householder, The Numerical Treatment of a Single Nonlinear Equation, McGraw-Hill, New York, (1970)

[7]

V.J. Katz, A History of Mathematics-An Introduction, Addison-Wesley, (1998)

[8]

K. Knopp, Theory and application of infinite series, Blackie & Son, London, (1951)

[9]

I. Newton, Methodus fluxionum et serierum infinitarum, (1664-1671)

[10]

P. Ribenboim, The new Book of Prime Number Records, Springer, (1996)