Mpmath basics

In interactive code examples that follow, it will be assumed that the main contents of the mpmath package have been imported with “import *“.

>>> from mpmath import *

Number types

Mpmath provides the following numerical types:

Class Description
mpf Real float
mpc Complex float
mpi Real interval
matrix Matrix



>>> mpf(4)
>>> mpf(2.5)
>>> mpf("1.25e6")
>>> mpf(mpf(2))
>>> mpf("inf")


>>> mpc(2,3)
>>> mpc(complex(2,3)).imag


>>> mpc(1j)**0.5

Setting the precision


>>> print mp
Mpmath settings:
  mp.prec = 53                [default: 53]
  mp.dps = 15                 [default: 15]
  mp.trap_complex = False     [default: False]



>>> mp.dps = 100
>>> mp.dps
>>> mp.prec


>>> mp.dps = 50
>>> mpf(1) / 6
>>> mp.dps = 25


>>> mp.dps = 10
>>> mpc(1,2) / 3


>>> mp.dps = 15

Or why not 1 googolplex:



Temporarily changing the precision


>>> mp.prec += 2
>>> # do_something()
>>> mp.prec -= 2


>>> from __future__ import with_statement
>>> with workdps(20):  
...     print mpf(1)/7
...     with extradps(10):
...         print mpf(1)/7
>>> mp.dps



>>> @workdps(6)
... def f():
...     return mpf(1)/3
>>> f()


>>> mp.dps = 15
>>> print exp(1)


>>> mp2 = mp.clone()
>>> mp.dps = 10
>>> mp2.dps = 20
>>> print mp.mpf(1) / 3
>>> print mp2.mpf(1) / 3


Providing correct input


>>> mp.dps = 30
>>> mpf(10.9)   # bad
>>> mpf('10.9')  # good
>>> mp.dps = 15






>>> mpf("3.14159")
>>> print mpf("3.14159")
>>> print mpc(1j)**0.5
(0.707106781186548 + 0.707106781186548j)


>>> mp.pretty = True
>>> mpf(0.6)
>>> mp.pretty = False
>>> mpf(0.6)

The number of digits with which numbers are printed by default is determined by the working precision. To specify the number of digits to show without changing the working precision, use the nstr and nprint functions:

>>> a = mpf(1) / 6
>>> a
>>> nstr(a, 8)
>>> nprint(a, 8)
>>> nstr(a, 50)

Fast low-precision arithmetic

Many of the algorithms used in mpmath work perfectly well with ordinary Python float and complex numbers, which use double precision (53 bits). Although the global functions (which are methods of the mp object) always convert inputs to mpmath numbers, the alternative fp object instead converts them to float or complex, and uses internal functions designed for double precision. When large amounts of function evaluations (numerical integration, plotting, etc) are required, and when fp arithmetic provides sufficient accuracy, this can give a significant speedup over mp arithmetic.

To take advantage of this feature, simply use fp.func instead of func or mp.func:

>>> u = fp.erfc(2.5)
>>> print u
>>> type(u)
<type 'float'>
>>> print mp.erfc(2.5)
>>> fp.matrix([[1,2],[3,4]]) ** 2
[['7.0', '10.0'],
 ['15.0', '22.0']])
>>> type(_[0,0])
<type 'float'>
>>> print fp.quad(fp.sin, [0, fp.pi])    # numerical integration

Due to intermediate rounding and cancellation errors, results computed with fp arithmetic may be much less accurate than those computed with mp using an equivalent precision (mp.prec = 53), since the latter often uses increased internal precision. The accuracy is highly problem-dependent: for some functions, fp almost always gives 14-15 correct digits; for others, results can be accurate to only 2-3 digits or even completely wrong. The recommended use for fp is therefore to speed up large-scale computations where accuracy can be verified in advance on a subset of the input set, or where results can be verified afterwards.


If matplotlib is available, the functions plot and cplot in mpmath can be used to plot functions respectively as x-y graphs and in the complex plane. Also, splot can be used to produce 3D surface plots.


Output of plot([cos, sin], [-4, 4])

