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



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



К статье есть исходник на Си

(c) Xavier Gourdon
web: numbers.computation.free.fr

log(2) = 0.69314718055994530941723212145817656807550013436025...

''(about the invention of logarithms) ... by shortening the labors, doubled the life of the astronomer''

Pierre-Simon de Laplace.


  1  Introduction



The story of logarithms really began with the Scot John Napier (1550-1617) in a work published in 1614. This treatise was in Latin and entitled ''Mirifici logarithmorum canonis descriptio'' [7]. This invention received an extraordinary welcome and George Gibson wrote during Napier Tercentenary Exhibition : ''... the invention of logarithms marks an epoch in the history of science''.

Very soon tables of logarithms were published. One of the first was due to the English mathematician Henry Briggs (1561-1631) who, in 1624, published his work in ''Arithmetica logarithmica'' [4]; the tables were to an accuracy of fourteen digits and containing the common logarithm of integers up to 20000 and from 90000 up to 101000. The remaining gap was completed [12] by the Dutchman Adrian Vlacq (1600-1667) in 1628.

The natural logarithm of a number x is now defined by :

Definition 1 Let x > 0
log(x)= у
х
x

1 
 dt

t
= у
х
x-1

0 
 dt

1+t

It may be interpreted as the area under the hyperbola y=1/t with t going from1 to x. (this geometric interpretation was showed by the Jesuit Grégoire de Saint-Vincent (1584-1667) in 1647).

Therefore, log(2) will be defined by


log(2)= у
х
2

1 
 dt

t
= у
х
1

0 
 dt

1+t
= у
х
1/2

0 
 dt

1-t
.

Theorem 2 log(2) is an irrational and transcendental number.


  2  Algorithms based on integral representation



From the integral representation of log(2) we may deduce some formulae. The first one is deduced form the rectangular approximation of the integral (here h=(b-a)/n)


у
х
b

a 
f(t) dt=h n
е
k=1 
f(a+kh)+O(h),
which, applied to f(t)=1/(1+t) on [0,1] gives


log(2)=
lim
n® Ґ 
ж
и
 1

n+1
+  1

n+2
+...+  1

2n
ц
ш
.

This sequence is increasing. Here are the numerical values of some partial sum for different value of n.


x10
=
0.6(687...)
x100
=
0.69(065...)
x1000
=
0.69(289...)

It's a very slow convergence, but it can be improved if we use the trapezoidal rule


у
х
b

a 
f(t)dt=  h

2
ж
и
f(a)+f(b)+2 n-1
е
k=1 
f(a+kh) ц
ш
+O(h2),
yielding
log(2)=
lim
n® Ґ 
ж
и
 3

4n
+ ж
и
 1

n+1
+  1

n+2
+ј+  1

2n-1
ц
ш
ц
ш

and some iterates are


x10
=
0.693(771...)
x100
=
0.6931(534...)
x1000
=
0.693147(243...)

The rate of convergence has been improved but it remains logarithmic. We now give a last improvement on the computation of the integral, known as the Simpson's rule


у
х
b

a 
f(t)dt=  h

6
n
е
k=0 
ж
и
f(a+kh)+f(a+(k-1)h)+4f ж
и
a+(k-  1

2
)h ц
ш
ц
ш
+O(h4),
giving
log(2)=
lim
n® Ґ 
n
е
k=0 
ж
и
 1

n+k
+  1

n+k-1
+  4

n+k-1/2
ц
ш

and the new iterates


x10
=
0.693147(374...)
x100
=
0.6931471805(794...)
x1000
=
0.69314718055994(726...)

From those examples we can notice that such formulae may only be used to compute a few digits of log(2), the rate of convergence is too slow to find constants with high precision.


  3  A simple limit



We know that the constant e may be defined by the famous limit


e=
lim
n® Ґ 
ж
и
1+  1

n
ц
ш
n

 
,

and it's interesting to note that a similar formula exists for log(2):


log(2)=
lim
n® Ґ 
n (21/n-1),

it's obviously deduced from the expansion


21/n=elog(2)/n=1+  log(2)

n
+O ж
и
 1

n2
ц
ш
.

With n=1000, we find


x1000=0.693(387...).

