:: алгоритмы  и методы :: :: олимпиадные задачи :: :: связь :: :: о сайте ::
Путь: Математика » Быстрые вычисления » Квадратный корень из 2
  Квадратный корень из 2



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

Статья любезно предоставлена:
© Xavier Gordon, numbers.computation.free.fr



К статье прилагается исходник на Си

Ц2 = 1.41421356237309504880168872420969807856967187537694...
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 p2=2q2 therefore p is even and p=2pў which leads to q2=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 en=|Ц2-pn/qn|, one easily see that en+1 < en/5, which involves a geometrical convergence (that is, one digit will be added at every N iterations).

Starting with p0=q0=1 we obtain :
к
к
к
к
к
к
к
x1=1.(500ј),
x3=1.41(666...),
x6=1.4142(011...),
x9=1.414213(624...),

x3=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=2p 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 a0=b0=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 a0=b0=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 un=pn+qn and vn=qn then the double sequence (1) is equivalent to this new one :
м
н
о
un+1=2un+un-1
vn+1=un

and un/vn=un/un-1=(pn+qn)/qn converge to Ц2+1. So starting with u0=1,u1=2 gives the set of un:
(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 un 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/a1 with a1=2.4142135...=2+1/a2 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+pk/(5qk), the error is about 1/qk2. Observe that here pk=qk-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 (mk,nk) 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)=x2-2. It gives the following sequence :
x0=1,

xn+1=xn-  f(xn)

fў(xn)
=  1

2
ж
и
xn+  2

xn
ц
ш
=  xn

2
+  1

xn
,

xn+1 can also be interpreted as the simple mean between the approximation xn and 2/xn 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 x0 sufficiently close to the root of the equation.

Computing Ц2 with this algorithm is convenient if one can compute easily the inverse of xn. 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/x2-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)=x2-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 x0 (observe the similarity with Halley's iteration).

By applying this general pattern for f(x)=1/x2-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 x0=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/x2-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 x0=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 x0=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 hn tends to 0 and the sequence xn 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)