# Calcuration of Euler Mascheroni constant using Euler Maclaurin Summation Formula package for Maxima CAS¶

Written by Yasuaki dot Honda at gmail dot com

GPL2.0

Application to Euler Mascheroni constant is taken from the following PDF:

http://people.csail.mit.edu/kuat/courses/euler-maclaurin.pdf

In [1]:
install_github("YasuakiHonda","euler-maclaurin-sum","master");

Out[1]:
$\tag{{\it \%o}_{1}}\left[ \mbox{ /root/quicklisp/dists/quicklisp/archives/euler-maclaurin-sum-master.gz } , \mbox{ /root/quicklisp/local-projects/YasuakiHonda-euler-maclaurin-sum-69eb859/ } \right]$
In [2]:
asdf_load_source("euler-maclaurin-sum");

Out[2]:
$\tag{{\it \%o}_{2}}\#$
In [3]:
ems;

Out[3]:
$\tag{{\it \%o}_{3}}\sum_{n=N}^{M}{f\left(n\right)}=\sum_{k=1}^{K-1}{\frac{\left(-1\right)^{k+1}\,B_{k+1}\,\left(\left.\frac{d^{k}}{d\,x^{k}}\,f\left(x\right)\right|_{x=M}-\left.\frac{d^{k}}{d\,x^{k}}\,f\left(x\right)\right|_{x=N}\right)}{\left(k+1\right)!}}+\frac{\left(-1\right)^{K+1}\,\int_{N}^{M}{\overline{B}_{K}\left(x\right)\,\left(\frac{d^{K}}{d\,x^{K}}\,f\left(x\right)\right)\;dx}}{K!}+\int_{N}^{M}{f\left(x\right)\;dx}+\frac{f\left(N\right)}{2}+\frac{f\left(M\right)}{2}$
In [4]:
assume(M>1);
ems,f(x):=1/x,N:1,K:1;

Out[4]:
$\tag{{\it \%o}_{4}}\left[ M>1 \right]$
Out[4]:
$\tag{{\it \%o}_{5}}\sum_{n=1}^{M}{\frac{1}{n}}=\int_{1}^{M}{\frac{d}{d\,x}\,\left(\frac{1}{x}\right)\,\overline{B}_{1}\left(x\right)\;dx}+\int_{1}^{M}{\frac{1}{x}\;dx}+\frac{1}{2\,M}+\frac{1}{2}$
In [5]:
%,nouns;

Out[5]:
$\tag{{\it \%o}_{6}}\sum_{n=1}^{M}{\frac{1}{n}}=-\int_{1}^{M}{\frac{\overline{B}_{1}\left(x\right)}{x^2}\;dx}+\log M+\frac{1}{2\,M}+\frac{1}{2}$
In [6]:
EuGamma:%-log(M);

Out[6]:
$\tag{{\it \%o}_{7}}\sum_{n=1}^{M}{\frac{1}{n}}-\log M=-\int_{1}^{M}{\frac{\overline{B}_{1}\left(x\right)}{x^2}\;dx}+\frac{1}{2\,M}+\frac{1}{2}$
In [7]:
plot2d([bernpoly(x-floor(x),1),bernpoly(x-floor(x),1)/x^2,1/x^2],[x,1,5])\$

In [8]:
apply1(first(rhs(EuGamma)), emsRemInt);

Out[8]:
$\tag{{\it \%o}_{9}}-\sum_{k=1}^{M-1}{\left(\log \left(k+1\right)-\log k+\frac{1}{2\,k+2}+\frac{k}{k+1}-\frac{1}{2\,k}-1\right)}$
In [9]:
rem:%,M:100,nouns,numer;

Out[9]:
$\tag{{\it \%o}_{10}}0.0772073316515306$
In [10]:
limit(rest(rhs(EuGamma)),M,inf)+rem,numer;

Out[10]:
$\tag{{\it \%o}_{11}}0.5772073316515306$
In [11]:
%gamma,numer;

Out[11]:
$\tag{{\it \%o}_{12}}0.5772156649015329$
In [ ]:


In [12]:
assume(M-10000>0);
ems,f(x):=1/x,N:10000,K:6,nouns;

Out[12]:
$\tag{{\it \%o}_{13}}\left[ M>10000 \right]$
Out[12]:
$\tag{{\it \%o}_{14}}\sum_{n=10000}^{M}{\frac{1}{n}}=-\int_{10000}^{M}{\frac{\overline{B}_{6}\left(x\right)}{x^7}\;dx}+\log M+\frac{1}{2\,M}+\frac{\frac{1}{100000000}-\frac{1}{M^2}}{12}-\frac{\frac{3}{5000000000000000}-\frac{6}{M^4}}{720}+\frac{\frac{3}{25000000000000000000000}-\frac{120}{M^6}}{30240}-\log 10000+\frac{1}{20000}$
In [13]:
EuGamma:%-log(M);

Out[13]:
$\tag{{\it \%o}_{15}}\sum_{n=10000}^{M}{\frac{1}{n}}-\log M=-\int_{10000}^{M}{\frac{\overline{B}_{6}\left(x\right)}{x^7}\;dx}+\frac{1}{2\,M}+\frac{\frac{1}{100000000}-\frac{1}{M^2}}{12}-\frac{\frac{3}{5000000000000000}-\frac{6}{M^4}}{720}+\frac{\frac{3}{25000000000000000000000}-\frac{120}{M^6}}{30240}-\log 10000+\frac{1}{20000}$
In [14]:
diff(bernpoly(x,6),x),factor;

Out[14]:
$\tag{{\it \%o}_{16}}\left(x-1\right)\,x\,\left(2\,x-1\right)\,\left(3\,x^2-3\,x-1\right)$
In [15]:
[bernpoly(0,6),bernpoly(1/2,6),bernpoly(1,6)];

Out[15]:
$\tag{{\it \%o}_{17}}\left[ \frac{1}{42} , -\frac{31}{1344} , \frac{1}{42} \right]$
In [16]:
[err:integrate(%[1]/x^7,x,10000,inf),ev(err,numer)];

Out[16]:
$\tag{{\it \%o}_{18}}\left[ \frac{1}{252000000000000000000000000} , 3.968253968253969 \times 10^{-27} \right]$
In [17]:
fpprec:30;

Out[17]:
$\tag{{\it \%o}_{19}}30$
In [18]:
bfloat(sum(1/n,n,1,9999)+limit(rest(rhs(EuGamma)),M,inf));

Out[18]:
$\tag{{\it \%o}_{20}}5.77215664901532860606512090082_B \times 10^{-1}$
In [19]:
bfloat(%gamma);

Out[19]:
$\tag{{\it \%o}_{21}}5.77215664901532860606512090082_B \times 10^{-1}$