|
Ц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 :
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
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 :
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 :
|
ж з
и
|
|
ц ч
ш
|
= |
ж з
и
|
|
ц ч
ш
|
|
ж з
и
|
|
ц ч
ш
|
|
|
so that
|
ж з
и
|
|
ц ч
ш
|
=An |
ж з
и
|
|
ц ч
ш
|
= |
ж з
и
|
|
ц ч
ш
|
n
|
|
ж з
и
|
|
ц ч
ш
|
. |
|
2.1.1 Quadratic convergence
This matrix representation may be used with n=2p and thanks to the
relation (exercise : show this by recurrence) :
comes the following algorithm
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
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 :
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 :
Easy manipulations show that
and that for the error
we have the bound
en+1 < en |
ж з
и
|
ж Ц
|
|
-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
which corresponds to A=50 and B=49 and the iterative system becomes :
as for the error, we have
en+1 < en |
ж з
и
|
ж Ц
|
|
-1 |
ц ч
ш
|
2
|
< en/9604. |
|
Here are the first iterates :
|
к к к к
к к к
|
|
|
|
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]).
It is more convenient to use the notation
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
(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
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 :
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 :
A nice improvement is to use the previous results on continuous fractions. We
can write :
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
|
|
ж Ц
|
|
, |
7
5
|
|
ж Ц
|
|
, |
17
12
|
|
ж Ц
|
|
, |
41
29
|
|
ж Ц
|
|
, |
99
70
|
|
ж Ц
|
|
, |
239
169
|
|
ж Ц
|
|
,ј |
ц ч
ш
|
. |
|
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 :
|
к к к к
к к к
|
|
|
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 :
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 :
|
к к к к
к к к
|
|
|
|
x4=1.41421356237(468...). |
|
|
|
|
Number D of correct digits after n iterations with the elementary Newton
iteration :
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 :
|
к к к к к
к к к к
|
|
|
|
|
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 :
|
к к к
к к
|
|
|
x3=1.414213562373095048(795...). |
|
|
|
|
Number D of correct digits after n iterations with Halley's iteration :
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 :
Number D of correct digits after n iterations with Householder's cubic
algorithm :
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
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
|
к к
к
|
|
x2=1.41421356237309(494...). |
|
|
|
|
Number D of correct digits after n iterations with the quartic algorithm :
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 :
|
к к
к
|
|
x2=1.414213562373095048801688(932...). |
|
|
|
|
Number D of correct digits after n iterations with the quintic algorithm :
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) )) , |
|
|
|
|
|
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 :
and
then the sequence hn tends to 0 and the sequence xn tends to
ЦA (A < 3).
Example 10
Let A=2, therefore
|
к к к к к
к к к к
|
|
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 :
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)
|