A better formula is given by


log(2)=
lim
n® Ґ 
 n

2
(21/n-2-1/n),

again, with n=1000


x1000=0.693147(236...),

the error is O(1/n2). To conclude this section, we give the unusual infinite product published by Seidel in 1871


 log(x)

x-1
= k=Ґ
Х
k=0 
 2

x1/2k+1
,

we apply it for x=2


log(2)=  2

21/2+1
 2

21/4+1
 2

21/8+1
...


  4  Taylor's expansion



In 1668, Nicolas Mercator (1620-1687) published in his Logarithmotechnia [6] the well known series expansion for the logarithmic function :

Theorem 3 (Mercator)
log(1+x)=x-  x2

2
+  x3

3
-...=
е
k і 1 
 (-1)k-1xk

k
.       -1 < x Ј 1

The Mercator series is very similar to another series studied by James Gregory (1638-1675) at the same time :


arctan(x)=x-  x3

3
+  x5

5
-ј

One of the first mathematician to use Mercator's series is Isaac Newton (1643-1727) who made various computation of logarithms in the celebrated ''De Methodis Serierum et Fluxionum'' written in 1671 but only published and translated from Latin to English in 1736, by John Colson. He made several accurate computations for small values of x, like ±0.1,±0.2 and was then able to compute logarithms of larger numbers using relations like log(2)=2log(1.2)-log(0.9)-log(0.8).

Applying the Mercator series for x=1 leads to the famous and very slowly convergent sequence
log(2)=1-  1

2
+  1

3
-  1

4
+  1

5
-...

which first values, taking respectively 10, 100, 1000 and 10000 terms are :


x10
=
0.6(456...)
x100
=
0.6(881...)
x1000
=
0.69(264...)
x10000
=
0.693(097...).

A much better idea is to apply the series with x=-1/2, this produces a formula that was often used for numerical computation:


log(2)=
е
k і 1 
 1

k2k
(1)

The rate of convergence is now geometric (p iterations are needed to obtain one more digits).


x1
=
0.(5)
x10
=
0.693(064...)
x100
=
0.6931471805599453094172321214581(688...)

Applying the formula up to k=1000 gives more than 300 digits for log(2). The number of iterations k needed to obtain d decimal digits is given by (with N=2 in the previous formula)


 1

10d
=  1

kNk

which gives for large values of k


k »  d

log10(N)

If we use the trivial decomposition 2=(4/3)(3/2) and take the logarithm of both side of the equality we find a new formula


log(2)=
е
k і 1 
 1

k
ж
и
 1

3k
+  1

4k
ц
ш

The partial sum of this series up to k=1000 will give more than 470 digits.

From the relation log(2)=-1/2log(1-1/4)+ arctanh(1/2), (see next paragraph for the definition of  arctanh) we deduce the relations


log(2)
=

е
k і 0 
ж
и
 1

8k+8
+  1

4k+2
ц
ш
 1

4k
,
log(2)
=
 2

3
+
е
k і 1 
ж
и
 1

2k
+  1

4k+1
+  1

8k+4
+  1

16k+12
ц
ш
 1

16k
.

Those formulae can also be used to compute directly the n-th binary digit of log(2) without computing the previous ones.


  5  Machin like formulae



Using Mercator series has improved a lot our ability to compute many digits for log(2), but like for the number p it's possible to go further. We recall the definition of the inverse hyperbolic tangent.

Definition 4 Let |x| < 1:
 arctanh(x)=  1

2
log ж
и
 1+x

1-x
ц
ш
=
е
k і 0 
 x2k+1

2k+1
.

For x=1/3, this leads easily to the basic formula log(2)=2 arctanh(1/3) which is similar to the formula p = arctan(1). We rewrite it as :


log(2)=  2

3

е
k і 0 
 1

(2k+1)9k

Two iterates are then :


x1
=
0.6(666...)
x10
=
0.6931471805(498...).

Applying the formula up to k=1000 gives more than 950 digits for log(2). Is it possible to do better? There is a famous formula for p,this formula is known as Machin's formula and is given by the celebrated relation :


 p

4
=4arctan ж
и
 1

5
ц
ш
-arctan ж
и
 1

239
ц
ш

