# Exponential integrals and error functions

Exponential integrals give closed-form solutions to a large class of commonly occurring transcendental integrals that cannot be evaluated using elementary functions. Integrals of this type include those with an integrand of the form \(t^a e^{t}\) or \(e^{-x^2}\), the latter giving rise to the Gaussian (or normal) probability distribution.

The most general function in this section is the incomplete gamma function, to which all others can be reduced. The incomplete gamma function, in turn, can be expressed using hypergeometric functions (see Hypergeometric functions).

## Incomplete gamma functions

mpmath.gammainc(z, a=0, b=None, regularized=False)

  `gammainc(z, a=0, b=inf)` computes the (generalized) incomplete gamma function with integration limits \([a, b]\):

  \[\Gamma(z,a,b) = \int_a^b t^{z-1} e^{-t} \, dt\]

  The generalized incomplete gamma function reduces to the following special cases when one or both endpoints are fixed:

  - \(\Gamma(z,0,\infty)\) is the standard (“complete”) gamma function, \(\Gamma(z)\) (available directly as the mpmath function `gamma()`)

  - \(\Gamma(z,a,\infty)\) is the “upper” incomplete gamma function, \(\Gamma(z,a)\)

  - \(\Gamma(z,0,b)\) is the “lower” incomplete gamma function, \(\gamma(z,b)\).

  Of course, we have \(\Gamma(z,0,x) + \Gamma(z,x,\infty) = \Gamma(z)\) for all \(z\) and \(x\).

  Note however that some authors reverse the order of the arguments when defining the lower and upper incomplete gamma function, so one should be careful to get the correct definition.

  If also given the keyword argument `regularized=True`, `gammainc()` computes the “regularized” incomplete gamma function

  \[P(z,a,b) = \frac{\Gamma(z,a,b)}{\Gamma(z)}.\]

  Examples

  We can compare with numerical quadrature to verify that `gammainc()` computes the integral in the definition:

      >>> from mpmath import (mp, gammainc, quad, exp, findroot, mpf, sqrt,
      ...                     erf, pi, identify, ei, lower_gamma, upper_gamma)
      >>> mp.dps = 25
      >>> mp.pretty = True
      >>> gammainc(2+3j, 4, 10)
      (0.00977212668627705160602312 - 0.0770637306312989892451977j)
      >>> quad(lambda t: t**(2+3j-1) * exp(-t), [4, 10])
      (0.00977212668627705160602312 - 0.0770637306312989892451977j)

  Argument symmetries follow directly from the integral definition:

      >>> gammainc(3, 4, 5) + gammainc(3, 5, 4)
      0.0
      >>> lower_gamma(3,2) + gammainc(3,2,4)
      1.523793388892911312363331
      >>> lower_gamma(3,4)
      1.523793388892911312363331
      >>> findroot(lambda z: gammainc(2,z,3), 1)
      3.0

  Evaluation for arbitrarily large arguments:

      >>> upper_gamma(10, 100)
      4.083660630910611272288592e-26
      >>> upper_gamma(10, 10000000000000000)
      5.290402449901174752972486e-4342944819032375
      >>> upper_gamma(3+4j, 1000000+1000000j)
      (-1.257913707524362408877881e-434284 + 2.556691003883483531962095e-434284j)

  Evaluation of a generalized incomplete gamma function automatically chooses the representation that gives a more accurate result, depending on which parameter is larger:

      >>> upper_gamma(10000000, 3) - upper_gamma(10000000, 2)   # Bad
      0.0
      >>> gammainc(10000000, 2, 3)   # Good
      1.755146243738946045873491e+4771204
      >>> lower_gamma(2, 100000001) - lower_gamma(2, 100000000)   # Bad
      0.0
      >>> gammainc(2, 100000000, 100000001)   # Good
      4.078258353474186729184421e-43429441

  The incomplete gamma functions satisfy simple recurrence relations:

      >>> mp.dps = 25
      >>> z, a = mpf(3.5), mpf(2)
      >>> upper_gamma(z+1, a)
      10.60130296933533459267329
      >>> z*upper_gamma(z,a) + a**z*exp(-a)
      10.60130296933533459267329
      >>> lower_gamma(z+1,a)
      1.030425427232114336470932
      >>> z*lower_gamma(z,a) - a**z*exp(-a)
      1.030425427232114336470932

  Evaluation at integers and poles:

      >>> gammainc(-3, -4, -5)
      (-0.2214577048967798566234192 + 0.0j)
      >>> lower_gamma(-3, 5)
      inf

  If \(z\) is an integer, the recurrence reduces the incomplete gamma function to \(P(a) \exp(-a) + Q(b) \exp(-b)\) where \(P\) and \(Q\) are polynomials:

      >>> upper_gamma(1, 2)
      0.1353352832366126918939995
      >>> exp(-2)
      0.1353352832366126918939995
      >>> mp.dps = 50
      >>> identify(gammainc(6, 1, 2), ['exp(-1)', 'exp(-2)'])
      '(326*exp(-1) + (-872)*exp(-2))'

  The incomplete gamma functions reduce to functions such as the exponential integral Ei and the error function for special arguments:

      >>> mp.dps = 25
      >>> upper_gamma(0, 4)
      0.00377935240984890647887486
      >>> -ei(-4)
      0.00377935240984890647887486
      >>> lower_gamma(0.5, 2)
      1.691806732945198336509541
      >>> sqrt(pi)*erf(sqrt(2))
      1.691806732945198336509541

  Related functions

  See also `lower_gamma()` and `upper_gamma()`.

