Euler-Mascheroni Constant

Admin 1
Viewed 907

In programming, I used only Integers. But this time for some calculations. I need to calculate Euler-Mascheroni Constant γ . up to n-th decimal.{Though n ∈ [30, 150] is enough for me.

enter image description here

  • [x] = gif(x) = math.floor(x)

But, I doubt the precision Numerical Algorithm

I need higher degree of accuracy using Python.

floating-accuracy precision math algorithm python
2 Answers

From the French Wikipedia discussion page, an approximation to 6 decimal places:

import math as m
EulerMascheroniApp = round( (1.-m.gamma(1+1.e-8))*1.e14 )*1.e-6
print(EulerMascheroniApp)
# 0.577216

This constant is also available in the sympy module, under the name EulerGamma:

>>> import sympy
>>> sympy.EulerGamma
EulerGamma
>>> sympy.EulerGamma.evalf()
0.577215664901533
>>> - sympy.polygamma(0,1)
EulerGamma
>>> sympy.stieltjes(0)
EulerGamma
>>> sympy.stieltjes(0, 1)
EulerGamma

Documentation:

On this last documentation link, you can find more information about how to evaluate the constant with more precision, if the default of .evalf() is not enough.

If you still want to compute the constant yourself as an exercise, I suggest comparing your results to sympy's constant, to check for accuracy and correctness.

user_30807230

You can calculate it using python Decimal built-in module to control how many decimals (https://docs.python.org/2/library/decimal.html) you are going to use.

a = 1/7
len(str(a))-2
Out[1] 17

using Decimal:

from decimal import *
getcontext().prec = 90 #90 decimals precision
a = Decimal(1) / Decimal(7)
len(str(a))-2
Out[2] 90

basically:

n = 100000
Euler_Mascheroni = -Decimal(log(Decimal(n))) + sum([Decimal(1)/Decimal(i) for i in range(1,n)])
Euler_Mascheroni

Out[3] Decimal('0.577210664893199330073570099082905499710324918344701101627529415938181982282214')

finally, you can "arbitrarily" increase precision:

from decimal import *
from math import log

def Euler_Mascheroni(n,nth_decimals = 80):
    getcontext().prec = nth_decimals
    
    SUM = Decimal(0)
    for i in range(1,n):
        SUM+=Decimal(1)/Decimal(i)
    
    return -Decimal(log(Decimal(n))) + SUM
Euler_Mascheroni(100000000,nth_decimals = 120)

which gives:

Decimal('0.5772156599015311156682000509495086978690376512201034388184221886576113026091829254475798266636558124658249350393045066')

Answering comment from @Stef

EM = Decimal(0.57721566490153286060651209008240243104215933593992)#reference taken from wikipedia

n = 100000000
Decimal(log(Decimal(n)))

getcontext().prec = 100
SUM = Decimal(0)
for i in range(1,n):
    SUM+=Decimal(1)/Decimal(i)

EM - (SUM-Decimal(log(Decimal(n))))

will give

Decimal('5.00000174 ... 85E-9')
user_111458200
Related