todo: updated

notes/docs/bibliography.bib

@@ -169,3 +169,13 @@
   year={1963},
   institution={Stanford Researhc INST Menlo Park CA}
 }
+
+@incollection{brent00,
+  title={Some new algorithms for high-precision computation of Euler’s
+         constant},
+  author={Brent, Richard P and McMillan, Edwin M},
+  booktitle={Pi: A Source Book},
+  pages={448--455},
+  year={2000},
+  publisher={Springer}
+}

notes/sections/2.md

@@ -214,11 +214,13 @@ 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 = \text{floor} \left( 2 + \frac{1}{4} \cdot \ln(10) \cdot D \right)
-$$
+  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.
 
-given in [@yee19], where floor returns the highest integer smaller than its
-argument. Results are shown in @tbl:fourth.
 
 --------- ------------------------------
 true:	     0.57721\ 56649\ 01532\ 75452
@@ -231,9 +233,6 @@ diff:	     0.00000\ 00000\ 00000\ 11102
 Table: $\gamma$ estimation with the fastest
 known convergence formula (@eq:faster). {#tbl:fourth}
 
-Due to roundoff errors, the best results was obtained for $N = 10$. Up to 15
-places were correctly computed.
-
 
 ### Arbitrary precision
 
@@ -248,71 +247,167 @@ 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. 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
-was verified to produce the correct digits up to 500 decimal places.
+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, $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
+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_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 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
 $$
 
-the following series, which is convergent for $y < 1$, can therefore be used:
-
+and therefore:
 $$
-  \ln \left( \frac{1 + y}{1 - y} \right) =
+  \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. Letting $S$ the value
-of the series [@riddle08]:
-
+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 $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$.
+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}
+$$
 
-The same holds for $\ln(2)$, which is given by the following series:
+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 ratio is $L = 1/2$.
+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:
+$$
+  
+$$

notes/todo

@@ -1,6 +1,6 @@
 - riscrivere il 2
 - riscrivere il readme del 6
 - riscrivere il readme del 7
-- incolonnare bene ex-5
+- aggiungere tempi di calcolo della γ
 
 On the Lambert W function, formula 4.19 Corless