Mathematical constants

Mpmath supports arbitrary-precision computation of various common (and less common) mathematical constants. These constants are implemented as lazy objects that can evaluate to any precision. Whenever the objects are used as function arguments or as operands in arithmetic operations, they automagically evaluate to the current working precision. A lazy number can be converted to a regular mpf using the unary + operator, or by calling it as a function:

>>> from mpmath import *
>>> mp.dps = 15
>>> pi
<pi: 3.14159~>
>>> 2*pi
mpf('6.2831853071795862')
>>> +pi
mpf('3.1415926535897931')
>>> pi()
mpf('3.1415926535897931')
>>> mp.dps = 40
>>> pi
<pi: 3.14159~>
>>> 2*pi
mpf('6.283185307179586476925286766559005768394338')
>>> +pi
mpf('3.141592653589793238462643383279502884197169')
>>> pi()
mpf('3.141592653589793238462643383279502884197169')

Exact constants

The predefined objects j (imaginary unit), inf (positive infinity) and nan (not-a-number) are shortcuts to mpc and mpf instances with these fixed values.

Pi (pi)

mp.pi

\pi, roughly equal to 3.141592654, represents the area of the unit circle, the half-period of trigonometric functions, and many other things in mathematics.

Mpmath can evaluate \pi to arbitrary precision:

>>> from mpmath import *
>>> mp.dps = 50; mp.pretty = True
>>> +pi
3.1415926535897932384626433832795028841971693993751

This shows digits 99991-100000 of \pi:

>>> mp.dps = 100000
>>> str(pi)[-10:]
'5549362464'

Possible issues

pi always rounds to the nearest floating-point number when used. This means that exact mathematical identities involving \pi will generally not be preserved in floating-point arithmetic. In particular, multiples of pi (except for the trivial case 0*pi) are not the exact roots of sin(), but differ roughly by the current epsilon:

>>> mp.dps = 15
>>> sin(pi)
1.22464679914735e-16

One solution is to use the sinpi() function instead:

>>> sinpi(1)
0.0

See the documentation of trigonometric functions for additional details.

Degree (degree)

mp.degree

Represents one degree of angle, 1^{\circ} = \pi/180, or about 0.01745329. This constant may be evaluated to arbitrary precision:

>>> from mpmath import *
>>> mp.dps = 50; mp.pretty = True
>>> +degree
0.017453292519943295769236907684886127134428718885417

The degree object is convenient for conversion to radians:

>>> sin(30 * degree)
0.5
>>> asin(0.5) / degree
30.0

Base of the natural logarithm (e)

mp.e

The transcendental number e = 2.718281828... is the base of the natural logarithm (ln()) and of the exponential function (exp()).

Mpmath can be evaluate e to arbitrary precision:

>>> from mpmath import *
>>> mp.dps = 50; mp.pretty = True
>>> +e
2.7182818284590452353602874713526624977572470937

This shows digits 99991-100000 of e:

>>> mp.dps = 100000
>>> str(e)[-10:]
'2100427165'

Possible issues

e always rounds to the nearest floating-point number when used, and mathematical identities involving e may not hold in floating-point arithmetic. For example, ln(e) might not evaluate exactly to 1.

In particular, don’t use e**x to compute the exponential function. Use exp(x) instead; this is both faster and more accurate.

Golden ratio (phi)

mp.phi

Represents the golden ratio \phi = (1+\sqrt 5)/2, approximately equal to 1.6180339887. To high precision, its value is:

>>> from mpmath import *
>>> mp.dps = 50; mp.pretty = True
>>> +phi
1.6180339887498948482045868343656381177203091798058

Formulas for the golden ratio include the following:

>>> (1+sqrt(5))/2
1.6180339887498948482045868343656381177203091798058
>>> findroot(lambda x: x**2-x-1, 1)
1.6180339887498948482045868343656381177203091798058
>>> limit(lambda n: fib(n+1)/fib(n), inf)
1.6180339887498948482045868343656381177203091798058

Euler’s constant (euler)

mp.euler

Euler’s constant or the Euler-Mascheroni constant \gamma = 0.57721566... is a number of central importance to number theory and special functions. It is defined as the limit

\gamma = \lim_{n\to\infty} H_n - \log n