mpmath.lower_gamma(z, b, regularized=False)

  `lower_gamma(z, b)` is the “lower” incomplete gamma function.

  \[\Gamma(z,0,b) = \int_0^b t^{z-1} e^{-t} \, dt\]

  See also `gammainc()`.

mpmath.upper_gamma(z, a, regularized=False)

  `upper_gamma(z, a)` is the “upper” incomplete gamma function.

  \[\Gamma(z,a,\infty) = \int_a^{\infty} t^{z-1} e^{-t} \, dt\]

  See also `gammainc()`.

## Exponential integrals

mpmath.ei(x, **kwargs)

  Computes the exponential integral or Ei-function, \(\mathrm{Ei}(x)\). The exponential integral is defined as

  \[\mathrm{Ei}(x) = \int_{-\infty\,}^x \frac{e^t}{t} \, dt.\]

  When the integration range includes \(t = 0\), the exponential integral is interpreted as providing the Cauchy principal value.

  For real \(x\), the Ei-function behaves roughly like \(\mathrm{Ei}(x) \approx \exp(x) + \log(|x|)\).

  The Ei-function is related to the more general family of exponential integral functions denoted by \(E_n\), which are available as `expint()`.

  Basic examples

  Some basic values and limits are:

      >>> from mpmath import (mp, ei, inf, quad, exp, chop, si, pi, j, chi,
      ...                     shi, hyper, euler, ln, ci)
      >>> mp.pretty = True
      >>> ei(0)
      -inf
      >>> ei(1)
      1.89511781635594
      >>> ei(inf)
      inf
      >>> ei(-inf)
      0.0

  For \(x < 0\), the defining integral can be evaluated numerically as a reference:

      >>> ei(-4)
      -0.00377935240984891
      >>> quad(lambda t: exp(t)/t, [-inf, -4])
      -0.00377935240984891

  `ei()` supports complex arguments and arbitrary precision evaluation:

      >>> mp.dps = 50
      >>> ei(pi)
      10.928374389331410348638445906907535171566338835056
      >>> mp.dps = 25
      >>> ei(3+4j)
      (-4.154091651642689822535359 + 4.294418620024357476985535j)

  Related functions

  The exponential integral is closely related to the logarithmic integral. See `li()` for additional information.

  The exponential integral is related to the hyperbolic and trigonometric integrals (see `chi()`, `shi()`, `ci()`, `si()`) similarly to how the ordinary exponential function is related to the hyperbolic and trigonometric functions:

      >>> mp.dps = 15
      >>> ei(3)
      9.93383257062542
      >>> chi(3) + shi(3)
      9.93383257062542
      >>> chop(ci(3j) - j*si(3j) - pi*j/2)
      9.93383257062542

  Beware that logarithmic corrections, as in the last example above, are required to obtain the correct branch in general. For details, see [1].

  The exponential integral is also a special case of the hypergeometric function \(\,_2F_2\):

      >>> z = 0.6
      >>> z*hyper([1,1],[2,2],z) + (ln(z)-ln(1/z))/2 + euler
      0.769881289937359
      >>> ei(z)
      0.769881289937359

  References

  1. [WolframFunctions] http://functions.wolfram.com/GammaBetaErf/ExpIntegralEi/27/01/

  2. [AbramowitzStegun], section 5:

  3. [Weisstein] http://mathworld.wolfram.com/En-Function.html

