analistica/notes/sections/2.md

# 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/gamma-area.png){#fig:gamma width=7cm}

## Computing the constant

### Definition

First, in order to have 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}$.


--------- -----------------------
true:	     0.57721\ 56649\ 01533

approx:    0.57721\ 56648\ 77325

diff:	     0.00000\ 00000\ 24207
--------- -----------------------

Table: First method results. {#tbl:first}

### Alternative formula

As a first alternative, the constant was computed through the identity which
relates $\gamma$ to the $\Gamma$ function as follow:

$$
  \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.

--------- -----------------------
true:	     0.57721\ 56649\ 01533

approx:    0.57225\ 72410\ 34058

diff:	     0.00495\ 84238\ 67473
--------- -----------------------

Table: Second method results. {#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 (even if this is not a series
and, consequently, it is not properly a "convergence") is slowed down by
the logarithmic factor.

### Reciprocal $\Gamma$ function

A better result was found using the well known reciprocal $\Gamma$ function
formula:

$$
  \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$, meaning that for this
value of $k$, the term of the product is equal to 1). 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$. As can be seen in @tbl:3_results, it's 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.

---------------------------------------------------------------
    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 the obtained values of $\gamma$ and the exact
one. {#tbl:3_results}

--------- -----------------------
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:  
(source: http://www.numberworld.org/y-cruncher/internals/formulas.html):

$$
  \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$ are computed till there is no difference between two
consecutive terms.  
The number of desired correct decimals is given in input and $N$ is
consequently computed through a formula given in the same article
above-mentioned. Results are shown in @tbl:fourth.

--------- ------------------------------
true:	     0.57721\ 56649\ 01532\ 75452

approx:    0.57721\ 56649\ 01532\ 86554

diff:	     0.00000\ 00000\ 00000\ 11102
--------- ------------------------------

Table: Fourth method results. {#tbl:fourth}

Due to roundoff errors, the best results is obtained for $N = 10$. Since up to
15 places were correctly computed, an approximation of $\gamma$ better than the
one reached with the definition in @eq:gamma was obtained.

### 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 has been implemented. Unlike the previously
mentioned programs, this one was more carefully optimized.

The $A$ and $B$ series are computed up to an arbitrary limit $k_{\text{max}}$.
Different values of $k_{\text{max}}$ were tested but, obviously, they all
eventually reach a point where the approximation cannot guarantee the
requested precision; the solution turned out to be to let $k_{\text{max}}$ depends on
$N$.  Consequently, $k_{\text{max}}$ was chosen to be $5N$, after it has been
verified to produce 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 $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 = b - 1 \thus N_0 = \frac{N}{2^{n - 1}}
$$

Then, by defining:

$$
  N_0 = \frac{1 + y}{1 - y} \thus y = \frac{N_0 - 1}{N_0 + 1} < 1
$$

and the following series (which is convergent for $y < 1$) can therefore be
used:

$$
  \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, meaning that no matter how large $k$ can be,
it will not affect that decimal place. The key lies in the following concept.
Letting $S$ the value of the series:

$$
  S_k + \frac{a_{k+1}}{1 - \frac{a_{k+1}}{a_k}} < S < S_k + a_k \frac{L}{1 -L}
$$

where $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$, for a given $S_k$, gives
the precision of the partial sum with respect to $S$. By imposing $\Delta S <
1/10^D$ where $D$ is the correct decimal place required, the desired precision
is obtained.  This is achieved by trials, checking at every step whether $\Delta
S$ is less or greater than $1/10^D$.

The same holds for $\ln(2)$, which is given by the following series:

$$
  \log(2) = \sum_{k=1}^{+ \infty} \frac{1}{k \cdot 2^k}
$$

In this case the ratio is $L = 1/2$.