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')
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, 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.
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
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.
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 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 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
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
>>> -diff(trigamma, 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 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*n))
>>> exp(nsum(f, [1,inf])/log(2))
2.6854520010653064453097148354817956938203822939945
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
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
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