mpmath.e1(x, **kwargs)

  Computes the exponential integral \(\mathrm{E}_1(z)\), given by

  \[\mathrm{E}_1(z) = \int_z^{\infty} \frac{e^{-t}}{t} dt.\]

  This is equivalent to `expint()` with \(n = 1\).

  Examples

  Two ways to evaluate this function:

      >>> from mpmath import mp, e1, expint, ei
      >>> mp.dps = 25
      >>> mp.pretty = True
      >>> e1(6.25)
      0.0002704758872637179088496194
      >>> expint(1,6.25)
      0.0002704758872637179088496194

  The E1-function is essentially the same as the Ei-function (`ei()`) with negated argument, except for an imaginary branch cut term:

      >>> e1(2.5)
      0.02491491787026973549562801
      >>> -ei(-2.5)
      0.02491491787026973549562801
      >>> e1(-2.5)
      (-7.073765894578600711923552 - 3.141592653589793238462643j)
      >>> -ei(2.5)
      -7.073765894578600711923552

mpmath.expint(n, z)

  `expint()` gives the generalized exponential integral or En-function,

  \[\mathrm{E}_n(z) = \int_1^{\infty} \frac{e^{-zt}}{t^n} dt,\]

  where \(n\) and \(z\) may both be complex numbers. The case with \(n = 1\) is also given by `e1()`.

  Examples

  Evaluation at real and complex arguments:

      >>> from mpmath import mp, expint, fac, exp, pi
      >>> mp.dps = 25
      >>> mp.pretty = True
      >>> expint(1, 6.25)
      0.0002704758872637179088496194
      >>> expint(-3, 2+3j)
      (0.00299658467335472929656159 + 0.06100816202125885450319632j)
      >>> expint(2+3j, 4-5j)
      (0.001803529474663565056945248 - 0.002235061547756185403349091j)

  At negative integer values of \(n\), \(E_n(z)\) reduces to a rational-exponential function:

      >>> f = lambda n, z: fac(n)*sum(z**k/fac(k-1) for k in range(1,n+2))/\
      ...     exp(z)/z**(n+2)
      >>> n = 3
      >>> z = 1/pi
      >>> expint(-n,z)
      584.2604820613019908668219
      >>> f(n,z)
      584.2604820613019908668219
      >>> n = 5
      >>> expint(-n,z)
      115366.5762594725451811138
      >>> f(n,z)
      115366.5762594725451811138

## Logarithmic integral