mpmath.plot(ctx, f, xlim=[, -5, 5], ylim=None, points=200, file=None, dpi=None, singularities=[], axes=None)

Shows a simple 2D plot of a function f(x) or list of functions [f_0(x), f_1(x), \ldots, f_n(x)] over a given interval specified by xlim. Some examples:

plot(lambda x: exp(x)*li(x), [1, 4])
plot([fresnels, fresnelc], [-4, 4])
plot([sqrt, cbrt], [-4, 4])
plot(lambda t: zeta(0.5+t*j), [-20, 20])
plot([floor, ceil, abs, sign], [-5, 5])

Points where the function raises a numerical exception or returns an infinite value are removed from the graph. Singularities can also be excluded explicitly as follows (useful for removing erroneous vertical lines):

plot(cot, ylim=[-5, 5])   # bad
plot(cot, ylim=[-5, 5], singularities=[-pi, 0, pi])  # good

For parts where the function assumes complex values, the real part is plotted with dashes and the imaginary part is plotted with dots.


This function requires matplotlib (pylab).


Output of fp.cplot(fp.gamma, points=100000)

mpmath.cplot(ctx, f, re=[, -5, 5], im=[, -5, 5], points=2000, color=None, verbose=False, file=None, dpi=None, axes=None)

Plots the given complex-valued function f over a rectangular part of the complex plane specified by the pairs of intervals re and im. For example:

cplot(lambda z: z, [-2, 2], [-10, 10])
cplot(zeta, [0, 1], [0, 50])

By default, the complex argument (phase) is shown as color (hue) and the magnitude is show as brightness. You can also supply a custom color function (color). This function should take a complex number as input and return an RGB 3-tuple containing floats in the range 0.0-1.0.

To obtain a sharp image, the number of points may need to be increased to 100,000 or thereabout. Since evaluating the function that many times is likely to be slow, the ‘verbose’ option is useful to display progress.


This function requires matplotlib (pylab).


Output of splot for the donut example.

mpmath.splot(ctx, f, u=[, -5, 5], v=[, -5, 5], points=100, keep_aspect=True, wireframe=False, file=None, dpi=None, axes=None)

Plots the surface defined by f.

If f returns a single component, then this plots the surface defined by z = f(x,y) over the rectangular domain with x = u and y = v.

If f returns three components, then this plots the parametric surface x, y, z = f(u,v) over the pairs of intervals u and v.

For example, to plot a simple function:

>>> from mpmath import *
>>> f = lambda x, y: sin(x+y)*cos(y)
>>> splot(f, [-pi,pi], [-pi,pi])    

Plotting a donut:

>>> r, R = 1, 2.5
>>> f = lambda u, v: [r*cos(u), (R+r*sin(u))*cos(v), (R+r*sin(u))*sin(v)]
>>> splot(f, [0, 2*pi], [0, 2*pi])    


This function requires matplotlib (pylab) or higher.

Utility functions

The following convenience functions are provided to simplify common tasks. (Mathematical functions with more specific applications are covered in later sections of the documentation.)

Conversion and utilities for numbers

mpmath.mpmathify(x, strings=True)

Converts x to an mpf or mpc. If x is of type mpf, mpc, int, float, complex, the conversion will be performed losslessly.

If x is a string, the result will be rounded to the present working precision. Strings representing fractions or complex numbers are permitted.

>>> from mpmath import *
>>> mp.dps = 15; mp.pretty = False
>>> mpmathify(3.5)
>>> mpmathify('2.1')
>>> mpmathify('3/4')
>>> mpmathify('2+3j')
mpmath.almosteq(s, t, rel_eps=None, abs_eps=None)

Determine whether the difference between s and t is smaller than a given epsilon, either relatively or absolutely.

Both a maximum relative difference and a maximum difference (‘epsilons’) may be specified. The absolute difference is defined as |s-t| and the relative difference is defined as |s-t|/\max(|s|, |t|).