where H_n = 1 + \frac{1}{2} + \ldots + \frac{1}{n} is a harmonic number (see harmonic()).

Evaluation of \gamma is supported at arbitrary precision:

>>> from mpmath import *
>>> mp.dps = 50; mp.pretty = True
>>> +euler
0.57721566490153286060651209008240243104215933593992

We can also compute \gamma directly from the definition, although this is less efficient:

>>> limit(lambda n: harmonic(n)-log(n), inf)
0.57721566490153286060651209008240243104215933593992

This shows digits 9991-10000 of \gamma:

>>> mp.dps = 10000
>>> str(euler)[-10:]
'4679858165'

Integrals, series, and representations for \gamma in terms of special functions include the following (there are many others):

>>> mp.dps = 25
>>> -quad(lambda x: exp(-x)*log(x), [0,inf])
0.5772156649015328606065121
>>> quad(lambda x,y: (x-1)/(1-x*y)/log(x*y), [0,1], [0,1])
0.5772156649015328606065121
>>> nsum(lambda k: 1/k-log(1+1/k), [1,inf])
0.5772156649015328606065121
>>> nsum(lambda k: (-1)**k*zeta(k)/k, [2,inf])
0.5772156649015328606065121
>>> -diff(gamma, 1)
0.5772156649015328606065121
>>> limit(lambda x: 1/x-gamma(x), 0)
0.5772156649015328606065121
>>> limit(lambda x: zeta(x)-1/(x-1), 1)
0.5772156649015328606065121
>>> (log(2*pi*nprod(lambda n:
...     exp(-2+2/n)*(1+2/n)**n, [1,inf]))-3)/2
0.5772156649015328606065121

For generalizations of the identities \gamma = -\Gamma'(1) and \gamma = \lim_{x\to1} \zeta(x)-1/(x-1), see psi() and stieltjes() respectively.

Catalan’s constant (catalan)

mp.catalan

Catalan’s constant K = 0.91596559... is given by the infinite series

K = \sum_{k=0}^{\infty} \frac{(-1)^k}{(2k+1)^2}.

Mpmath can evaluate it to arbitrary precision:

>>> from mpmath import *
>>> mp.dps = 50; mp.pretty = True
>>> +catalan
0.91596559417721901505460351493238411077414937428167

One can also compute K directly from the definition, although this is significantly less efficient:

>>> nsum(lambda k: (-1)**k/(2*k+1)**2, [0, inf])
0.91596559417721901505460351493238411077414937428167

This shows digits 9991-10000 of K:

>>> mp.dps = 10000
>>> str(catalan)[-10:]
'9537871503'

Catalan’s constant has numerous integral representations:

>>> mp.dps = 50
>>> quad(lambda x: -log(x)/(1+x**2), [0, 1])
0.91596559417721901505460351493238411077414937428167
>>> quad(lambda x: atan(x)/x, [0, 1])
0.91596559417721901505460351493238411077414937428167
>>> quad(lambda x: ellipk(x**2)/2, [0, 1])
0.91596559417721901505460351493238411077414937428167
>>> quad(lambda x,y: 1/(1+(x*y)**2), [0, 1], [0, 1])
0.91596559417721901505460351493238411077414937428167

As well as series representations:

>>> pi*log(sqrt(3)+2)/8 + 3*nsum(lambda n:
...  (fac(n)/(2*n+1))**2/fac(2*n), [0, inf])/8
0.91596559417721901505460351493238411077414937428167
>>> 1-nsum(lambda n: n*zeta(2*n+1)/16**n, [1,inf])
0.91596559417721901505460351493238411077414937428167

Apery’s constant (apery)

mp.apery

Represents Apery’s constant, which is the irrational number approximately equal to 1.2020569 given by

\zeta(3) = \sum_{k=1}^\infty\frac{1}{k^3}.

The calculation is based on an efficient hypergeometric series. To 50 decimal places, the value is given by:

>>> from mpmath import *
>>> mp.dps = 50; mp.pretty = True
>>> +apery
1.2020569031595942853997381615114499907649862923405

Other ways to evaluate Apery’s constant using mpmath include:

>>> zeta(3)
1.2020569031595942853997381615114499907649862923405
>>> -psi(2,1)/2
1.2020569031595942853997381615114499907649862923405
>>> 8*nsum(lambda k: 1/(2*k+1)**3, [0,inf])/7
1.2020569031595942853997381615114499907649862923405
>>> f = lambda k: 2/k**3/(exp(2*pi*k)-1)
>>> 7*pi**3/180 - nsum(f, [1,inf])
1.2020569031595942853997381615114499907649862923405

This shows digits 9991-10000 of Apery’s constant:

>>> mp.dps = 10000
>>> str(apery)[-10:]
'3189504235'

Khinchin’s constant (khinchin)

mp.khinchin

Khinchin’s constant K = 2.68542... is a number that appears in the theory of continued fractions. Mpmath can evaluate it to arbitrary precision:

>>> from mpmath import *
>>> mp.dps = 50; mp.pretty = True
>>> +khinchin
2.6854520010653064453097148354817956938203822939945

An integral representation is:

>>> I = quad(lambda x: log((1-x**2)/sincpi(x))/x/(1+x), [0, 1])
>>> 2*exp(1/log(2)*I)
2.6854520010653064453097148354817956938203822939945

The computation of khinchin is based on an efficient implementation of the following series:

>>> f = lambda n: (zeta(2*n)-1)/n*sum((-1)**(k+1)/mpf(k)
...     for k in range(1,2*int(n)))
>>> exp(nsum(f, [1,inf])/log(2))
2.6854520010653064453097148354817956938203822939945

Glaisher’s constant (glaisher)

mp.glaisher

Glaisher’s constant A, also known as the Glaisher-Kinkelin constant, is a number approximately equal to 1.282427129 that sometimes appears in formulas related to gamma and zeta functions. It is also related to the Barnes G-function (see barnesg()).

The constant is defined as A = \exp(1/12-\zeta'(-1)) where \zeta'(s) denotes the derivative of the Riemann zeta function (see zeta()).

Mpmath can evaluate Glaisher’s constant to arbitrary precision:

>>> from mpmath import *
>>> mp.dps = 50; mp.pretty = True
>>> +glaisher
1.282427129100622636875342568869791727767688927325

We can verify that the value computed by glaisher is correct using mpmath’s facilities for numerical differentiation and arbitrary evaluation of the zeta function:

>>> exp(mpf(1)/12 - diff(zeta, -1))
1.282427129100622636875342568869791727767688927325

Here is an example of an integral that can be evaluated in terms of Glaisher’s constant:

>>> mp.dps = 15
>>> quad(lambda x: log(gamma(x)), [1, 1.5])
-0.0428537406502909
>>> -0.5 - 7*log(2)/24 + log(pi)/4 + 3*log(glaisher)/2
-0.042853740650291

Mpmath computes Glaisher’s constant by applying Euler-Maclaurin summation to a slowly convergent series. The implementation is reasonably efficient up to about 10,000 digits. See the source code for additional details.

References: http://mathworld.wolfram.com/Glaisher-KinkelinConstant.html

Mertens constant (mertens)

mp.mertens

Represents the Mertens or Meissel-Mertens constant, which is the prime number analog of Euler’s constant:

B_1 = \lim_{N\to\infty}
    \left(\sum_{p_k \le N} \frac{1}{p_k} - \log \log N \right)

Here p_k denotes the k-th prime number. Other names for this constant include the Hadamard-de la Vallee-Poussin constant or the prime reciprocal constant.

The following gives the Mertens constant to 50 digits:

>>> from mpmath import *
>>> mp.dps = 50; mp.pretty = True
>>> +mertens
0.2614972128476427837554268386086958590515666482612

References: http://mathworld.wolfram.com/MertensConstant.html

Twin prime constant (twinprime)

mp.twinprime

Represents the twin prime constant, which is the factor C_2 featuring in the Hardy-Littlewood conjecture for the growth of the twin prime counting function,

\pi_2(n) \sim 2 C_2 \frac{n}{\log^2 n}.

It is given by the product over primes

C_2 = \prod_{p\ge3} \frac{p(p-2)}{(p-1)^2} \approx 0.66016

Computing C_2 to 50 digits:

>>> from mpmath import *
>>> mp.dps = 50; mp.pretty = True
>>> +twinprime
0.66016181584686957392781211001455577843262336028473

References: http://mathworld.wolfram.com/TwinPrimesConstant.html