mpmath.li(z, offset=False)

  Computes the logarithmic integral or li-function \(\mathrm{li}(x)\), defined by

  \[\mathrm{li}(x) = \int_0^x \frac{1}{\log t} \, dt\]

  The logarithmic integral has a singularity at \(x = 1\).

  Alternatively, `li(x, offset=True)` computes the offset logarithmic integral (used in number theory)

  \[\mathrm{Li}(x) = \int_2^x \frac{1}{\log t} \, dt.\]

  These two functions are related via the simple identity \(\mathrm{Li}(x) = \mathrm{li}(x) - \mathrm{li}(2)\).

  The logarithmic integral should also not be confused with the polylogarithm (also denoted by Li), which is implemented as `polylog()`.

  Examples

  Some basic values and limits:

      >>> from mpmath import mp, li, findroot, ei, log, quad, inf, ln
      >>> mp.dps = 30
      >>> mp.pretty = True
      >>> li(0)
      0.0
      >>> li(1)
      -inf
      >>> li(1)
      -inf
      >>> li(2)
      1.04516378011749278484458888919
      >>> findroot(li, 2)
      1.45136923488338105028396848589
      >>> li(inf)
      inf
      >>> li(2, offset=True)
      0.0
      >>> li(1, offset=True)
      -inf
      >>> li(0, offset=True)
      -1.04516378011749278484458888919
      >>> li(10, offset=True)
      5.12043572466980515267839286347

  The logarithmic integral can be evaluated for arbitrary complex arguments:

      >>> mp.dps = 20
      >>> li(3+4j)
      (3.1343755504645775265 + 2.6769247817778742392j)

  The logarithmic integral is related to the exponential integral:

      >>> ei(log(3))
      2.1635885946671919729
      >>> li(3)
      2.1635885946671919729

  The logarithmic integral grows like \(O(x/\log(x))\):

      >>> mp.dps = 15
      >>> x = 10**100
      >>> x/log(x)
      4.34294481903252e+97
      >>> li(x)
      4.3619719871407e+97

  The prime number theorem states that the number of primes less than \(x\) is asymptotic to \(\mathrm{Li}(x)\) (equivalently \(\mathrm{li}(x)\)). For example, it is known that there are exactly 1,925,320,391,606,803,968,923 prime numbers less than \(10^{23}\) [1]. The logarithmic integral provides a very accurate estimate:

      >>> li(10**23, offset=True)
      1.92532039161405e+21

  A definite integral is:

      >>> quad(li, [0, 1])
      -0.693147180559945
      >>> -ln(2)
      -0.693147180559945

  References

  1. [Weisstein] http://mathworld.wolfram.com/PrimeCountingFunction.html

  2. [Weisstein] http://mathworld.wolfram.com/LogarithmicIntegral.html

## Trigonometric integrals

mpmath.ci(z)

  Computes the cosine integral,

  \[\mathrm{Ci}(x) = -\int_x^{\infty} \frac{\cos t}{t}\,dt = \gamma + \log x + \int_0^x \frac{\cos t - 1}{t}\,dt\]

  Examples

  Some values and limits:

      >>> from mpmath import (mp, ci, pi, inf, chop, sinc, limit, findroot,
      ...                     cos, quadosc, fac, nsum, euler, j)
      >>> mp.dps = 25
      >>> mp.pretty = True
      >>> ci(0)
      -inf
      >>> ci(1)
      0.3374039229009681346626462
      >>> ci(pi)
      0.07366791204642548599010096
      >>> ci(inf)
      0.0
      >>> ci(-inf)
      (0.0 + 3.141592653589793238462643j)
      >>> ci(2+3j)
      (1.408292501520849518759125 - 2.983617742029605093121118j)

  The cosine integral behaves roughly like the sinc function (see `sinc()`) for large real \(x\):

      >>> ci(10**10)
      -4.875060251748226537857298e-11
      >>> sinc(10**10)
      -4.875060250875106915277943e-11
      >>> chop(limit(ci, inf))
      0.0

  It has infinitely many roots on the positive real axis:

      >>> findroot(ci, 1)
      0.6165054856207162337971104
      >>> findroot(ci, 2)
      3.384180422551186426397851

  Evaluation is supported for \(z\) anywhere in the complex plane:

      >>> ci(10**6*(1+j))
      (4.449410587611035724984376e+434287 + 9.75744874290013526417059e+434287j)

  We can evaluate the defining integral as a reference:

      >>> mp.dps = 15
      >>> -quadosc(lambda t: cos(t)/t, [5, inf], omega=1)
      -0.190029749656644
      >>> ci(5)
      -0.190029749656644

  Some infinite series can be evaluated using the cosine integral:

      >>> nsum(lambda k: (-1)**k/(fac(2*k)*(2*k)), [1,inf])
      -0.239811742000565
      >>> ci(1) - euler
      -0.239811742000565