If only one epsilon is given, both are set to the same value. If none is given, both epsilons are set to 2^{-p+m} where p is the current working precision and m is a small integer. The default setting typically allows almosteq() to be used to check for mathematical equality in the presence of small rounding errors.


>>> from mpmath import *
>>> mp.dps = 15
>>> almosteq(3.141592653589793, 3.141592653589790)
>>> almosteq(3.141592653589793, 3.141592653589700)
>>> almosteq(3.141592653589793, 3.141592653589700, 1e-10)
>>> almosteq(1e-20, 2e-20)
>>> almosteq(1e-20, 2e-20, rel_eps=0, abs_eps=0)

Return True if x is integer-valued; otherwise return False:

>>> from mpmath import *
>>> isint(3)
>>> isint(mpf(3))
>>> isint(3.2)
>>> isint(inf)

Optionally, Gaussian integers can be checked for:

>>> isint(3+0j)
>>> isint(3+2j)
>>> isint(3+2j, gaussian=True)

Return True if the absolute value of x is infinite; otherwise return False:

>>> from mpmath import *
>>> isinf(inf)
>>> isinf(-inf)
>>> isinf(3)
>>> isinf(3+4j)
>>> isinf(mpc(3,inf))
>>> isinf(mpc(inf,3))

Return True if x is a NaN (not-a-number), or for a complex number, whether either the real or complex part is NaN; otherwise return False:

>>> from mpmath import *
>>> isnan(3.14)
>>> isnan(nan)
>>> isnan(mpc(3.14,2.72))
>>> isnan(mpc(3.14,nan))


mpmath.fsum(terms, absolute=False, squared=False)

Calculates a sum containing a finite number of terms (for infinite series, see nsum()). The terms will be converted to mpmath numbers. For len(terms) > 2, this function is generally faster and produces more accurate results than the builtin Python function sum().

>>> from mpmath import *
>>> mp.dps = 15; mp.pretty = False
>>> fsum([1, 2, 0.5, 7])

With squared=True each term is squared, and with absolute=True the absolute value of each term is used.


Calculates a product containing a finite number of factors (for infinite products, see nprod()). The factors will be converted to mpmath numbers.

>>> from mpmath import *
>>> mp.dps = 15; mp.pretty = False
>>> fprod([1, 2, 0.5, 7])
mpmath.fdot(A, B=None)

Computes the dot product of the iterables A and B,

\sum_{k=0} A_k B_k.

Alternatively, fdot() accepts a single iterable of pairs. In other words, fdot(A,B) and fdot(zip(A,B)) are equivalent. The elements are automatically converted to mpmath numbers.

With conjugate=True, the elements in the second vector will be conjugated:

\sum_{k=0} A_k \overline{B_k}


>>> from mpmath import *
>>> mp.dps = 15; mp.pretty = False
>>> A = [2, 1.5, 3]
>>> B = [1, -1, 2]
>>> fdot(A, B)
>>> zip(A, B)
[(2, 1), (1.5, -1), (3, 2)]
>>> fdot(_)
>>> A = [2, 1.5, 3j]
>>> B = [1+j, 3, -1-j]
>>> fdot(A, B)
mpc(real='9.5', imag='-1.0')
>>> fdot(A, B, conjugate=True)
mpc(real='3.5', imag='-5.0')

Number generation

mpmath.fraction(p, q)

Given Python integers (p, q), returns a lazy mpf representing the fraction p/q. The value is updated with the precision.

>>> from mpmath import *
>>> mp.dps = 15
>>> a = fraction(1,100)
>>> b = mpf(1)/100
>>> print a; print b
>>> mp.dps = 30
>>> print a; print b      # a will be accurate
>>> mp.dps = 15
Returns an mpf with value chosen randomly from [0, 1). The number of randomly generated bits in the mantissa is equal to the working precision.

This is a generalized version of Python’s range() function that accepts fractional endpoints and step sizes and returns a list of mpf instances. Like range(), arange() can be called with 1, 2 or 3 arguments:

[0, 1, 2, \ldots, x]
arange(a, b)
[a, a+1, a+2, \ldots, x]
arange(a, b, h)
[a, a+h, a+h, \ldots, x]