The idea is to look for relations like this one for log(2). If we note (a1,a2,...,an) a set of rational numbers and (p1,p2,...,pn) a set of increasing integers we are looking for relations like


log(2)= n
е
i=1 
ai arctanh ж
и
 1

pi
ц
ш

We also want the formula to be efficient. A good relation is a compromise between a small number of terms (n must be small) and large values for the pi. According to Lehmer, we define the efficiency E by

Definition 5 (Lehmer's measure). Let pi > 0 a set of n numbers, then :
E= n
е
i=1 
 1

log10(pi2)
.

For example the efficiency of log(2)=2 arctanh(1/3) (n=1,p1=3) is 1.05. With the same definition the efficiency of Machin's formula is 0.926 (n=2,p1=5,p2=239), so it's a little bit more efficient than the formula for log(2).

It's possible to show the following decomposition theorem :

Theorem 6 Let p > 1:
arctanh ж
и
 1

p
ц
ш
=
arctanh ж
и
 1

2p-1
ц
ш
arctanh ж
и
 1

2p+1
ц
ш
arctanh ж
и
 1

p
ц
ш
=
arctanh ж
и
 1

2p-1
ц
ш
- arctanh ж
и
 1

2p2-1
ц
ш
arctanh ж
и
 1

p
ц
ш
=
arctanh ж
и
 1

2p+1
ц
ш
arctanh ж
и
 1

2p2-1
ц
ш
arctanh ж
и
 1

p
ц
ш
=
arctanh ж
и
 1

2p
ц
ш
arctanh ж
и
 1

4p3-3p
ц
ш

Applying this theorem for p=3 gives relations with n=2 :

Corollary 7
log(2)
=
arctanh ж
и
 1

5
ц
ш
+2 arctanh ж
и
 1

7
ц
ш
       Euler 1748
(2)
log(2)
=
arctanh ж
и
 1

5
ц
ш
-2 arctanh ж
и
 1

17
ц
ш
(3)
log(2)
=
arctanh ж
и
 1

7
ц
ш
+2 arctanh ж
и
 1

17
ц
ш
(4)
log(2)
=
arctanh ж
и
 1

6
ц
ш
+2 arctanh ж
и
 1

99
ц
ш
(5)

The efficiencies are respectively E=1.307,E=1.121,E=0.998 and E=0.893. Using other decomposition formulas and computer calculation it's possible to find many other relations of this nature (given in [8]). To illustrate this, we select another formula with 3 terms and two others with 4 terms.

Theorem 8 (Sebah 1997). The constant log(2) may be obtained by one of the following fast converging series :
18 arctanh ж
и
 1

26
ц
ш
-arctanh ж
и
 1

4801
ц
ш
+8 arctanh ж
и
 1

8749
ц
ш
,
(6)
144 arctanh ж
и
 1

251
ц
ш
+54 arctanh ж
и
 1

449
ц
ш
-38 arctanh ж
и
 1

4801
ц
ш
+62 arctanh ж
и
 1

8749
ц
ш
,
(7)
72 arctanh ж
и
 1

127
ц
ш
+54 arctanh ж
и
 1

449
ц
ш
+34 arctanh ж
и
 1

4801
ц
ш
-10 arctanh ж
и
 1

8749
ц
ш
.
(8)

The efficiency are respectively (E=0.616,E=0.659,E=0.689). The two last formulae were used to compute more than 108 digits for log(2). One was used for the computation and the other for the verification and , as you can see, only the computation of 5 arctanh were necessary because of the similitude of the 2 formulae. Thanks to relation (6), the record was increased up to more than 5.108 a few years later.

Note that it's also possible to perform separately the computation of each term, making those relations parallelisable.

We give, for example, the first iterates of the second formula (the one with 251):


x1
=
0.69314(394...)
x2
=
0.6931471805(304...)
x3
=
0.69314718055994(497...)
x4
=
0.69314718055994530941(317...)

and each iteration will add about 4.8 digits.

5.1  Formula with rational numbers

It may be convenient to look for relations with rational numbers as arguments of the arctanh function, that is identities of the form :
log(2)= n
е
i=1 
ai arctanh ж
и
 mi

pi
ц
ш

with mi being integers eventually different of 1. Lehmer's measure must be adapted to take this in account
E= n
е
i=1 
 1

2log10(pi/mi)
,

but it only gives a rough estimation.

Theorem 9 (Sebah 1997).
log(2)
=
6 arctanh ж
и
 1

9
ц
ш
+2 arctanh ж
и
 3

253
ц
ш
,
(9)
log(2)
=
10 arctanh ж
и
 1

17
ц
ш
+4 arctanh ж
и
 13

449
ц
ш
,
(10)

with respective efficiency E=0.784 and E=0.722.

Note that relation (10) was used for several high precision computation of the constant.


  6  Variation with hypergeometric sequences



Gauss studied the following family of series


F(a,b;c;z)=1+  a.b

c
 z

1!
+  a(a+1).b(b+1)

c(c+1)
 z2

2!
+...

where a,b,c are real numbers. The series converges in |z| < 1. On the circle |z|=1, the series converges when c > a+b. Those functions are called hypergeometric functions [1]. Many usual functions can be represented as hypergeometric functions with suitable values for a,b,c.

Example 10
ln(1+x)
=
xF(1,1;2;-x)
 arctanh(x)
=
xF(  1

2
,1;  3

2
;x2)
arctan(x)
=
xF(  1

2
,1;  3

2
;-x2)
 1

1-x
=
F(1,1;1;x)

In 1797, Johann Friederich Pfaff (1765-1825), a friend of Gauss, gave the interesting relation


F(a,b;c;z)=(1-z)-bF(c-a,b;c;  z

z-1
).

If we apply this identity to the arctanh function


arctanh(x)=xF(  1

2
,1;  3

2
;x2),

we find


arctanh(x)=  x

(1-x2)
F(1,1;  3

2
;  x2

x2-1
),

giving a new sequence for this function


arctanh(x)=  x

(1-x2)
ж
и
1-  2

3
ж
и
 x2

1-x2
ц
ш
+  2.4

3.5
ж
и
 x2

1-x2
ц
ш
2

 
-  2.4.6

3.5.7
ж
и
 x2

1-x2
ц
ш
3

 
+... ц
ш

and with log(2)=2 arctanh(1/3); after some easy manipulations comes the series :

Corollary 11
log(2)=  3

4
ж
и
1-  1

12
+  1.2

12.20
-  1.2.3

12.20.28
+... ц
ш
.

Each new term of this series is multiplied by -k/(8k+4) where k=1,2,... The same kind of formula was given by Euler to compute p. One nice application of this corollary is to make it feasible to write a tiny code to compute log(2) just like it was done for p.


  7  log(2) and the AGM



The rate of convergence of all the previous series were logarithmic or geometric. It's possible to find for log(2) just as for p, a sequence that will have quadratic convergence, based on the AGM (Arithmetic-Geometric Mean) (see [2]).

The most classical among the AGM based quadratic convergence algorithms for log(2) is the following. Starting from a0 and b0 > 0, we consider the AGM iteration
an+1=  an+bn

2
,       bn+1=
Ц
 

anbn
 
,
and we use the notation
R(a0,b0)= ж
з
и
 1

1-
е
n і 0 
2n-1(an2-bn2)
ц
ч
ш
.

Theorem 12 Let N і 3 and 1/2 Ј x Ј 1,
| log(x)-R(1,10-N)+R(1,10-Nx)| Ј  N

102(N-2)
.
(11)

By choosing x=1/2, this algorithm gives about 2N decimal digits of log(2). The convergence is quadratic and should be stopped when 2n-1(an2-bn2) is less than 10-2N, which occurs for n » log2(N).

Notice that this algorithm is not specific to the computation of log(2) and can be used to evaluate log(x) for any real number x. Using FFT multiplication, its complexity is O(nlog2(n)) to obtain n digits of x.


  8  Computation of log(p)



Starting with a good approximation of log(2), it's easy to find the logarithm of all the integers thanks to the relation


2 arctanh ж
и
 1

2p+1
ц
ш
=log ж
з
з
з
з
и
1+  1

2p+1

1-  1

2p+1
ц
ч
ч
ч
ч
ш
=log ж
и
 p+1

p
ц
ш
,

thus we deduce that :

Theorem 13 Let p > 0
log(p+1)=log(p)+2 ж
и
 1

2p+1
+  1

3
 1

(2p+1)3
+  1

5
 1

(2p+1)5
+... ц
ш
.

This may be useful when you need to compute the logarithms of many successive integers.

Example 14 With p=2, p=4 and p=7
log(3)
=
log(2)+  2

5
ж
и
1+  1

3
 1

25
+  1

5
 1

252
+... ц
ш
,
log(5)
=
2log(2)+  2

9
ж
и
1+  1

3
 1

81
+  1

5
 1

812
+... ц
ш
,
log(7)
=
3log(2)-  2

15
ж
и
1+  1

3
 1

225
+  1

5
 1

2252
+... ц
ш
.


  9  Records of computation for log(2)



Number of digits When Who Notes
16 ~ 1671 I. Newton He used log(2)=2log(1.2)-log(0.9)-log(0.8) and Mercator's series for log.
25 1748 L. Euler Relation (2) was used. Euler also computed the natural logarithm of integers from to 3 to 10 with 25 digits [5].
137 1853 W. Shanks Shanks also computed log(3), log(5) and log(10) with 137 digits. His work was published with a p calculation  [9].
273 1878 Adams Adams also computed log(3), log(5) and log(7) with the same precision.
330 1940 H. S. Uhler Uhler also computed log(3), log(5), log(7) and log(17) with the same precision [11].
3 683 1962 D.W. Sweeney The computation of log(2) was necessary for a computation of Euler's constant up to 3566 digits in the same article). The formula used for log2 was the expansion (1) (see [10]).
2 015 926 1997 P. Demichel The computation used a classical approach on formula (2). It took 16 hours on a HP9000/871 at 160MHz.
5 039 926 1997 P. Demichel The same technique was used and the computation took 130 hours on the same computer.
10 079 926 1997 P. Demichel Same technique, 400 hours on the same computer.
29 243 200 1997 Dec X. Gourdon The AGM formulae (11) was used (with FFT techniques). The computation was done on a SGI R10000 workstation and took 20 hours 58 minutes.
58 484 499 1997 Dec X. Gourdon Same AGM technique, 83 hours on the same computer.
108 000 000 1998 Sep X. Gourdon The four term formula (7) was used together with a binary splitting process. Formula (8) was used for the verification. The computation took 47 hours on a PII 450 with 256 Mo.
200 001 000 2001 Sep X. Gourdon & S. Kondo Formulae (6) and (5) were used. The computation took 8 hours on an Athlon 1300 with 2 Go.
240 000 000 2001 Sep X. Gourdon & P. Sebah Formulae (6) and (10) were used. The computation took 14 hours on a PIV 1400 with 1024 Mo.
500 000 999 2001 Sep X. Gourdon & S. Kondo Formulae (6) and (10) were used. The computation took 28 hours on a PIV 1400 with 1024 Mo.


  References



[1]
M. Abramowitz and I. Stegun, Handbook of Mathematical Functions, Dover, New York, (1964)

[2]
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)

[3]
R.P. Brent, Fast multiple-Precision evaluation of elementary functions, J. Assoc. Comput. Mach., (1976), vol. 23, p.242-251

[4]
H.Briggs, Arithmetica logarithmica sive logarithmorum Chiliades Triginta, London, (1624)

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

[6]
N. Mercator, Logarithmotechnia, London, (1668)

[7]
J. Napier, Mirifici logarithmorum canonis descriptio, Edimboug, (1614)

[8]
P. Sebah, Machin like formulae for logarithm, Unpublished, (1997).

[9]
W. Shanks, Contributions to Mathematics Comprising Chiefly the Rectification of the Circle to 607 Places of Decimals, G. Bell, London, (1853)

[10]
D.W. Sweeney, On the Computation of Euler's Constant, Mathematics of Computation, (1963), p. 170-178

[11]
H.S. Uhler, Recalculation and extension of the modulus and of the logarithms of 2, 3, 5, 7 and 17, Proc. Nat. Acad. Sci., (1940), vol. 26, p. 205-212

[12]
A. Vlacq, Arithmetica logarithmica, Gouda,28)