mpmath.si(z)

  Computes the sine integral,

  \[\mathrm{Si}(x) = \int_0^x \frac{\sin t}{t}\,dt.\]

  The sine integral is thus the antiderivative of the sinc function (see `sinc()`).

  Examples

  Some values and limits:

      >>> from mpmath import si, mp, pi, inf, j, quad, sinc, nsum, fac
      >>> mp.dps = 25
      >>> mp.pretty = True
      >>> si(0)
      0.0
      >>> si(1)
      0.9460830703671830149413533
      >>> si(-1)
      -0.9460830703671830149413533
      >>> si(pi)
      1.851937051982466170361053
      >>> si(inf)
      1.570796326794896619231322
      >>> si(-inf)
      -1.570796326794896619231322
      >>> si(2+3j)
      (4.547513889562289219853204 + 1.399196580646054789459839j)

  The sine integral approaches \(\pi/2\) for large real \(x\):

      >>> si(10**10)
      1.570796326707584656968511
      >>> pi/2
      1.570796326794896619231322

  Evaluation is supported for \(z\) anywhere in the complex plane:

      >>> si(10**6*(1+j))
      (-9.75744874290013526417059e+434287 + 4.449410587611035724984376e+434287j)

  We can evaluate the defining integral as a reference:

      >>> mp.dps = 15
      >>> quad(sinc, [0, 5])
      1.54993124494467
      >>> si(5)
      1.54993124494467

  Some infinite series can be evaluated using the sine integral:

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

## Hyperbolic integrals

mpmath.chi(z)

  Computes the hyperbolic cosine integral, defined in analogy with the cosine integral (see `ci()`) as

  \[\mathrm{Chi}(x) = -\int_x^{\infty} \frac{\cosh t}{t}\,dt = \gamma + \log x + \int_0^x \frac{\cosh t - 1}{t}\,dt\]

  Some values and limits:

      >>> from mpmath import mp, chi, inf, findroot, j
      >>> mp.dps = 25
      >>> mp.pretty = True
      >>> chi(0)
      -inf
      >>> chi(1)
      0.8378669409802082408946786
      >>> chi(inf)
      inf
      >>> findroot(chi, 0.5)
      0.5238225713898644064509583
      >>> chi(2+3j)
      (-0.1683628683277204662429321 + 2.625115880451325002151688j)

  Evaluation is supported for \(z\) anywhere in the complex plane:

      >>> chi(10**6*(1+j))
      (4.449410587611035724984376e+434287 - 9.75744874290013526417059e+434287j)

mpmath.shi(z)

  Computes the hyperbolic sine integral, defined in analogy with the sine integral (see `si()`) as

  \[\mathrm{Shi}(x) = \int_0^x \frac{\sinh t}{t}\,dt.\]

  Some values and limits:

      >>> from mpmath import mp, shi, inf, j
      >>> mp.dps = 25
      >>> mp.pretty = True
      >>> shi(0)
      0.0
      >>> shi(1)
      1.057250875375728514571842
      >>> shi(-1)
      -1.057250875375728514571842
      >>> shi(inf)
      inf
      >>> shi(2+3j)
      (-0.1931890762719198291678095 + 2.645432555362369624818525j)

  Evaluation is supported for \(z\) anywhere in the complex plane:

      >>> shi(10**6*(1+j))
      (4.449410587611035724984376e+434287 - 9.75744874290013526417059e+434287j)

## Error functions

mpmath.erf(z)

  Computes the error function, \(\mathrm{erf}(x)\). The error function is the normalized antiderivative of the Gaussian function \(\exp(-t^2)\). More precisely,

  \[\mathrm{erf}(x) = \frac{2}{\sqrt \pi} \int_0^x \exp(-t^2) \,dt\]

  Basic examples

  Simple values and limits include:

      >>> from mpmath import mp, erf, inf, nprint, chop, taylor, j
      >>> mp.pretty = True
      >>> erf(0)
      0.0
      >>> erf(1)
      0.842700792949715
      >>> erf(-1)
      -0.842700792949715
      >>> erf(inf)
      1.0
      >>> erf(-inf)
      -1.0

  For large real \(x\), \(\mathrm{erf}(x)\) approaches 1 very rapidly:

      >>> erf(3)
      0.999977909503001
      >>> erf(5)
      0.999999999998463

  The error function is an odd function:

      >>> nprint(chop(taylor(erf, 0, 5)))
      [0.0, 1.12838, 0.0, -0.376126, 0.0, 0.112838]

  `erf()` implements arbitrary-precision evaluation and supports complex numbers:

      >>> mp.dps = 50
      >>> erf(0.5)
      0.52049987781304653768274665389196452873645157575796
      >>> mp.dps = 25
      >>> erf(1+j)
      (1.316151281697947644880271 + 0.1904534692378346862841089j)

  Evaluation is supported for large arguments:

      >>> mp.dps = 25
      >>> erf('1e1000')
      1.0
      >>> erf('-1e1000')
      -1.0
      >>> erf('1e-1000')
      1.128379167095512573896159e-1000
      >>> erf('1e7j')
      (0.0 + 8.593897639029319267398803e+43429448190317j)
      >>> erf('1e7+1e7j')
      (0.9999999858172446172631323 + 3.728805278735270407053139e-8j)

  Related functions

  See also `erfc()`, which is more accurate for large \(x\), and `erfi()` which gives the antiderivative of \(\exp(t^2)\).

  The Fresnel integrals `fresnels()` and `fresnelc()` are also related to the error function.

