""" This module contains various functions that are special cases
of incomplete gamma functions. It should probably be renamed. """
from __future__ import print_function, division
from sympy.core import Add, S, sympify, cacheit, pi, I
from sympy.core.function import Function, ArgumentIndexError
from sympy.core.symbol import Symbol
from sympy.functions.combinatorial.factorials import factorial
from sympy.functions.elementary.integers import floor
from sympy.functions.elementary.miscellaneous import sqrt, root
from sympy.functions.elementary.exponential import exp, log
from sympy.functions.elementary.complexes import polar_lift
from sympy.functions.elementary.hyperbolic import cosh, sinh
from sympy.functions.elementary.trigonometric import cos, sin, sinc
from sympy.functions.special.hyper import hyper, meijerg
from sympy.core.compatibility import range
# TODO series expansions
# TODO see the "Note:" in Ei
###############################################################################
################################ ERROR FUNCTION ###############################
###############################################################################
[docs]class erf(Function):
r"""
The Gauss error function. This function is defined as:
.. math ::
\mathrm{erf}(x) = \frac{2}{\sqrt{\pi}} \int_0^x e^{-t^2} \mathrm{d}t.
Examples
========
>>> from sympy import I, oo, erf
>>> from sympy.abc import z
Several special values are known:
>>> erf(0)
0
>>> erf(oo)
1
>>> erf(-oo)
-1
>>> erf(I*oo)
oo*I
>>> erf(-I*oo)
-oo*I
In general one can pull out factors of -1 and I from the argument:
>>> erf(-z)
-erf(z)
The error function obeys the mirror symmetry:
>>> from sympy import conjugate
>>> conjugate(erf(z))
erf(conjugate(z))
Differentiation with respect to z is supported:
>>> from sympy import diff
>>> diff(erf(z), z)
2*exp(-z**2)/sqrt(pi)
We can numerically evaluate the error function to arbitrary precision
on the whole complex plane:
>>> erf(4).evalf(30)
0.999999984582742099719981147840
>>> erf(-4*I).evalf(30)
-1296959.73071763923152794095062*I
See Also
========
erfc: Complementary error function.
erfi: Imaginary error function.
erf2: Two-argument error function.
erfinv: Inverse error function.
erfcinv: Inverse Complementary error function.
erf2inv: Inverse two-argument error function.
References
==========
.. [1] http://en.wikipedia.org/wiki/Error_function
.. [2] http://dlmf.nist.gov/7
.. [3] http://mathworld.wolfram.com/Erf.html
.. [4] http://functions.wolfram.com/GammaBetaErf/Erf
"""
unbranched = True
def fdiff(self, argindex=1):
if argindex == 1:
return 2*exp(-self.args[0]**2)/sqrt(S.Pi)
else:
raise ArgumentIndexError(self, argindex)
def inverse(self, argindex=1):
"""
Returns the inverse of this function.
"""
return erfinv
@classmethod
def eval(cls, arg):
if arg.is_Number:
if arg is S.NaN:
return S.NaN
elif arg is S.Infinity:
return S.One
elif arg is S.NegativeInfinity:
return S.NegativeOne
elif arg is S.Zero:
return S.Zero
if isinstance(arg, erfinv):
return arg.args[0]
if isinstance(arg, erfcinv):
return S.One - arg.args[0]
if isinstance(arg, erf2inv) and arg.args[0] is S.Zero:
return arg.args[1]
# Try to pull out factors of I
t = arg.extract_multiplicatively(S.ImaginaryUnit)
if t is S.Infinity or t is S.NegativeInfinity:
return arg
# Try to pull out factors of -1
if arg.could_extract_minus_sign():
return -cls(-arg)
@staticmethod
@cacheit
def taylor_term(n, x, *previous_terms):
if n < 0 or n % 2 == 0:
return S.Zero
else:
x = sympify(x)
k = floor((n - 1)/S(2))
if len(previous_terms) > 2:
return -previous_terms[-2] * x**2 * (n - 2)/(n*k)
else:
return 2*(-1)**k * x**n/(n*factorial(k)*sqrt(S.Pi))
def _eval_conjugate(self):
return self.func(self.args[0].conjugate())
def _eval_is_real(self):
return self.args[0].is_real
def _eval_rewrite_as_uppergamma(self, z):
from sympy import uppergamma
return sqrt(z**2)/z*(S.One - uppergamma(S.Half, z**2)/sqrt(S.Pi))
def _eval_rewrite_as_fresnels(self, z):
arg = (S.One - S.ImaginaryUnit)*z/sqrt(pi)
return (S.One + S.ImaginaryUnit)*(fresnelc(arg) - I*fresnels(arg))
def _eval_rewrite_as_fresnelc(self, z):
arg = (S.One - S.ImaginaryUnit)*z/sqrt(pi)
return (S.One + S.ImaginaryUnit)*(fresnelc(arg) - I*fresnels(arg))
def _eval_rewrite_as_meijerg(self, z):
return z/sqrt(pi)*meijerg([S.Half], [], [0], [-S.Half], z**2)
def _eval_rewrite_as_hyper(self, z):
return 2*z/sqrt(pi)*hyper([S.Half], [3*S.Half], -z**2)
def _eval_rewrite_as_expint(self, z):
return sqrt(z**2)/z - z*expint(S.Half, z**2)/sqrt(S.Pi)
def _eval_rewrite_as_tractable(self, z):
return S.One - _erfs(z)*exp(-z**2)
def _eval_rewrite_as_erfc(self, z):
return S.One - erfc(z)
def _eval_rewrite_as_erfi(self, z):
return -I*erfi(I*z)
def _eval_as_leading_term(self, x):
from sympy import Order
arg = self.args[0].as_leading_term(x)
if x in arg.free_symbols and Order(1, x).contains(arg):
return 2*x/sqrt(pi)
else:
return self.func(arg)
def as_real_imag(self, deep=True, **hints):
if self.args[0].is_real:
if deep:
hints['complex'] = False
return (self.expand(deep, **hints), S.Zero)
else:
return (self, S.Zero)
if deep:
x, y = self.args[0].expand(deep, **hints).as_real_imag()
else:
x, y = self.args[0].as_real_imag()
sq = -y**2/x**2
re = S.Half*(self.func(x + x*sqrt(sq)) + self.func(x - x*sqrt(sq)))
im = x/(2*y) * sqrt(sq) * (self.func(x - x*sqrt(sq)) -
self.func(x + x*sqrt(sq)))
return (re, im)
[docs]class erfc(Function):
r"""
Complementary Error Function. The function is defined as:
.. math ::
\mathrm{erfc}(x) = \frac{2}{\sqrt{\pi}} \int_x^\infty e^{-t^2} \mathrm{d}t
Examples
========
>>> from sympy import I, oo, erfc
>>> from sympy.abc import z
Several special values are known:
>>> erfc(0)
1
>>> erfc(oo)
0
>>> erfc(-oo)
2
>>> erfc(I*oo)
-oo*I
>>> erfc(-I*oo)
oo*I
The error function obeys the mirror symmetry:
>>> from sympy import conjugate
>>> conjugate(erfc(z))
erfc(conjugate(z))
Differentiation with respect to z is supported:
>>> from sympy import diff
>>> diff(erfc(z), z)
-2*exp(-z**2)/sqrt(pi)
It also follows
>>> erfc(-z)
-erfc(z) + 2
We can numerically evaluate the complementary error function to arbitrary precision
on the whole complex plane:
>>> erfc(4).evalf(30)
0.0000000154172579002800188521596734869
>>> erfc(4*I).evalf(30)
1.0 - 1296959.73071763923152794095062*I
See Also
========
erf: Gaussian error function.
erfi: Imaginary error function.
erf2: Two-argument error function.
erfinv: Inverse error function.
erfcinv: Inverse Complementary error function.
erf2inv: Inverse two-argument error function.
References
==========
.. [1] http://en.wikipedia.org/wiki/Error_function
.. [2] http://dlmf.nist.gov/7
.. [3] http://mathworld.wolfram.com/Erfc.html
.. [4] http://functions.wolfram.com/GammaBetaErf/Erfc
"""
unbranched = True
def fdiff(self, argindex=1):
if argindex == 1:
return -2*exp(-self.args[0]**2)/sqrt(S.Pi)
else:
raise ArgumentIndexError(self, argindex)
def inverse(self, argindex=1):
"""
Returns the inverse of this function.
"""
return erfcinv
@classmethod
def eval(cls, arg):
if arg.is_Number:
if arg is S.NaN:
return S.NaN
elif arg is S.Infinity:
return S.Zero
elif arg is S.Zero:
return S.One
if isinstance(arg, erfinv):
return S.One - arg.args[0]
if isinstance(arg, erfcinv):
return arg.args[0]
# Try to pull out factors of I
t = arg.extract_multiplicatively(S.ImaginaryUnit)
if t is S.Infinity or t is S.NegativeInfinity:
return -arg
# Try to pull out factors of -1
if arg.could_extract_minus_sign():
return S(2) - cls(-arg)
@staticmethod
@cacheit
def taylor_term(n, x, *previous_terms):
if n == 0:
return S.One
elif n < 0 or n % 2 == 0:
return S.Zero
else:
x = sympify(x)
k = floor((n - 1)/S(2))
if len(previous_terms) > 2:
return -previous_terms[-2] * x**2 * (n - 2)/(n*k)
else:
return -2*(-1)**k * x**n/(n*factorial(k)*sqrt(S.Pi))
def _eval_conjugate(self):
return self.func(self.args[0].conjugate())
def _eval_is_real(self):
return self.args[0].is_real
def _eval_rewrite_as_tractable(self, z):
return self.rewrite(erf).rewrite("tractable", deep=True)
def _eval_rewrite_as_erf(self, z):
return S.One - erf(z)
def _eval_rewrite_as_erfi(self, z):
return S.One + I*erfi(I*z)
def _eval_rewrite_as_fresnels(self, z):
arg = (S.One - S.ImaginaryUnit)*z/sqrt(pi)
return S.One - (S.One + S.ImaginaryUnit)*(fresnelc(arg) - I*fresnels(arg))
def _eval_rewrite_as_fresnelc(self, z):
arg = (S.One-S.ImaginaryUnit)*z/sqrt(pi)
return S.One - (S.One + S.ImaginaryUnit)*(fresnelc(arg) - I*fresnels(arg))
def _eval_rewrite_as_meijerg(self, z):
return S.One - z/sqrt(pi)*meijerg([S.Half], [], [0], [-S.Half], z**2)
def _eval_rewrite_as_hyper(self, z):
return S.One - 2*z/sqrt(pi)*hyper([S.Half], [3*S.Half], -z**2)
def _eval_rewrite_as_uppergamma(self, z):
from sympy import uppergamma
return S.One - sqrt(z**2)/z*(S.One - uppergamma(S.Half, z**2)/sqrt(S.Pi))
def _eval_rewrite_as_expint(self, z):
return S.One - sqrt(z**2)/z + z*expint(S.Half, z**2)/sqrt(S.Pi)
def _eval_expand_func(self, **hints):
return self.rewrite(erf)
def _eval_as_leading_term(self, x):
from sympy import Order
arg = self.args[0].as_leading_term(x)
if x in arg.free_symbols and Order(1, x).contains(arg):
return S.One
else:
return self.func(arg)
def as_real_imag(self, deep=True, **hints):
if self.args[0].is_real:
if deep:
hints['complex'] = False
return (self.expand(deep, **hints), S.Zero)
else:
return (self, S.Zero)
if deep:
x, y = self.args[0].expand(deep, **hints).as_real_imag()
else:
x, y = self.args[0].as_real_imag()
sq = -y**2/x**2
re = S.Half*(self.func(x + x*sqrt(sq)) + self.func(x - x*sqrt(sq)))
im = x/(2*y) * sqrt(sq) * (self.func(x - x*sqrt(sq)) -
self.func(x + x*sqrt(sq)))
return (re, im)
[docs]class erfi(Function):
r"""
Imaginary error function. The function erfi is defined as:
.. math ::
\mathrm{erfi}(x) = \frac{2}{\sqrt{\pi}} \int_0^x e^{t^2} \mathrm{d}t
Examples
========
>>> from sympy import I, oo, erfi
>>> from sympy.abc import z
Several special values are known:
>>> erfi(0)
0
>>> erfi(oo)
oo
>>> erfi(-oo)
-oo
>>> erfi(I*oo)
I
>>> erfi(-I*oo)
-I
In general one can pull out factors of -1 and I from the argument:
>>> erfi(-z)
-erfi(z)
>>> from sympy import conjugate
>>> conjugate(erfi(z))
erfi(conjugate(z))
Differentiation with respect to z is supported:
>>> from sympy import diff
>>> diff(erfi(z), z)
2*exp(z**2)/sqrt(pi)
We can numerically evaluate the imaginary error function to arbitrary precision
on the whole complex plane:
>>> erfi(2).evalf(30)
18.5648024145755525987042919132
>>> erfi(-2*I).evalf(30)
-0.995322265018952734162069256367*I
See Also
========
erf: Gaussian error function.
erfc: Complementary error function.
erf2: Two-argument error function.
erfinv: Inverse error function.
erfcinv: Inverse Complementary error function.
erf2inv: Inverse two-argument error function.
References
==========
.. [1] http://en.wikipedia.org/wiki/Error_function
.. [2] http://mathworld.wolfram.com/Erfi.html
.. [3] http://functions.wolfram.com/GammaBetaErf/Erfi
"""
unbranched = True
def fdiff(self, argindex=1):
if argindex == 1:
return 2*exp(self.args[0]**2)/sqrt(S.Pi)
else:
raise ArgumentIndexError(self, argindex)
@classmethod
def eval(cls, z):
if z.is_Number:
if z is S.NaN:
return S.NaN
elif z is S.Zero:
return S.Zero
elif z is S.Infinity:
return S.Infinity
# Try to pull out factors of -1
if z.could_extract_minus_sign():
return -cls(-z)
# Try to pull out factors of I
nz = z.extract_multiplicatively(I)
if nz is not None:
if nz is S.Infinity:
return I
if isinstance(nz, erfinv):
return I*nz.args[0]
if isinstance(nz, erfcinv):
return I*(S.One - nz.args[0])
if isinstance(nz, erf2inv) and nz.args[0] is S.Zero:
return I*nz.args[1]
@staticmethod
@cacheit
def taylor_term(n, x, *previous_terms):
if n < 0 or n % 2 == 0:
return S.Zero
else:
x = sympify(x)
k = floor((n - 1)/S(2))
if len(previous_terms) > 2:
return previous_terms[-2] * x**2 * (n - 2)/(n*k)
else:
return 2 * x**n/(n*factorial(k)*sqrt(S.Pi))
def _eval_conjugate(self):
return self.func(self.args[0].conjugate())
def _eval_is_real(self):
return self.args[0].is_real
def _eval_rewrite_as_tractable(self, z):
return self.rewrite(erf).rewrite("tractable", deep=True)
def _eval_rewrite_as_erf(self, z):
return -I*erf(I*z)
def _eval_rewrite_as_erfc(self, z):
return I*erfc(I*z) - I
def _eval_rewrite_as_fresnels(self, z):
arg = (S.One + S.ImaginaryUnit)*z/sqrt(pi)
return (S.One - S.ImaginaryUnit)*(fresnelc(arg) - I*fresnels(arg))
def _eval_rewrite_as_fresnelc(self, z):
arg = (S.One + S.ImaginaryUnit)*z/sqrt(pi)
return (S.One - S.ImaginaryUnit)*(fresnelc(arg) - I*fresnels(arg))
def _eval_rewrite_as_meijerg(self, z):
return z/sqrt(pi)*meijerg([S.Half], [], [0], [-S.Half], -z**2)
def _eval_rewrite_as_hyper(self, z):
return 2*z/sqrt(pi)*hyper([S.Half], [3*S.Half], z**2)
def _eval_rewrite_as_uppergamma(self, z):
from sympy import uppergamma
return sqrt(-z**2)/z*(uppergamma(S.Half, -z**2)/sqrt(S.Pi) - S.One)
def _eval_rewrite_as_expint(self, z):
return sqrt(-z**2)/z - z*expint(S.Half, -z**2)/sqrt(S.Pi)
def _eval_expand_func(self, **hints):
return self.rewrite(erf)
def as_real_imag(self, deep=True, **hints):
if self.args[0].is_real:
if deep:
hints['complex'] = False
return (self.expand(deep, **hints), S.Zero)
else:
return (self, S.Zero)
if deep:
x, y = self.args[0].expand(deep, **hints).as_real_imag()
else:
x, y = self.args[0].as_real_imag()
sq = -y**2/x**2
re = S.Half*(self.func(x + x*sqrt(sq)) + self.func(x - x*sqrt(sq)))
im = x/(2*y) * sqrt(sq) * (self.func(x - x*sqrt(sq)) -
self.func(x + x*sqrt(sq)))
return (re, im)
[docs]class erf2(Function):
r"""
Two-argument error function. This function is defined as:
.. math ::
\mathrm{erf2}(x, y) = \frac{2}{\sqrt{\pi}} \int_x^y e^{-t^2} \mathrm{d}t
Examples
========
>>> from sympy import I, oo, erf2
>>> from sympy.abc import x, y
Several special values are known:
>>> erf2(0, 0)
0
>>> erf2(x, x)
0
>>> erf2(x, oo)
-erf(x) + 1
>>> erf2(x, -oo)
-erf(x) - 1
>>> erf2(oo, y)
erf(y) - 1
>>> erf2(-oo, y)
erf(y) + 1
In general one can pull out factors of -1:
>>> erf2(-x, -y)
-erf2(x, y)
The error function obeys the mirror symmetry:
>>> from sympy import conjugate
>>> conjugate(erf2(x, y))
erf2(conjugate(x), conjugate(y))
Differentiation with respect to x, y is supported:
>>> from sympy import diff
>>> diff(erf2(x, y), x)
-2*exp(-x**2)/sqrt(pi)
>>> diff(erf2(x, y), y)
2*exp(-y**2)/sqrt(pi)
See Also
========
erf: Gaussian error function.
erfc: Complementary error function.
erfi: Imaginary error function.
erfinv: Inverse error function.
erfcinv: Inverse Complementary error function.
erf2inv: Inverse two-argument error function.
References
==========
.. [1] http://functions.wolfram.com/GammaBetaErf/Erf2/
"""
def fdiff(self, argindex):
x, y = self.args
if argindex == 1:
return -2*exp(-x**2)/sqrt(S.Pi)
elif argindex == 2:
return 2*exp(-y**2)/sqrt(S.Pi)
else:
raise ArgumentIndexError(self, argindex)
@classmethod
def eval(cls, x, y):
I = S.Infinity
N = S.NegativeInfinity
O = S.Zero
if x is S.NaN or y is S.NaN:
return S.NaN
elif x == y:
return S.Zero
elif (x is I or x is N or x is O) or (y is I or y is N or y is O):
return erf(y) - erf(x)
if isinstance(y, erf2inv) and y.args[0] == x:
return y.args[1]
#Try to pull out -1 factor
sign_x = x.could_extract_minus_sign()
sign_y = y.could_extract_minus_sign()
if (sign_x and sign_y):
return -cls(-x, -y)
elif (sign_x or sign_y):
return erf(y)-erf(x)
def _eval_conjugate(self):
return self.func(self.args[0].conjugate(), self.args[1].conjugate())
def _eval_is_real(self):
return self.args[0].is_real and self.args[1].is_real
def _eval_rewrite_as_erf(self, x, y):
return erf(y) - erf(x)
def _eval_rewrite_as_erfc(self, x, y):
return erfc(x) - erfc(y)
def _eval_rewrite_as_erfi(self, x, y):
return I*(erfi(I*x)-erfi(I*y))
def _eval_rewrite_as_fresnels(self, x, y):
return erf(y).rewrite(fresnels) - erf(x).rewrite(fresnels)
def _eval_rewrite_as