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

Written by Yasuaki dot Honda at gmail dot com

Copyright 2020 Yasuaki Honda

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])$
Gnuplot Produced by GNUPLOT 5.2 patchlevel 8 x (-floor(x))+x-1/2 (-floor(x))+x-1/2 ((-floor(x))+x-1/2)/x2 ((-floor(x))+x-1/2)/x2 1/x2 1/x2 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1 1.5 2 2.5 3 3.5 4 4.5 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}\]