mpmath.erfc(z)

  Computes the complementary error function, \(\mathrm{erfc}(x) = 1-\mathrm{erf}(x)\). This function avoids cancellation that occurs when naively computing the complementary error function as `1-erf(x)`:

      >>> from mpmath import mp, erf, erfc
      >>> mp.pretty = True
      >>> 1 - erf(10)
      0.0
      >>> erfc(10)
      2.08848758376254e-45

  `erfc()` works accurately even for ludicrously large arguments:

      >>> erfc(10**10)
      4.3504398860243e-43429448190325182776

  Complex arguments are supported:

      >>> erfc(500+50j)
      (1.19739830969552e-107492 + 1.46072418957528e-107491j)

mpmath.erfi(z)

  Computes the imaginary error function, \(\mathrm{erfi}(x)\). The imaginary error function is defined in analogy with the error function, but with a positive sign in the integrand:

  \[\mathrm{erfi}(x) = \frac{2}{\sqrt \pi} \int_0^x \exp(t^2) \,dt\]

  Whereas the error function rapidly converges to 1 as \(x\) grows, the imaginary error function rapidly diverges to infinity. The functions are related as \(\mathrm{erfi}(x) = -i\,\mathrm{erf}(ix)\) for all complex numbers \(x\).

  Examples

  Basic values and limits:

      >>> from mpmath import mp, erfi, inf, erf
      >>> mp.pretty = True
      >>> erfi(0)
      0.0
      >>> erfi(1)
      1.65042575879754
      >>> erfi(-1)
      -1.65042575879754
      >>> erfi(inf)
      inf
      >>> erfi(-inf)
      -inf

  Note the symmetry between erf and erfi:

      >>> erfi(3j)
      (0.0 + 0.999977909503001j)
      >>> erf(3)
      0.999977909503001
      >>> erf(1+2j)
      (-0.536643565778565 - 5.04914370344703j)
      >>> erfi(2+1j)
      (-5.04914370344703 - 0.536643565778565j)

  Large arguments are supported:

      >>> erfi(1000)
      1.71130938718796e+434291
      >>> erfi(10**10)
      7.3167287567024e+43429448190325182754
      >>> erfi(-10**10)
      -7.3167287567024e+43429448190325182754
      >>> erfi(1000-500j)
      (2.49895233563961e+325717 + 2.6846779342253e+325717j)
      >>> erfi(100000j)
      (0.0 + 1.0j)
      >>> erfi(-100000j)
      (0.0 - 1.0j)

mpmath.erfinv(x)

  Computes the inverse error function, satisfying

  \[\mathrm{erf}(\mathrm{erfinv}(x)) = \mathrm{erfinv}(\mathrm{erf}(x)) = x.\]

  This function is defined only for \(-1 \le x \le 1\).

  Examples

  Special values include:

      >>> from mpmath import mp, erfinv, erf, quad, sqrt, pi
      >>> mp.pretty = True
      >>> erfinv(0)
      0.0
      >>> erfinv(1)
      inf
      >>> erfinv(-1)
      -inf

  The domain is limited to the standard interval:

      >>> erfinv(2)
      Traceback (most recent call last):
        ...
      ValueError: erfinv(x) is defined only for -1 <= x <= 1

  It is simple to check that `erfinv()` computes inverse values of `erf()` as promised:

      >>> erf(erfinv(0.75))
      0.75
      >>> erf(erfinv(-0.995))
      -0.995

  `erfinv()` supports arbitrary-precision evaluation:

      >>> mp.dps = 50
      >>> x = erf(2)
      >>> x
      0.99532226501895273416206925636725292861089179704006
      >>> erfinv(x)
      2.0

  A definite integral involving the inverse error function:

      >>> mp.dps = 15
      >>> quad(erfinv, [0, 1])
      0.564189583547756
      >>> 1/sqrt(pi)
      0.564189583547756

  The inverse error function can be used to generate random numbers with a Gaussian distribution (although this is a relatively inefficient algorithm):

      >>> nprint([erfinv(2*rand()-1) for n in range(6)])
      [-0.586747, 1.10233, -0.376796, 0.926037, -0.708142, -0.732012]

