414 lines
14 KiB
Markdown
414 lines
14 KiB
Markdown
# Exercise 2
|
||
|
||
|
||
## Euler-Mascheroni constant
|
||
|
||
The Euler-Mascheroni constant is defined as the limiting difference between the
|
||
partial sums of the harmonic series and the natural logarithm:
|
||
$$
|
||
\gamma = \lim_{n \rightarrow +\infty} \left( \sum_{k=1}^{n} \frac{1}{k}
|
||
- \ln(n) \right)
|
||
$$ {#eq:gamma}
|
||
|
||
and represents the limiting blue area in @fig:gamma. The first 30 digits of
|
||
$\gamma$ are:
|
||
$$
|
||
\gamma = 0.57721\ 56649\ 01532\ 86060\ 65120\ 90082 \dots
|
||
$$ {#eq:exact}
|
||
|
||
In complex analysis, this constant is related to many functions and can be
|
||
evaluated through a variety of identities. In this work, five methods were
|
||
implemented and their results discussed. In fact, evaluating $\gamma$ with a
|
||
limited precision due to floating points number representation entails limited
|
||
precision on the estimation of $\gamma$ itself due to roundoff errors. All the
|
||
known methods involve sums, subtractions or products of very big or small
|
||
numbers, packed in series, partial sums or infinite products. Thus, the
|
||
efficiency of the methods lies on how quickly they converge to their limit.
|
||
|
||
![The area of the blue region converges to the Euler–Mascheroni
|
||
constant.](images/2-gamma-area.png){#fig:gamma width=7cm}
|
||
|
||
|
||
## Computing the constant
|
||
|
||
|
||
### Definition
|
||
|
||
First, in order to give a quantitative idea of how hard it is to reach a good
|
||
estimation of $\gamma$, it was naively computed using the definition given in
|
||
@eq:gamma. The difference was computed for increasing value of $n$, with
|
||
$n_{i+1} = 10 \cdot n_i$ and $n_1 = 20$, till the approximation starts getting
|
||
worse, namely:
|
||
$$
|
||
| \gamma(n_{i+1}) - \gamma | > | \gamma(n_i) - \gamma|
|
||
$$
|
||
|
||
and $\gamma (n_i)$ was selected as the result (see @tbl:1_results).
|
||
|
||
-----------------------------------------------
|
||
n sum $|\gamma(n)-\gamma|$
|
||
----------- ------------- ---------------------
|
||
\SI{2e1}{} \SI{2.48e-02}{}
|
||
|
||
\SI{2e2}{} \SI{2.50e-03}{}
|
||
|
||
\SI{2e3}{} \SI{2.50e-04}{}
|
||
|
||
\SI{2e4}{} \SI{2.50e-05}{}
|
||
|
||
\SI{2e5}{} \SI{2.50e-06}{}
|
||
|
||
\SI{2e6}{} \SI{2.50e-07}{}
|
||
|
||
\SI{2e7}{} \SI{2.50e-08}{}
|
||
|
||
\SI{2e8}{} \SI{2.50e-09}{}
|
||
|
||
\SI{2e9}{} \SI{2.55e-10}{}
|
||
|
||
\SI{2e10}{} \SI{2.42e-11}{}
|
||
|
||
\SI{2e11}{} \SI{1.44e-08}{}
|
||
-----------------------------------------------
|
||
|
||
Table: Partial results using the definition of $\gamma$ with double
|
||
precision. {#tbl:1_results}
|
||
|
||
The convergence is logarithmic: to fix the first $d$ decimal places, about
|
||
$10^d$ terms are needed. The double precision runs out at the
|
||
10\textsuperscript{th} place, $n=\SI{2e10}{}$.
|
||
Since all the number are given with double precision, there can be at best 15
|
||
correct digits but only 10 were correctly computed: this means that when the
|
||
terms of the series start being smaller than the smallest representable double,
|
||
the sum of all the remaining terms give a number $\propto 10^{-11}$.
|
||
Best result in @tbl:first.
|
||
|
||
--------- -----------------------
|
||
true: 0.57721\ 56649\ 01533
|
||
|
||
approx: 0.57721\ 56648\ 77325
|
||
|
||
diff: 0.00000\ 00000\ 24207
|
||
--------- -----------------------
|
||
|
||
Table: First method best result. From the top down: true value, best estimation
|
||
and difference between them. {#tbl:first}
|
||
|
||
|
||
### Alternative formula
|
||
|
||
As a first alternative, the constant was computed through the identity which
|
||
relates $\gamma$ to the $\Gamma$ function as follow [@davis59]:
|
||
$$
|
||
\gamma = \lim_{M \rightarrow + \infty} \sum_{k = 1}^{M}
|
||
\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
|
||
@tbl:second). It went sour: the convergence is worse than using the definition
|
||
itself. Only two places were correctly computed (#@tbl:second).
|
||
|
||
--------- -----------------------
|
||
true: 0.57721\ 56649\ 01533
|
||
|
||
approx: 0.57225\ 72410\ 34058
|
||
|
||
diff: 0.00495\ 84238\ 67473
|
||
--------- -----------------------
|
||
|
||
Table: Best estimation of $\gamma$ using
|
||
the alternative formula. {#tbl:second}
|
||
|
||
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
|
||
represented. Furthermore, the convergence is slowed down by the logarithmic
|
||
factor.
|
||
|
||
|
||
### Reciprocal $\Gamma$ function
|
||
|
||
A better result was found using the well known reciprocal $\Gamma$ function
|
||
formula [@bak91]:
|
||
$$
|
||
\frac{1}{\Gamma(z)} = z e^{yz} \prod_{k = 1}^{+ \infty}
|
||
\left( 1 + \frac{z}{k} \right) e^{-z/k}
|
||
$$
|
||
|
||
which gives:
|
||
$$
|
||
\gamma = - \frac{1}{z} \ln \left( z \Gamma(z)
|
||
\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
|
||
of the infinite product (it happens for $k = 456565794 \sim \SI{4.6e8}{}$,
|
||
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}
|
||
= 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 |$
|
||
----- ------------------------ ------ ------------------------
|
||
1 \SI{9.712e-9}{} 8.95 \SI{9.770e-9}{}
|
||
|
||
3 \SI{9.320e-9}{} 8.96 \SI{9.833e-9}{}
|
||
|
||
5 \SI{9.239e-9}{} 8.97 \SI{9.622e-9}{}
|
||
|
||
7 \SI{9.391e-9}{} 8.98 \SI{9.300e-9}{}
|
||
|
||
9 \SI{8.482e-9}{} 8.99 \SI{9.059e-9}{}
|
||
|
||
11 \SI{9.185e-9}{} 9.00 \SI{8.482e-9}{}
|
||
|
||
13 \SI{9.758e-9}{} 9.01 \SI{9.564e-9}{}
|
||
|
||
15 \SI{9.747e-9}{} 9.02 \SI{9.260e-9}{}
|
||
|
||
17 \SI{9.971e-9}{} 9.03 \SI{9.264e-9}{}
|
||
|
||
19 \SI{10.084e-9}{} 9.04 \SI{9.419e-9}{}
|
||
---------------------------------------------------------------
|
||
|
||
Table: Differences between some obtained values of $\gamma$ and
|
||
the exact one found with the reciprocal $\Gamma$ function formula.
|
||
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
|
||
one ($z = 9.00$) are listed. {#tbl:3_results}
|
||
|
||
As can be seen in @tbl:3_results, the best value for $z$ is only by chance,
|
||
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.
|
||
|
||
--------- -----------------------
|
||
true: 0.57721\ 56649\ 01533
|
||
|
||
approx: 0.57721\ 56564\ 18607
|
||
|
||
diff: 0.00000\ 00084\ 82925
|
||
--------- -----------------------
|
||
|
||
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
|
||
|
||
The fastest known convergence belongs to the following formula [@yee19]:
|
||
$$
|
||
\gamma = \frac{A(N)}{B(N)} -\frac{C(N)}{B^2(N)} - \ln(N)
|
||
$$ {#eq:faster}
|
||
|
||
with:
|
||
\begin{align*}
|
||
&A(N) = \sum_{k=1}^{+ \infty} \frac{N^k}{k!} \cdot H(k)
|
||
\with H(k) = \sum_{j=1}^{k} \frac{1}{j} \\
|
||
&B(N) = \sum_{k=1}^{+ \infty} \frac{N^k}{k!} \\
|
||
&C(N) = \frac{1}{4N} \sum_{k=0}^{2N}
|
||
\frac{((2k)!)^3}{(k!)^4 \cdot (16k)^2k} \\
|
||
\end{align*}
|
||
|
||
The series $A$ and $B$ were computed till there is no difference between two
|
||
consecutive terms. The number of desired correct decimals $D$ was given in
|
||
input and $N$ was consequently computed through the formula:
|
||
$$
|
||
N = \left\lfloor 2 + \frac{1}{4} \cdot \ln(10) \cdot D \right\rfloor
|
||
$$ {#eq:NeD}
|
||
|
||
given in [@brent00]. Results are shown in @tbl:fourth.
|
||
Due to roundoff errors, the best results was obtained for $N = 10$. Up to 15
|
||
places were correctly computed.
|
||
|
||
|
||
--------- ------------------------------
|
||
true: 0.57721\ 56649\ 01532\ 75452
|
||
|
||
approx: 0.57721\ 56649\ 01532\ 86554
|
||
|
||
diff: 0.00000\ 00000\ 00000\ 11102
|
||
--------- ------------------------------
|
||
|
||
Table: $\gamma$ estimation with the fastest
|
||
known convergence formula (@eq:faster). {#tbl:fourth}
|
||
|
||
|
||
### Arbitrary precision
|
||
|
||
To overcome the issues related to the double representation, one can resort to a
|
||
representation with arbitrary precision. In the GMP library (which stands for
|
||
GNU Multiple Precision arithmetic), for example, real numbers can be
|
||
approximated by a ratio of two integers (fraction) with arbitrary precision:
|
||
this means a check is performed on every operation and in case of an integer
|
||
overflow, additional memory is requested and used to represent the larger result.
|
||
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
|
||
precision. Thus, a program that computes the Euler-Mascheroni constant within
|
||
a user controllable precision was implemented.
|
||
|
||
According to [@brent00], the $A(N)$, $B(N)$ and $C(N)$ series awere computed up
|
||
to $k_{\text{max}} = 5N$, which was verified to compute the correct digits up to
|
||
500 decimal places.
|
||
The GMP library offers functions to perform some operations such as addition,
|
||
multiplication, division, etc. However, the logarithm function is not
|
||
implemented. Thus, most of the code carries out the $\ln(N)$ computation.
|
||
First, it should be noted that the logarithm of only some special numbers can
|
||
be computed with arbitrary precision, namely the ones of which a converging
|
||
series is known. 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)
|
||
$$
|
||
|
||
Since a fast converging series for $\ln(2)$ is known (it will be shwn shortly),
|
||
$b = 2$ was chosen. As well as for the scientific notation, in order to get the
|
||
mantissa $1 \leqslant N_0 < 2$, the number of binary digits of $N$ must be
|
||
computed (conveniently, a dedicated function `mpz_sizeinbase()` can be found in
|
||
GMP). If the digits are
|
||
$n$:
|
||
$$
|
||
e = n - 1 \thus N = N_0 \cdot 2^{n-1} \thus N_0 = \frac{N}{2^{n - 1}}
|
||
$$
|
||
|
||
The logarithm od whichever number $N_0$ can be computed using the series of
|
||
$\text{atanh}(x)$, which converges for $|x| < 1$:
|
||
$$
|
||
\text{atanh}(x) = \sum_{k = 0}^{+ \infty} \frac{x^{2k + 1}}{2x + 1}
|
||
$$
|
||
|
||
In fact:
|
||
$$
|
||
\text{tanh}(x) = \frac{e^x - e^{-x}}{e^x + e^{-x}}
|
||
= \frac{e^{2x} - 1}{e^{2x} + 1}
|
||
$$
|
||
|
||
with the change of variable $z = e^{2x}$, one gets:
|
||
$$
|
||
\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)
|
||
$$
|
||
|
||
Then, by defining:
|
||
$$
|
||
y = \frac{z - 1}{z + 1} \thus z = \frac{1 + y}{1 - y}
|
||
$$
|
||
|
||
from which:
|
||
$$
|
||
\ln \left( \frac{1 + y}{1 - y} \right) = 2 \, \text{atanh}(y)
|
||
= 2 \, \sum_{k = 0}^{+ \infty} \frac{y^{2k + 1}}{2y + 1}
|
||
$$
|
||
|
||
which is true if $|y| < 1$. The idea is to set:
|
||
$$
|
||
N_0 = \frac{1 + y}{1 - y} \thus y = \frac{N_0 - 1}{N_0 + 1} < 1
|
||
$$
|
||
|
||
and therefore:
|
||
$$
|
||
\ln(N_0) = \ln \left( \frac{1 + y}{1 - y} \right) =
|
||
2 \sum_{k = 0}^{+ \infty} \frac{y^{2k + 1}}{2k + 1}
|
||
$$
|
||
|
||
But when to stop computing the series?
|
||
Given a partial sum $S_k$ of the series, it is possible to know when a digit is
|
||
definitely correct. The key lies in the following concept [@riddle08]. Letting
|
||
$S$ be the value of the series:
|
||
$$
|
||
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
|
||
ratio of the series terms, which must be $< 1$ in order for it to converge (in
|
||
this case, it is easy to prove that $L = y^2$). The width $\Delta S$ of the
|
||
interval containing $S$ gives the precision of the estimate $\tilde{S}$ if this
|
||
last is assumed to be the middle value of it, namely:
|
||
$$
|
||
\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}}
|
||
\right) \et
|
||
\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
|
||
evaluated. With a bit of math:
|
||
$$
|
||
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
|
||
$$
|
||
|
||
hence:
|
||
$$
|
||
\frac{a_{k+1}}{1 - \frac{a_{k+1}}{a_k}}
|
||
= \frac{a_k}{\frac{2k + 3}{2k + 1}y^{-2} - 1}
|
||
$$
|
||
|
||
from which:
|
||
$$
|
||
\Delta S = a_k \left(
|
||
\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}
|
||
$$
|
||
|
||
By imposing $\Delta S < 1/10^D$, where $D$ is the correct decimal place required,
|
||
$k$ at which to stop can be obtained. This is achieved by trials, checking for
|
||
every $k$ whether $\Delta S_k$ is less or greater than $1/10^D$.
|
||
|
||
The same considerations holds for $\ln(2)$. It could be computed as the Taylor
|
||
expansion of $\ln(1 + x)$ with $x = 1$, but it would be too slow. A better idea
|
||
is to use the series with $x = -1/2$, which leads to a much more fast series:
|
||
$$
|
||
\log(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 =
|
||
1/2$ and therefore:
|
||
$$
|
||
\Delta S = \frac{1}{k(k + 2) 2^{k-1}}
|
||
$$
|
||
|
||
Once computed the estimations of both logarithms and the series $A$, $B$ and
|
||
$C$ with the required precision, the program can be used to estimate $\gamma$
|
||
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}
|
||
2666mt al secondo
|
||
\textcolor{red}{Scrivere specifice pc che fa i conti}
|
||
|
||
### Faster implemented formula
|
||
|
||
When implemented, @eq:faster, which is the most fast known convergence formula,
|
||
is not fast at all.
|
||
The procedure can be accelerated by removing the $C(N)$ correction and keeping
|
||
only the two series $A(N)$ and $B(N)$. In this case, the equation can be
|
||
rewritten in a most convenient way [@brent00]:
|
||
$$
|
||
\gamma_k = \frac{U_k(N)}{V_k(N)} \with
|
||
U_k(N) = \sum_{j = 0}^k A_j(N) \et V_k(N) = \sum_{j = 0}^k B_j(N)
|
||
$$
|
||
|
||
where:
|
||
$$
|
||
A_j = \frac{1}{j} \left( \frac{A_{j-1} N^2}{j} + B_j \right)
|
||
\with A_0 = - \ln(N)
|
||
$$
|
||
|
||
and:
|
||
$$
|
||
B_j = \frac{B_{j-1} N^2}{j^2}
|
||
\with B_0 = - \ln(N)
|
||
$$
|
||
|
||
This way, the error is given by:
|
||
$$
|
||
|\gamma - \gamma_k| = \pi e^{4N}
|
||
$$
|
||
|
||
Then, in order to avoid computing the logarithm of a generic number $N$,
|
||
instead of using @eq:NeD, a number $N = 2^p$ of digits were computed, where $p$
|
||
was chosen in order to satisfy:
|
||
$$
|
||
|
||
$$
|