Euler-Mascheroni Constant

Question

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.

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

But, I doubt the precision Numerical Algorithm

I need higher degree of accuracy using Python.

Oneill

Since it is a constant, if you need it for some other calculations why not just include the first couple hundred digits in your source code? On the other hand, if you need to approximate it yourself (e.g. as a homework project) then the decimal module is the way to go. – Heroin

@JohnColeman (so you teach maths and gave me Homework 💖) Thanks! the difficulty is the convergence. – Oneill

Although the constant is probably an irrational number, you may also want to consider using Python's fractions module if you want to compute the value yourself. This would alleviate the need to pick a specific decimal precision. – Knudson

I am too lazy to write the Python code, but the site from Xavier Gourdon and Pascal Sebah seems fantastic. Specifically the 6th paragraph of this page gives C code for more than 1600 decimals... – Morbid

@DarshanPatil: How is “the convergence” any difficulty in including the first 150 digits in your source code? Including the first 150 digits in your source code means you simply have to get the first 150 digits from any reference source, or calculate it once with special-purpose software such as Maple, and then you are done. You would not have to deal with “the convergence,” whatever you mean by that. – Mundy

@SergeBallesta Thanks! This website is like a pond for people with an unquenchable thirst for knowledge in number theory. – Oneill

Answer 1

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:

• math.gamma;
• sympy.EulerGamma;
• sympy.functions.special;
• sympy: numerical evaluation.

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.

Crabbe

how can the precision in EulerMascheroniApp = round( (1.-m.gamma(1+1.e-8))*1.e14 )*1.e-6, can be increased? – K

I tried playing with the round function, and gained two places, but not more: g = 0.57721566490153286060651209008240243104215933593992 g - round( (1.-m.gamma(1+1.e-8))*1.e14 )*1.e-6 # -3.350984670857926e-07 g - round( (1.-m.gamma(1+1.e-8))*1.e15 )*1.e-7 # 6.490153292570966e-08 g - round( (1.-m.gamma(1+1.e-8))*1.e16 )*1.e-8 # 3.4901532885989184e-08 g - round( (1.-m.gamma(1+1.e-8))*1.e17 )*1.e-9 # 3.390153280324881e-08 – Crabbe

@Crabbe Finally I installed sympy and can accept your answer😃. – Oneill

Answer 2

If you are using numpy, you can use numpy.euler_gamma:

>>> import numpy as np
>>> np.euler_gamma
0.5772156649015329

Sauers

Answer 3

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')

K

Your result differs from sympy.EulerGamma.evalf() starting at the 6th decimal: print(sympy.EulerGamma.evalf() - Euler_Mascheroni) results in 5.00000833358882e-6. I checked by comparing with the 50-places value given at Wikipedia and it also differs starting at the 6th decimal. So, it looks like tuning the number of decimals of module Decimal is not enough to fix the approximation error. – Crabbe

In fact, the approximation round( (1.-m.gamma(1+1.e-8))*1.e14 )*1.e-6 is more accurate than -Decimal(log(Decimal(n))) + sum([Decimal(1)/Decimal(i) for i in range(1,n)]) – Crabbe

Also, which log are you using? I used math.log to try your code; your only import is from decimal import *, but module decimal doesn't have a log. I suggest using import decimal as d rather than from decimal import *, so that readers can more easily figure out which functions your are using. – Crabbe

You are right, more decimal points alone don't solve the problem, we need more iterations. A edited my answer with the improvement. Also, i used log from math. And the question is note about getting the actual value of the constant but o calculate it, using arbitrarily (30 to 80 decimals) more precision. – K

You're using 100000000 iterations to get 5 decimal places of precision. That's a good sign that the method is not efficient, and a better method is required. Using Decimal instead of float will only be necessary when the precision of your method is better than the precision of float, which is far from being the situation at hand. – Crabbe

it is 9 digits of precision, you are comparing it to a 6 digit precision you give, not a true reference. If you use the one from wikipedia: 0.57721566490153286060651209008240243104215933593992, you will get the error to the right order of magnitude. Although, it is really not efficient, but can arbitrarily increase precision, which is the question. – K