## The normal distribution

mpmath.npdf(x, mu=0, sigma=1)

  `npdf(x, mu=0, sigma=1)` evaluates the probability density function of a normal distribution with mean value \(\mu\) and variance \(\sigma^2\).

  Elementary properties of the probability distribution can be verified using numerical integration:

      >>> from mpmath import mp, quad, inf, npdf
      >>> mp.pretty = True
      >>> quad(npdf, [-inf, inf])
      1.0
      >>> quad(lambda x: npdf(x, 3), [3, inf])
      0.5
      >>> quad(lambda x: npdf(x, 3, 2), [3, inf])
      0.5

  See also `ncdf()`, which gives the cumulative distribution.

mpmath.ncdf(x, mu=0, sigma=1)

  `ncdf(x, mu=0, sigma=1)` evaluates the cumulative distribution function of a normal distribution with mean value \(\mu\) and variance \(\sigma^2\).

  See also `npdf()`, which gives the probability density.

  Elementary properties include:

      >>> from mpmath import mp, ncdf, pi, inf, diff, npdf
      >>> mp.pretty = True
      >>> ncdf(pi, mu=pi)
      0.5
      >>> ncdf(-inf)
      0.0
      >>> ncdf(+inf)
      1.0

  The cumulative distribution is the integral of the density function having identical mu and sigma:

      >>> mp.dps = 15
      >>> diff(ncdf, 2)
      0.053990966513188
      >>> npdf(2)
      0.053990966513188
      >>> diff(lambda x: ncdf(x, 1, 0.5), 0)
      0.107981933026376
      >>> npdf(0, 1, 0.5)
      0.107981933026376

## Fresnel integrals

mpmath.fresnels(z)

  Computes the Fresnel sine integral

  \[S(x) = \int_0^x \sin\left(\frac{\pi t^2}{2}\right) \,dt\]

  Note that some sources define this function without the normalization factor \(\pi/2\).

  Examples

  Some basic values and limits:

      >>> from mpmath import mp, fresnels, inf, quad, sin, pi
      >>> mp.dps = 25
      >>> mp.pretty = True
      >>> fresnels(0)
      0.0
      >>> fresnels(inf)
      0.5
      >>> fresnels(-inf)
      -0.5
      >>> fresnels(1)
      0.4382591473903547660767567
      >>> fresnels(1+2j)
      (36.72546488399143842838788 + 15.58775110440458732748279j)

  Comparing with the definition:

      >>> fresnels(3)
      0.4963129989673750360976123
      >>> quad(lambda t: sin(pi*t**2/2), [0,3])
      0.4963129989673750360976123

mpmath.fresnelc(z)

  Computes the Fresnel cosine integral

  \[C(x) = \int_0^x \cos\left(\frac{\pi t^2}{2}\right) \,dt\]

  Note that some sources define this function without the normalization factor \(\pi/2\).

  Examples

  Some basic values and limits:

      >>> from mpmath import mp, fresnelc, inf, quad, cos, pi
      >>> mp.dps = 25
      >>> mp.pretty = True
      >>> fresnelc(0)
      0.0
      >>> fresnelc(inf)
      0.5
      >>> fresnelc(-inf)
      -0.5
      >>> fresnelc(1)
      0.7798934003768228294742064
      >>> fresnelc(1+2j)
      (16.08787137412548041729489 - 36.22568799288165021578758j)

  Comparing with the definition:

      >>> fresnelc(3)
      0.6057207892976856295561611
      >>> quad(lambda t: cos(pi*t**2/2), [0,3])
      0.6057207892976856295561611