mpmath
¶Mpmath
is a Python library for arbitrary-precision floating-point arithmetic. For general information about mpmath
, see the project website.
On top of its support for arbitrary-precision arithmetic, mpmath
provides extensive support for transcendental functions, evaluation of sums, integrals, limits, roots, and so on. It also does many standard mathematical tasks like,
from math import sqrt
print(sqrt(2))
1.4142135623730951
import numpy as np
print(np.sqrt(2))
1.4142135623730951
print(np.sqrt(2, dtype=np.float32))
1.4142135
import mpmath.libmp
mpmath.libmp.BACKEND
'gmpy'
import mpmath as mp
from mpmath import *
mp.dps=10
print("Sqrt of 2 (10 decimal places): ",mpf(2) ** mpf('0.5'))
print("Twice of pi (10 decimal places)", 2*pi)
Sqrt of 2 (10 decimal places): 1.414213562 Twice of pi (10 decimal places) 6.283185307
mp.dps=25
print("Sqrt of 2 (10 decimal places): ",mpf(2) ** mpf('0.5'))
print("Twice of pi (10 decimal places)", 2*pi)
Sqrt of 2 (10 decimal places): 1.414213562373095048801689 Twice of pi (10 decimal places) 6.283185307179586476925287
mp
context shows the current settings¶print(mp)
Mpmath settings: mp.prec = 86 [default: 53] mp.dps = 25 [default: 15] mp.trap_complex = False [default: False]
mpf
instances instead of regular float¶An mpf instance holds a real-valued floating-point number. They work analogously to Python floats, but support arbitrary-precision arithmetic.
mpf(12.3)
mpf('12.30000000000000071054273576')
mpf('12.3')
mpf('12.29999999999999999999999996')
mp.pretty = True
x=mpf('12.3')
print(x)
12.3
x**mpf(0.5)
3.507135583350036383363493
12.3**0.5
3.5071355833500366
%%timeit -n10 -r100
fac(10**5)
15.9 µs ± 2.79 µs per loop (mean ± std. dev. of 100 runs, 10 loops each)
from math import factorial
%%timeit -n10 -r100
factorial(10**5)
170 ms ± 11.7 ms per loop (mean ± std. dev. of 100 runs, 10 loops each)
f_mpmath = fac(10**5)
f_python = factorial(10**5)
print(f_mpmath)
print(float(f_python))
2.824229407960347874293422e+456573
--------------------------------------------------------------------------- OverflowError Traceback (most recent call last) <ipython-input-36-1475943975fc> in <module> 2 f_python = factorial(10**5) 3 print(f_mpmath) ----> 4 print(float(f_python)) OverflowError: int too large to convert to float
mpmathify('3/4')
0.75
mpmathify('3/4') + mpmathify('2/7')
1.035714285714285714285714
mpmathify('2+3j')
(2.0 + 3.0j)
mpmathify('2+3j')**(mpf('0.5'))
(1.674149228035540040448039 + 0.8959774761298381247157338j)
import numpy as np
A = np.random.normal(size=100)
B = np.random.normal(size=100)
%%timeit -n10 -r100
A.dot(B)
The slowest run took 96.61 times longer than the fastest. This could mean that an intermediate result is being cached. 2.01 µs ± 6.12 µs per loop (mean ± std. dev. of 100 runs, 10 loops each)
%%timeit -n10 -r100
fdot(A,B)
443 µs ± 159 µs per loop (mean ± std. dev. of 100 runs, 10 loops each)
plot([cos, sin], [-4, 4])
plot([fresnels, fresnelc], [-4, 4])
cplot(lambda z: z, [-2, 2], [-10, 10])
cplot(gamma, points=10000)
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, float(2*pi)], [0, float(2*pi)],points=200,dpi=120)
for r in [root(3+4j, 7, k) for k in range(7)]:
print("{}".format(r))
(1.247472705895526221523927 + 0.1662271241773525444423184j) (0.6478249113010029918267916 + 1.078954351705588446708222j) (-0.4396482547230982620878282 + 1.179206945741718717440614j) (-1.196057317750688078827009 + 0.3914926581963046418346134j) (-1.051810825389031340038652 - 0.691023585965792972199872j) (-0.1155293284786681638014006 - 1.253184975583352195227716j) (0.9077481091449566314041713 - 0.8716725182718191829981798j)
%%timeit -n5 -r20
for n in range(25,35):
[binomial(mpf('50'),mpf(str(k))) for k in range(n+1)]
13.4 ms ± 3.11 ms per loop (mean ± std. dev. of 20 runs, 5 loops each)
import scipy.special
binomial(10**20, 10**20-5)
8.3333333333333333325e+97
Scipy
cannot do calculation with these large arguments¶scipy.special.binom(int(10**20), int(10**20-5))
--------------------------------------------------------------------------- TypeError Traceback (most recent call last) <ipython-input-128-aacea9fae810> in <module> ----> 1 scipy.special.binom(int(10**20), int(10**20-5)) TypeError: ufunc 'binom' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''
binomial(2.5,1.2)
2.585287397240184524055334
binomial(mpf('-12/5'),mpf('2/3'))
-0.4150035108034928041587776
scipy.special.binom(-2.4,0.667)
-0.417154501641134
mpmath
¶%%timeit -n5 -r20
hyperfac(200)
341 µs ± 135 µs per loop (mean ± std. dev. of 20 runs, 5 loops each)
%%timeit -n5 -r20
prod=1
for n in range(1,201):
prod*=n**n
11.5 ms ± 827 µs per loop (mean ± std. dev. of 20 runs, 5 loops each)
hyperfac(200)
1.142560024245405712989778e+41908
mp.dps=25
fib(20)
6765.0
fib(20**5)
1.256941107429143065410115e+668760
fib(pi)
2.117027057916099911684598
fib(3+4j)
(-5248.511307283720182781954 - 14195.96228835301088580928j)
polyval([3, 0, -2, 5, 0, 0, 0,-3], 0.27)
-2.975983920803909993806701
polyroots([3, 0, -2, 5, 0, 0, 0,-3])
[-1.295821155626593414551295, -0.8973346781104871960058213, 0.8700356165162543117430376, (-0.1021871544563401270795588 - 0.8141076161287989904493385j), (-0.1021871544563401270795588 + 0.8141076161287989904493385j), (0.7637472630667532764865984 - 0.9407317375480711196626355j), (0.7637472630667532764865984 + 0.9407317375480711196626355j)]
def func1(x):
return exp(-0.1*x)*sin(x**2)
plot(func1,[-2,5])
for x0 in range(-2,6):
r = findroot(func1,x0=x0,tol=1e-10,solver="newton")
print(f"Solution near {x0}: {r}")
Solution near -2: -1.772453850905516027298168 Solution near -1: -0.0000005430976634644979218348307 Solution near 0: 0.0 Solution near 1: 0.0000002920122051525101216309755 Solution near 2: 1.772453850905516027298167 Solution near 3: 3.069980123839465465438655 Solution near 4: 3.963327297606011013345029 Solution near 5: 5.01325654926200100483153
print(diff(lambda x: x**2 + x, 1.0))
print(diff(lambda x: x**2 + x, 1.0, 2))
print(diff(lambda x: x**2 + x, 1.0, 3))
3.0 2.0 0.0
diff(lambda x,y: 3*x*y + 2*y - x, (0.25, 0.5), (0,1))
2.75
$\int_{0}^{\infty}\frac{1}{1+x^2}=\frac{\pi}{2}$
$\int_{-\infty}^{\infty}\text{exp}(-x^2)=\sqrt{\pi}$
mp.dps=50
quad(lambda x: 1/(x**2+1), [0, inf])
1.5707963267948966192313216916397514420985846996876
quad(lambda x: exp(-x**2), [-inf, inf])
1.7724538509055160272981674833411451827975494561224
$\int_{0}^{1}\int_{0}^{1}\frac{x-1}{(1-xy).\text{log}(xy)}$
mp.dps = 30
f = lambda x, y: (x-1)/((1-x*y)*log(x*y))
quad(f, [0, 1], [0, 1])
0.577215664901532860606512090082
$y'(x)=y(x), y(0)=1$
$\rightarrow y(x)=\text{exp}(x)$
f = odefun(lambda x, y: y, 0, 1)
f(0.0)
1.0
f(1)
2.71828182845904523536028747135
$y'(x)=y(x), y(0)=\pi/2$
$\rightarrow y(x)=2\text{tan}^{-1}\text{exp}(x^2/2)$
f = odefun(lambda x, y: x*sin(y), 0, pi/2)
f(1)
2.0511774059286260678362079816