ex-2: complete review

- general review and proofread
- complet the explanation of Lambert W based solutions
- add section on 1 million digits computation
- add references to GMP and mpmath
This commit is contained in:
Michele Guerini Rocco 2020-05-27 15:29:59 +02:00
parent fcbde10f9b
commit f2be522089
3 changed files with 267 additions and 218 deletions

View File

@ -198,3 +198,22 @@
computation of Euler's constant}, computation of Euler's constant},
year={2017} year={2017}
} }
@manual{mpmath13,
key={mpmath},
author={Fredrik Johansson and others},
title={mpmath: a {P}ython library for arbitrary-precision floating-point arithmetic (version 0.18)},
note={{\tt http://mpmath.org/}},
month={December},
year={2013},
}
@manual{gmp20,
key={gmp},
author={Torbjörn Granlund and the GMP development team},
title={GNU MP: The GNU Multiple Precision Arithmetic Library,
Edition 6.2.0 (version 0.18)},
note={{\tt https://gmplib.org/}},
month={January},
year={2020},
}

View File

@ -43,11 +43,11 @@ $$
| \gamma(n_{i+1}) - \gamma | > | \gamma(n_i) - \gamma| | \gamma(n_{i+1}) - \gamma | > | \gamma(n_i) - \gamma|
$$ $$
and $\gamma (n_i)$ was selected as the best result (see @tbl:1_results). and $\gamma (n_i)$ was selected as the best result (see @tbl:naive-errs).
--------------------------------------------- ----------------------------
n $|\gamma(n)-\gamma|$ n $|γ(n)-γ|$
---------------------- ---------------------- ----------- ----------------
\num{2e1} \num{2.48e-02} \num{2e1} \num{2.48e-02}
\num{2e2} \num{2.50e-03} \num{2e2} \num{2.50e-03}
@ -69,14 +69,14 @@ n $|\gamma(n)-\gamma|$
\num{2e10} \num{2.42e-11} \num{2e10} \num{2.42e-11}
\num{2e11} \num{1.44e-08} \num{2e11} \num{1.44e-08}
--------------------------------------------- ----------------------------
Table: Partial results using the definition of $\gamma$ with double Table: Partial results using the definition of $\gamma$ with double
precision. {#tbl:1_results} precision. {#tbl:naive-errs}
The convergence is logarithmic: to fix the first $d$ decimal places, about The convergence is logarithmic: to fix the first $d$ decimal places, about
$10^d$ terms of the armonic series are needed. The double precision runs out at $10^d$ terms of the harmonic series are needed. The double precision runs out
the $10^{\text{th}}$ place, at $n=\num{2e10}$. at the $10^{\text{th}}$ place, at $n=\num{2e10}$.
Since all the number are given with double precision, there can be at best 16 Since all the number are given with double precision, there can be at best 16
correct digits, since for a double 64 bits are allocated in memory: 1 for the correct digits, since for a double 64 bits are allocated in memory: 1 for the
sign, 8 for the exponent and 55 for the mantissa: sign, 8 for the exponent and 55 for the mantissa:
@ -87,43 +87,40 @@ $$
Only 10 digits were correctly computed: this means that when the terms of the Only 10 digits were correctly computed: this means that when the terms of the
series start being smaller than the smallest representable double, the sum of series start being smaller than the smallest representable double, the sum of
all the remaining terms gives a number $\propto 10^{-11}$. The best result is all the remaining terms gives a number $\propto 10^{-11}$. The best result is
shown in @tbl:first. shown in @tbl:naive-res.
--------- ----------------------- ------- --------------------
true: 0.57721\ 56649\ 01533 exact 0.57721 56649 01533
approx 0.57721 56648 77325
approx: 0.57721\ 56648\ 77325 diff 0.00000 00000 24207
------- --------------------
diff: 0.00000\ 00000\ 24207
--------- -----------------------
Table: First method best result. From the top down: true value, best estimation Table: First method best result. From the top down: true value, best estimation
and difference between them. {#tbl:first} and difference between them. {#tbl:naive-res}
### Alternative formula ### Alternative formula
As a first alternative, the constant was computed through the identity which As a first alternative, the constant was computed through the identity
relates $\gamma$ to the $\Gamma$ function as follow [@davis59]: [@davis59] which relates $\gamma$ to the $\Gamma$ function as follow:
$$ $$
\gamma = \lim_{M \rightarrow + \infty} \sum_{k = 1}^{M} \gamma = \lim_{M \rightarrow + \infty} \sum_{k = 1}^{M}
\binom{M}{k} \frac{(-1)^k}{k} \ln(\Gamma(k + 1)) \binom{M}{k} \frac{(-1)^k}{k} \ln(\Gamma(k + 1))
$$ $$
Varying $M$ from 1 to 100, the best result was obtained for $M = 41$ (see Varying $M$ from 1 to 100, the best result was obtained for $M = 41$ (see
@tbl:second). It went sour: the convergence is worse than using the definition @tbl:limit-res). This approximation gave an underwhelming result: the
itself. Only two places were correctly computed (@tbl:second). convergence is actually worse than the definition itself. Only two decimal
places were correctly computed:
--------- ----------------------- -------- --------------------
true: 0.57721\ 56649\ 01533 exact 0.57721 56649 01533
approx 0.57225 72410 34058
approx: 0.57225\ 72410\ 34058 diff 0.00495 84238 67473
-------- --------------------
diff: 0.00495\ 84238\ 67473
--------- -----------------------
Table: Best estimation of $\gamma$ using Table: Best estimation of $\gamma$ using
the alternative formula. {#tbl:second} the alternative formula. {#tbl:limit-res}
Here, the problem lies in the binomial term: computing the factorial of a Here, the problem lies in the binomial term: computing the factorial of a
number greater than 18 goes over 15 places and so cannot be correctly number greater than 18 goes over 15 places and so cannot be correctly
@ -146,15 +143,15 @@ $$
\prod_{k = 1}^{+ \infty} \left( 1 + \frac{z}{k} \right) e^{-z/k} \right) \prod_{k = 1}^{+ \infty} \left( 1 + \frac{z}{k} \right) e^{-z/k} \right)
$$ $$
The execution stops when there is no difference between two consecutive therms The execution stops when there is no difference between two consecutive terms
of the infinite product (it happens for $k = 456565794 \sim \num{4.6e8}$, of the infinite product (it happens for $k = 456565794 \sim \num{4.6e8}$,
meaning that for this value of $k$ the term of the product is equal to 1 in meaning that for this value of $k$ the term of the product is equal to 1 in
terms of floating points). Different values of $z$ were checked, with $z_{i+1} terms of floating points). Different values of $z$ were checked, with $z_{i+1}
= z_i + 0.01$ ranging from 0 to 20, and the best result was found for $z = 9$. = z_i + 0.01$ ranging from 0 to 20, and the best result was found for $z = 9$.
--------------------------------------------------------------- -----------------------------------------------
z $|\gamma(z) - \gamma |$ z $|\gamma(z) - \gamma |$ z $|γ(z) - γ|$ z $|γ(z) - γ|$
----- ------------------------ ------ ------------------------ ----- ---------------- ------ ----------------
1 \num{9.712e-9} 8.95 \num{9.770e-9} 1 \num{9.712e-9} 8.95 \num{9.770e-9}
3 \num{9.320e-9} 8.96 \num{9.833e-9} 3 \num{9.320e-9} 8.96 \num{9.833e-9}
@ -174,41 +171,38 @@ terms of floating points). Different values of $z$ were checked, with $z_{i+1}
17 \num{9.971e-9} 9.03 \num{9.264e-9} 17 \num{9.971e-9} 9.03 \num{9.264e-9}
19 \num{10.084e-9} 9.04 \num{9.419e-9} 19 \num{10.084e-9} 9.04 \num{9.419e-9}
--------------------------------------------------------------- -----------------------------------------------
Table: Differences between some obtained values of $\gamma$ and Table: Differences between some obtained values of $\gamma$ and
the exact one found with the reciprocal $\Gamma$ function formula. the exact one found with the reciprocal $\Gamma$ function formula.
The values on the left are shown to give an idea of the $z$ The values on the left are shown to give an idea of the $z$
large-scale behaviour; on the right, the values around the best large-scale behaviour; on the right, the values around the best
one ($z = 9.00$) are listed. {#tbl:3_results} one ($z = 9.00$) are listed. {#tbl:recip-errs}
As can be seen in @tbl:3_results, the best value for $z$ is only by chance, As can be seen in @tbl:recip-errs, the best value for $z$ is only by chance,
since all $|\gamma(z) - \gamma |$ are of the same order of magnitude. The best since all $|\gamma(z) - \gamma |$ are of the same order of magnitude. The best
one is compared with the exact value of $\gamma$ in @tbl:third. one is compared with the exact value of $\gamma$ in @tbl:recip-res.
--------- ----------------------- ------- --------------------
true: 0.57721\ 56649\ 01533 true 0.57721 56649 01533
approx 0.57721 56564 18607
diff 0.00000 00084 82925
------- --------------------
approx: 0.57721\ 56564\ 18607 Table: Third method results for z = 9.00. {#tbl:recip-res}
diff: 0.00000\ 00084\ 82925 In this approximation, the convergence of the infinite product is fast enough
--------- ----------------------- to reach the $8^{th}$ decimal place.
Table: Third method results for z = 9.00. {#tbl:third}
This time, the convergence of the infinite product is fast enough to ensure the
$8^{th}$ place.
### Fastest convergence formula ### Fastest convergence formula
The fastest known convergence belongs to the following formula, known as refined The fastest known convergence to $\gamma$ belongs to the following formula,
Brent-McMillan formula[@yee19]: known as refined Brent-McMillan formula [@yee19]:
$$ $$
\gamma_N = \frac{A(N)}{B(N)} -\frac{C(N)}{B^2(N)} - \ln(N) \gamma_N = \frac{A(N)}{B(N)} -\frac{C(N)}{B^2(N)} - \ln(N)
$$ {#eq:faster} $$ {#eq:faster}
where
with:
\begin{align*} \begin{align*}
&A(N) = \sum_{k=1}^{+ \infty} \frac{N^k}{k!} \cdot H(k) &A(N) = \sum_{k=1}^{+ \infty} \frac{N^k}{k!} \cdot H(k)
\with H(k) = \sum_{j=1}^{k} \frac{1}{j} \\ \with H(k) = \sum_{j=1}^{k} \frac{1}{j} \\
@ -217,95 +211,119 @@ with:
\frac{((2k)!)^3}{(k!)^4 \cdot (16k)^2k} \\ \frac{((2k)!)^3}{(k!)^4 \cdot (16k)^2k} \\
\end{align*} \end{align*}
and the difference between $\gamma_N$ and $\gamma$ is given by [@demailly17]: The asymptotic error of this estimation, given in [@demailly17], is:
$$ $$
|\gamma_N - \gamma| < \frac{5 \sqrt{2 \pi}}{12 \sqrt{N}} e^{-8N} |\gamma_N - \gamma| \sim \frac{5 \sqrt{2 \pi}}{12 \sqrt{N}} e^{-8N} = \Delta_N
$$ {#eq:NeD} $$ {#eq:NeD}
This gives the value of $N$ which is to be used into the formula in order to get The error bound gives the value of $N$ to be used to get $D$ correct decimal
$D$ correct decimal digits of $\gamma$. In fact, by imposing: digits of $\gamma$. In fact, this is done by imposing:
$$ $$
\frac{5 \sqrt{2 \pi}}{12 \sqrt{N}} e^{-8N} < 10^{-D} \frac{5 \sqrt{2 \pi}}{12 \sqrt{N}} e^{-8N} < 10^{-D}
$$ $$
With a bit of math, it can be found that the smallest integer which satisfies The inequality with the equal sign can be solved with the use of the Lambert
the inequality is given by: $W$ function. For a real number $x$, $W(x)$ is defined as the inverse function
of $x\exp(x)$, so the idea is to recast the equation into this form and take
$W$ both sides.
\begin{align*}{3}
\frac{5 \sqrt{2 \pi}}{12 \sqrt{x}} e^{-8x} = 10^{-D}
& \thus \left(\frac{12}{5}\right)^2 \frac{x}{2 \pi} e^{16x} = 10^{2D} \\
& \thus 16x e^{16x} = \left(\frac{5 \pi}{9}\right)^2 10^{2D + 1} \\
& \thus x = \frac{1}{16} W\left(\left(\frac{5 \pi}{9}\right)^2 10^{2D + 1}\right)
\end{align*}
The smallest integer which satisfies the inequality is then
$$ $$
N = 1 + \left\lfloor \frac{1}{16} N = 1 + \left\lfloor \frac{1}{16}
W \left( \frac{5 \pi}{9} 10^{2D + 1}\right) \right\rfloor W \left( \frac{5 \pi}{9} 10^{2D + 1}\right) \right\rfloor
$$ $$ {#eq:refined-n}
The series $A$ and $B$ were computed till there is no difference between two The series $A$ and $B$ were computed till there is no difference between two
consecutive terms. Results are shown in @tbl:fourth. consecutive terms. Results are shown in @tbl:fourth.
--------- ------------------------------ ------- --------------------------
true: 0.57721\ 56649\ 01532\ 75452 exact 0.57721 56649 01532 75452
approx 0.57721 56649 01532 53248
diff 0.00000 00000 00000 33062
------- --------------------------
approx: 0.57721\ 56649\ 01532\ 53248 Table: $\gamma$ estimation with the fastest known convergence formula
(@eq:faster). {#tbl:fourth}
diff: 0.00000\ 00000\ 00000\ 33062 Due to roundoff errors, the best result was obtained only for $D = 15$, for
--------- ------------------------------
Table: $\gamma$ estimation with the fastest
known convergence formula (@eq:faster). {#tbl:fourth}
Because of roundoff errors, the best result was obtained only for $D = 15$, for
which the code accurately computes 15 digits and gives an error of which the code accurately computes 15 digits and gives an error of
\num{3.3e16}. For $D > 15$, the requested can't be fulfilled. \num{3.3e16}. For $D > 15$, the requested can't be fulfilled.
### Arbitrary precision ### Arbitrary precision
To overcome the issues related to the double representation, one can resort to a To overcome the issues related to the limited precision of the machine floating
representation with arbitrary precision. In the GMP library (which stands for points, one can resort to a software implementation with arbitrary precision.
GNU Multiple Precision arithmetic), for example, real numbers can be
approximated by a ratio of two integers (a fraction) with arbitrary precision: In the GMP [@gmp20] library (GNU Multiple Precision arithmetic), for
this means a check is performed on every operation and in case of an integer example, reals can be approximated by a rational or floating point numbers
overflow, additional memory is requested and used to represent the larger result. with arbitrary precision. For integer and fractions this means a check is
Additionally, the library automatically switches to the optimal algorithm performed on every operation that can overflow and additional memory is
to compute an operation based on the size of the operands. requested and used to represent the larger result. For floating points it
means the size of the mantissa can be chosen before initialising the number.
Additionally, the library automatically switches to the optimal algorithm to
compute an operation based on the size of the operands.
The terms in @eq:faster can therefore be computed with arbitrarily large The terms in @eq:faster can therefore be computed with arbitrarily large
precision. Thus, a program that computes the Euler-Mascheroni constant within precision. Thus, a program that computes the Euler-Mascheroni constant within a
a user controllable precision was implemented. user controllable precision was implemented.
To compute $N$ by @eq:refined-n, a trick must be used to avoid overflows in the
exponentiation, particularly when computing more than a few hundreds digits are
to be computed. The details are explained in @sec:optimised, below. Once the
number $N$ has been fixed, the series $A(N)$ and $B(N)$ are to be evaluated.
With rational numbers, a different criterion for the truncation must be
considered, because two consecutive term are always different. Brent and
McMillan[@brent00] prove that to reduce the partial sum to less than $\Delta_N$
it is sufficient to compute $\alpha N$ terms, where $\alpha$ is the solution to
the equation:
\begin{align*}
\alpha\ln(\alpha) = 3 + \alpha &&
\alpha = \exp(W(3e) - 1) \approx 4.97
\end{align*}
According to [@brent00], the $A(N)$, $B(N)$ and $C(N)$ series were computed up
to $k_{\text{max}} = 5N$ \textcolor{red}{sure?}, since it guarantees the accuracy
od the result up to the $D$ decimal digit.
The GMP library offers functions to perform some operations such as addition, The GMP library offers functions to perform some operations such as addition,
multiplication, division, etc. However, the logarithm function is not multiplication, division, etc. However, the logarithm function is not
implemented. Thus, most of the code carries out the $\ln(N)$ computation. implemented. Thus, most of the code carries out the $\ln(N)$ computation.
This is because the logarithm of only some special numbers can be computed The reason is that the logarithm of only some special numbers can be computed
with arbitrary precision, namely the ones of which a converging series is known. with arbitrary precision, namely the ones of which a converging series is known.
This forces $N$ to be rewritten in the following way: This forces $N$ to be rewritten in the following way:
$$ $$
N = N_0 \cdot b^e \thus \ln(N) = \ln(N_0) + e \cdot \ln(b) N = N_0 \cdot b^e \thus \ln(N) = \ln(N_0) + e \cdot \ln(b)
$$ $$
Since a fast converging series for $\ln(2)$ is known (it will be shwn shortly), Since a fast converging series for $\ln(2)$ is known (it will be shown shortly),
$b = 2$ was chosen. In case $N$ is a power of two, only $\ln(2)$ is to be $b = 2$ was chosen. If $N$ is a power of two, $N_0$ is 1 and only $\ln(2)$ is to be
computed. If not, it must be handled as follows. computed. More generally, the problem reduces to the calculation of $\ln(N_0)$.
As well as for the scientific notation, in order to get the mantissa $1
Regarding the scientific notation, to find the mantissa $1
\leqslant N_0 < 2$, the number of binary digits of $N$ must be computed \leqslant N_0 < 2$, the number of binary digits of $N$ must be computed
(conveniently, a dedicated function `mpz_sizeinbase()` is found in GMP). If the (conveniently, a dedicated function `mpz_sizeinbase()` exists in GMP). If the
digits are digits are
$n$: $n$:
$$ $$
e = n - 1 \thus N = N_0 \cdot 2^{n-1} \thus N_0 = \frac{N}{2^{n - 1}} e = n - 1 \thus N = N_0 \cdot 2^{n-1} \thus N_0 = \frac{N}{2^{n - 1}}
$$ $$
The logarithm of whichever number $N_0$ can be computed using the series of The logarithm of $N_0$ can be computed from the Taylor series of the hyperbolic
$\text{atanh}(x)$, which converges for $|x| < 1$: tangent, which is convergent for $|x| < 1$:
$$ $$
\text{atanh}(x) = \sum_{k = 0}^{+ \infty} \frac{x^{2k + 1}}{2x + 1} \text{atanh}(x) = \sum_{k = 0}^{+ \infty} \frac{x^{2k + 1}}{2x + 1}
$$ $$
In fact: The relation with the logarithm follows from the definition
$$ $$
\text{tanh}(x) = \frac{e^x - e^{-x}}{e^x + e^{-x}} \text{tanh}(x) = \frac{e^x - e^{-x}}{e^x + e^{-x}}
= \frac{e^{2x} - 1}{e^{2x} + 1} = \frac{e^{2x} - 1}{e^{2x} + 1}
$$ $$
and the change of variable $z = e^{2x}$:
with the change of variable $z = e^{2x}$, one gets:
$$ $$
\text{tanh} \left( \frac{\ln(z)}{2} \right) = \frac{z - 1}{z + 1} \text{tanh} \left( \frac{\ln(z)}{2} \right) = \frac{z - 1}{z + 1}
\thus \ln(z) = 2 \, \text{atanh} \left( \frac{z - 1}{z + 1} \right) \thus \ln(z) = 2 \, \text{atanh} \left( \frac{z - 1}{z + 1} \right)
@ -322,19 +340,19 @@ $$
= 2 \sum_{k = 0}^{+ \infty} \frac{y^{2k + 1}}{2k + 1} = 2 \sum_{k = 0}^{+ \infty} \frac{y^{2k + 1}}{2k + 1}
$$ $$
But when to stop computing the series? To estimate a series with a given precision some care must be taken:
Given a partial sum $S_k$ of the series, it is possible to know when a digit is different techniques apply to different series and are explained in
definitely correct. The key lies in the following concept [@riddle08]. Letting the paper [@riddle08]. In this case, letting $S$ be the value of the
$S$ be the value of the series: series and $S_k$ the $k$-th partial sum the following bounds can be found:
$$ $$
S_k + a_k \frac{L}{1 -L} < S < S_k + \frac{a_{k+1}}{1 - \frac{a_{k+1}}{a_k}} S_k + a_k \frac{L}{1 -L} < S < S_k + \frac{a_{k+1}}{1 - \frac{a_{k+1}}{a_k}}
$$ $$
where $a_k$ is the $k^{\text{th}}$ term of the series and $L$ is the limiting where $a_k$ is the $k^{\text{th}}$ term of the series and $L$ is the limiting
ratio of the series terms, which must be $< 1$ in order for it to converge (in ratio of the series terms, which must be $\le 1$ in order for it to converge (in
this case, it is easy to prove that $L = y^2 < 1$). The width $\Delta S$ of the this case, it is easy to prove that $L = y^2 < 1$). The width $\Delta S$ of the
interval containing $S$ gives the precision of the estimate $\tilde{S}$ if this interval containing $S$ gives the precision of the estimate $\tilde{S}$ if this
last is assumed to be the middle value of it, namely: last is assumed to be its middle value, namely:
$$ $$
\tilde{S} = S_k + \frac{1}{2} \left( \tilde{S} = S_k + \frac{1}{2} \left(
a_k \frac{L}{1 -L} + \frac{a_{k+1}}{1 - \frac{a_{k+1}}{a_k}} a_k \frac{L}{1 -L} + \frac{a_{k+1}}{1 - \frac{a_{k+1}}{a_k}}
@ -342,169 +360,185 @@ $$
\Delta S = \frac{a_{k+1}}{1 - \frac{a_{k+1}}{a_k}} - a_k \frac{L}{1 -L} \Delta S = \frac{a_{k+1}}{1 - \frac{a_{k+1}}{a_k}} - a_k \frac{L}{1 -L}
$$ $$
In order to know when to stop computing the series, $\Delta S$ must be For this series:
evaluated. With a bit of math: \begin{align*}
$$ a_k = \frac{y^{2k + 1}}{2k + 1} &&
a_k = \frac{y^{2k + 1}}{2k + 1} \thus
a_{k + 1} = \frac{y^{2k + 3}}{2k + 3} = a_k \, \frac{2k + 1}{2k + 3} \, y^2 a_{k + 1} = \frac{y^{2k + 3}}{2k + 3} = a_k \, \frac{2k + 1}{2k + 3} \, y^2
$$ \end{align*}
so, the right bound results:
hence:
$$ $$
\frac{a_{k+1}}{1 - \frac{a_{k+1}}{a_k}} \frac{a_{k+1}}{1 - \frac{a_{k+1}}{a_k}}
= \frac{a_k}{\frac{2k + 3}{2k + 1}y^{-2} - 1} = \frac{a_k}{\frac{2k + 3}{2k + 1}y^{-2} - 1}
$$ $$
and it finally follows:
from which:
$$ $$
\Delta S = a_k \left( \Delta S_k = a_k \left(
\frac{1}{\frac{2k + 3}{2k + 1}y^{-2} - 1} - \frac{y^2}{1 - y^2} \frac{1}{\frac{2k + 3}{2k + 1}y^{-2} - 1} - \frac{y^2}{1 - y^2}
\right) = \Delta S_k \with a_k = \frac{y^{2k + 1}}{2k + 1} \right)
$$ $$
By imposing $\Delta S < 10^{-D}$, where $D$ is the last correct decimal place By imposing $\Delta S_k < 10^{-D}$, where $D$ is the number decimal places
required, $k^{\text{max}}$ at which to stop can be obtained. This is achieved by required, $k^{\text{max}}$ at which to stop the summation can be obtained. This
trials, checking for every $k$ whether $\Delta S_k$ is less or greater than is achieved by trials, checking for every $k$ whether $\Delta S_k$ is less or
$10^{-D}$. greater than $10^{-D}$.
The same considerations holds for $\ln(2)$. It could be computed as the Taylor Similar considerations can be made for $\ln(2)$. The number could be computed
expansion of $\ln(1 + x)$ with $x = 1$, but it would be too slow. A better idea from the MacLaurin series of $\ln(1 + x)$ with $x = 1$, but that yields a
is to use the series with $x = -1/2$, which leads to a much more fast series: slowly convergent series with alternating sign. A trick is to use $x = -1/2$,
which leads to a much faster series with constant sign:
$$ $$
\log(2) = \sum_{k=1}^{+ \infty} \frac{1}{k \cdot 2^k} \ln(2) = \sum_{k=1}^{+ \infty} \frac{1}{k \cdot 2^k}
$$ $$
In this case the math siyplifies a lot: the ratio turns out to be $L = In this case the series ratio is found to be $L = 1/2$ and the error
1/2$ and therefore:
$$ $$
\Delta S = \frac{1}{k(k + 2) 2^{k-1}} \Delta S_k = \frac{1}{k(k + 2) 2^{k-1}}
$$ $$
Once computed the estimations of both logarithms and the series $A$, $B$ and Once the logarithms and the terms $A$, $B$ and $C$ have been computed @eq:faster
$C$ with the required precision, the program can be used to estimate $\gamma$ can finally be used to obtain $\gamma$ up to the given decimal place.
up to whatever decimal digit. To give an idea of the time it takes to compute
it, 1000 digits were computed in \SI{0.63}{s}, with a \SI{3.9}{GHz} The program was implemented with no particular care for performance but
2666mt al secondo was found to be relatively fast. On a \SI{3.9}{GHz} desktop computer with
\textcolor{red}{Scrivere specifice pc che fa i conti} \SI{2666}{MT\per\s} memory, it takes \SI{0.63}{\second} to compute the first
1000 digits. However, the quadratic complexity of the algorithm makes this
quite unpractical for computing more than a few thousands digits. For this
reason a more optimised algorithm was implemented.
### Faster implemented formula ### Optimised implementation {#sec:optimised}
When implemented, @eq:faster, which is the most fast known convergence formula, The refined Brent-McMillan formula (@eq:faster) is theoretically the fastest
is not fast at all. but is difficult to implement efficiently. In practice it turns out to be
The procedure can be accelerated by removing the $C(N)$ correction and keeping slower than the standard formula even if its asymptotic error is better.
only the two series $A(N)$ and $B(N)$. In this case, the equation can be
rewritten in a most convenient way [@brent00]: The standard formula [@brent00] ignores the correction $C(N)$ and is rewritten
in a more convenient way:
$$ $$
\gamma_k = \frac{U_k(N)}{V_k(N)} \with \gamma_k = \frac{\sum_{k = 0}^{k_\text{max}} A_k(N)}
U_k(N) = \sum_{j = 0}^k A_j(N) \et V_k(N) = \sum_{j = 0}^k B_j(N) {\sum_{k = 0}^{k_\text{max}} B_k(N)}
$$ $$
where: where:
\begin{align*}
A_k &= \frac{1}{k} \left(\frac{A_{k-1} N^2}{k} + B_k \right) & A_0 &= - \ln(N) \\
B_k &= \frac{B_{k-1} N^2}{k^2} & B_0 &= 1 \\
\end{align*}
As said, the asymptotic error of the formula decreases slower
$$ $$
A_j = \frac{1}{j} \left( \frac{A_{j-1} N^2}{j} + B_j \right) |\gamma - \gamma_k| \sim \pi e^{-4N} = \Delta_N
\with A_0 = - \ln(N)
$$ $$
and: In order to avoid the expensive computation of $\ln(N)$, N was chosen to be a
power of 2, $N = 2^p$. As before:
$$ $$
B_j = \frac{B_{j-1} N^2}{j^2} \Delta_p = \pi e^{-2^{p+2}} < 10^{-D} \thus
\with B_0 = 1 p > \log_2 \left( \ln(\pi) + D \ln(10) \right) - 2
$$ $$
This way, the error is given by: and the smaller integer is found by taking the floor:
$$ $$
|\gamma - \gamma_k| = \pi e^{-4N} p = \left\lfloor \log_2 \left( \ln(\pi) + D \ln(10) \right) \right\rfloor + 1
$$ $$
which is greater with respect to the refined formula. Then, in order to avoid The number of terms to compute, $k_\text{max}$ is still proportional to $N$
computing the logarithm of a generic number $N$, instead of using @eq:NeD, a but by a different factor, $\beta$ such that:
number $N = 2^p$ of digits were computed, where $p$ was chosen to satisfy: \begin{align*}
$$ \beta\ln(\beta) = 1 + \beta &&
\pi e^{-4N} = \pi e^{-42^p} < 10^{-D} \thus \beta = \exp(W(e^{-1}) + 1) \approx 3.59
p > \log_2 \left[ \ln(\pi) + D \ln(10) \right] -2 \end{align*}
$$
and the smaller integer greater than it was selected, namely: As regards the numbers representation, this time the GMP type `mpf_t` was
$$ employed. It consists of an arbitrary-precision floating point with a fixed
p = \left\lfloor \log_2 \left[ \ln(\pi) + D \ln(10) \right] \right\rfloor + 1 exponent of 32 or 64 bits but arbitrary mantissa length.
$$
As regards the precision of the used numbers, this time the type `pmf_t` of GMP
was emploied. It consists of an arbitrary-precision floating point whose number
of bits is to be defined when the variable is initialized.
If a number $D$ of digits is required, the number of bits can be computed from: If a number $D$ of digits is required, the number of bits can be computed from:
$$ $$
10^D = 2^b \thus b = D \log_2 (10) 10^D = 2^b \thus b = D \log_2 (10)
$$ $$
in order to the roundoff errors to not affect the final result, a security in order to avoid roundoff errors affecting the final result, a security range
range of 64 bits was added to $b$. of 64 bits was added to $b$.
Furthermore, the computation of the logarithm was made even faster by solving Furthermore, the computation of $\ln(2)$ was optimized by solving
@eq:delta instead @eq:delta analytically, instead of finding $k^{\text{max}}$ by trials.
of finding $k^{\text{max}}$ by trials. The series estimation error was
$$ $$
\Delta S = \frac{1}{k (k + 2) 2^{k - 1}} < 10^{-D} \Delta S_k = \frac{1}{k (k + 2) 2^{k - 1}} < 10^{-D}
$$ {#eq:delta} $$ {#eq:delta}
This was accomplished as follow: and the following condition should hold:
$$ $$
k (k + 2) 2^{k - 1} > 10^D \thus 2^{k - 1} [(k + 1)^2 - 1] > 10^D k (k + 2) 2^{k - 1} > 10^D \thus 2^{k - 1} [(k + 1)^2 - 1] > 10^D
$$ $$
thus, for example, the inequality is satisfied for $k = x$ such that: Asking for the equality with $k = x$ yields
$$ $$
2^{x - 1} (x + 1)^2 = 10^D 2^{x - 1} (x + 1)^2 = 10^D
$$ $$
which can be again solved by the Lambert $W$ function as follows
In order to lighten the notation: $q := 10^D$, from which:
\begin{align*} \begin{align*}
2^{x - 1} (x + 1)^2 = q 2^{x - 1} (x + 1)^2 = 10^D
&\thus \exp \left( \ln(2) \frac{x - 1}{2} \right) (x + 1) &\thus \exp \left( \ln(2) \frac{x - 1}{2} \right) (x + 1)
= \sqrt{q} \\ = 10^{D/2} \\
&\thus \exp \left( \frac{\ln(2)}{2} (x + 1)\right) e^{- \ln(2)} (x + 1) &\thus \exp \left( \frac{\ln(2)}{2} (x + 1)\right) e^{- \ln(2)} (x + 1)
= \sqrt{q} \\ = 10^{D/2} \\
&\thus \exp \left( \frac{\ln(2)}{2} (x + 1)\right) \frac{1}{2} (x + 1)
= \sqrt{q} \\
&\thus \exp \left( \frac{\ln(2)}{2} (x + 1)\right) \frac{\ln(2)}{2} (x + 1) &\thus \exp \left( \frac{\ln(2)}{2} (x + 1)\right) \frac{\ln(2)}{2} (x + 1)
= \ln (2) \sqrt{q} = \ln (2) 10^{D/2}
\end{align*} \end{align*}
Having rearranged the equation this way allows the use of the Lambert $W$ Taking $W$ both sides of the equation:
function:
$$ $$
W(x \cdot e^x) = x \frac{\ln(2)}{2} (x + 1) = W ( \ln(2) 10^{D/2} )
$$ $$
Hence, by computing the Lambert W function of both sides of the equation:
$$
\frac{\ln(2)}{2} (x + 1) = W ( \ln(2) \sqrt{q} )
$$
which leads to: which leads to:
$$ $$
x = \frac{2}{\ln(2)} \, W ( \ln(2) \sqrt{q} ) - 1 x = \frac{2}{\ln(2)} \, W ( \ln(2) 10^{D/2} ) - 1
$$ $$
Finally, the smallest integer greater than $x$ can be found as: Finally, the smallest integer greater than $x$ can be found as:
$$ $$
k^{\text{max}} = k^{\text{max}} =
\left\lfloor \frac{2}{\ln(2)} \, W ( \ln(2) \sqrt{q} ) \right\rfloor \left\lfloor \frac{2}{\ln(2)} \, W (\ln(2) 10^{D/2}) \right\rfloor
\with 10^D
$$ $$
When a large number $D$ of digits is to be computed (but it also works for few When a large number $D$ of digits is to be computed, the exponentiation can
digits), the Lambert W function can be approximated by the first five terms of easily overflow if working in double precision. However, $W$ grows like a
its asymptotic value [@corless96]: logarithm, so intuitively this operation could be optimized out.
In fact this can be done by using the asymptotic expansion [@corless96]
at large $x$ of $W(x)$ :
$$ $$
W(x) = L_1 - L_2 + \frac{L_2}{L_1} + \frac{L_2 (L_2 - 2)}{2 L_1^2} W(x) = L_1 - L_2 + \frac{L_2}{L_1} + \frac{L_2 (L_2 - 2)}{2 L_1^2}
+ \frac{L_2 (36L_2 - 22L_2^2 + 3L_2^3 -12)}{12L_1^4} + \frac{L_2 (3L_2^3 - 22L_2^2 + 36L_2 - 12)}{12L_1^4}
+ O\left(\left(\frac{L_2}{L_1}\right)^5\right)
$$ $$
where where
$$ $$
L_1 = \ln(x) \et L_2 = \ln(\ln(x)) L_1 = \ln(x) \et L_2 = \ln(\ln(x))
$$ $$
\textcolor{red}{Tempi del risultato e specifiche del pc usato} Now the usual properties of the logarithm can be applied to write:
\begin{align*}
L_1 &= \ln(\ln(2) 10^{D/2}) = \ln(\ln(2)) + \frac{D}{2}\ln(10) \\
L_2 &= \ln(L_1)
\end{align*}
In this way $k^{\text{max}}$ can be easily computed with standard double
precision floating points, with no risk of overflow or need of arbitrary
precision arithmetic.
### 1 million digits computation
On the same desktop computer, the optimised program computes the first 1000
digits of $\gamma$ in \SI{4.6}{\milli\second}: a $137\times$ improvement over
the previous program. This make it suitable for a more intensive computation.
As a proof of concept, the program was then used to compute a million digits of
$\gamma$, which took \SI{1.26}{\hour} of running time.
The result was verified in \SI{6.33}{\hour} by comparing it with the
`mpmath.euler()` function from the mpmath [@mpmath13] multiple-precision
library with the Python integer backend. The library also uses the standard
Brent-McMillan algorithm and variable precision floating points but is not
based on GMP arithmetic types.

View File

@ -1,4 +0,0 @@
- aggiungere tempi di calcolo della γ
simpy 1mln cifre: 6h 20'
fast 1mld cifre: 1h 37'