where b-1 \le x < b (in the third case, b-h \le x < b).

Like Python’s range(), the endpoint is not included. To produce ranges where the endpoint is included, linspace() is more convenient.


>>> from mpmath import *
>>> mp.dps = 15; mp.pretty = False
>>> arange(4)
[mpf('0.0'), mpf('1.0'), mpf('2.0'), mpf('3.0')]
>>> arange(1, 2, 0.25)
[mpf('1.0'), mpf('1.25'), mpf('1.5'), mpf('1.75')]
>>> arange(1, -1, -0.75)
[mpf('1.0'), mpf('0.25'), mpf('-0.5')]
mpmath.linspace(*args, **kwargs)

linspace(a, b, n) returns a list of n evenly spaced samples from a to b. The syntax linspace(mpi(a,b), n) is also valid.

This function is often more convenient than arange() for partitioning an interval into subintervals, since the endpoint is included:

>>> from mpmath import *
>>> mp.dps = 15; mp.pretty = False
>>> linspace(1, 4, 4)
[mpf('1.0'), mpf('2.0'), mpf('3.0'), mpf('4.0')]

You may also provide the keyword argument endpoint=False:

>>> linspace(1, 4, 4, endpoint=False)
[mpf('1.0'), mpf('1.75'), mpf('2.5'), mpf('3.25')]

Rounding and printing

mpmath.chop(x, tol=None)

Chops off small real or imaginary parts, or converts numbers close to zero to exact zeros. The input can be a single number or an iterable:

>>> from mpmath import *
>>> mp.dps = 15; mp.pretty = False
>>> chop(5+1e-10j, tol=1e-9)
>>> nprint(chop([1.0, 1e-20, 3+1e-18j, -4, 2]))
[1.0, 0.0, 3.0, -4.0, 2.0]

The tolerance defaults to 100*eps.

mpmath.nstr(x, n=6, **kwargs)

Convert an mpf or mpc to a decimal string literal with n significant digits. The small default value for n is chosen to make this function useful for printing collections of numbers (lists, matrices, etc).

If x is a list or tuple, nstr() is applied recursively to each element. For unrecognized classes, nstr() simply returns str(x).

The companion function nprint() prints the result instead of returning it.

>>> from mpmath import *
>>> nstr([+pi, ldexp(1,-500)])
'[3.14159, 3.05494e-151]'
>>> nprint([+pi, ldexp(1,-500)])
[3.14159, 3.05494e-151]
mpmath.nprint(x, n=6, **kwargs)
Equivalent to print nstr(x, n).


mpmath.monitor(f, input='print', output='print')

Returns a wrapped copy of f that monitors evaluation by calling input with every input (args, kwargs) passed to f and output with every value returned from f. The default action (specify using the special string value 'print') is to print inputs and outputs to stdout, along with the total evaluation count:

>>> from mpmath import *
>>> mp.dps = 5; mp.pretty = False
>>> diff(monitor(exp), 1)   # diff will eval f(x-h) and f(x+h)
in  0 (mpf('0.99999999906867742538452148'),) {}
out 0 mpf('2.7182818259274480055282064')
in  1 (mpf('1.0000000009313225746154785'),) {}
out 1 mpf('2.7182818309906424675501024')

To disable either the input or the output handler, you may pass None as argument.

Custom input and output handlers may be used e.g. to store results for later analysis:

>>> mp.dps = 15
>>> input = []
>>> output = []
>>> findroot(monitor(sin, input.append, output.append), 3.0)
>>> len(input)  # Count number of evaluations
>>> print input[3], output[3]
((mpf('3.1415076583334066'),), {}) 8.49952562843408e-5
>>> print input[4], output[4]
((mpf('3.1415928201669122'),), {}) -1.66577118985331e-7
mpmath.timing(f, *args, **kwargs)

Returns time elapsed for evaluating f(). Optionally arguments may be passed to time the execution of f(*args, **kwargs).

If the first call is very quick, f is called repeatedly and the best time is returned.