Coverage for sympy/functions/special/error

Hot-keys on this page

r m x p toggle line displays

j k next/prev highlighted chunk

0 (zero) top of page

1 (one) first highlighted chunk

""" 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 ###############################

###############################################################################

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)

>>> erf(oo)

>>> 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 arg.func is erfinv:

return arg.args[0]

if arg.func is erfcinv:

return S.One - arg.args[0]

if arg.func is 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)

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)

>>> erfc(oo)

>>> erfc(-oo)

>>> 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 arg.func is erfinv:

return S.One - arg.args[0]

if arg.func is 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_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)

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)

>>> erfi(oo)

>>> erfi(-oo)

-oo

>>> erfi(I*oo)

>>> 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 nz.func is erfinv:

return I*nz.args[0]

if nz.func is erfcinv:

return I*(S.One - nz.args[0])

if nz.func is 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 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)

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)

>>> erf2(x, x)

>>> 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 y.func is 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_fresnelc(self, x, y):

return erf(y).rewrite(fresnelc) - erf(x).rewrite(fresnelc)

def _eval_rewrite_as_meijerg(self, x, y):

return erf(y).rewrite(meijerg) - erf(x).rewrite(meijerg)

def _eval_rewrite_as_hyper(self, x, y):

return erf(y).rewrite(hyper) - erf(x).rewrite(hyper)

def _eval_rewrite_as_uppergamma(self, x, y):

from sympy import uppergamma

return (sqrt(y**2)/y*(S.One - uppergamma(S.Half, y**2)/sqrt(S.Pi)) -

sqrt(x**2)/x*(S.One - uppergamma(S.Half, x**2)/sqrt(S.Pi)))

def _eval_rewrite_as_expint(self, x, y):

return erf(y).rewrite(expint) - erf(x).rewrite(expint)

class erfinv(Function):

r"""

Inverse Error Function. The erfinv function is defined as:

.. math ::

\mathrm{erf}(x) = y \quad \Rightarrow \quad \mathrm{erfinv}(y) = x

Examples

========

>>> from sympy import I, oo, erfinv

>>> from sympy.abc import x

Several special values are known:

>>> erfinv(0)

>>> erfinv(1)

Differentiation with respect to x is supported:

>>> from sympy import diff

>>> diff(erfinv(x), x)

sqrt(pi)*exp(erfinv(x)**2)/2

We can numerically evaluate the inverse error function to arbitrary precision

on [-1, 1]:

>>> erfinv(0.2).evalf(30)

0.179143454621291692285822705344

See Also

========

erf: Gaussian error function.

erfc: Complementary error function.

erfi: Imaginary error function.

erf2: Two-argument error function.

erfcinv: Inverse Complementary error function.

erf2inv: Inverse two-argument error function.

References

==========

.. [1] http://en.wikipedia.org/wiki/Error_function#Inverse_functions

.. [2] http://functions.wolfram.com/GammaBetaErf/InverseErf/

"""

def fdiff(self, argindex =1):

if argindex == 1:

return sqrt(S.Pi)*exp(self.func(self.args[0])**2)*S.Half

else :

raise ArgumentIndexError(self, argindex)

def inverse(self, argindex=1):

"""

Returns the inverse of this function.

"""

return erf

@classmethod

def eval(cls, z):

if z is S.NaN:

return S.NaN

elif z is S.NegativeOne:

return S.NegativeInfinity

elif z is S.Zero:

return S.Zero

elif z is S.One:

return S.Infinity

if (z.func is erf) and z.args[0].is_real:

return z.args[0]

# Try to pull out factors of -1

nz = z.extract_multiplicatively(-1)

if nz is not None and ((nz.func is erf) and (nz.args[0]).is_real):

return -nz.args[0]

def _eval_rewrite_as_erfcinv(self, z):

return erfcinv(1-z)

class erfcinv (Function):

r"""

Inverse Complementary Error Function. The erfcinv function is defined as:

.. math ::

\mathrm{erfc}(x) = y \quad \Rightarrow \quad \mathrm{erfcinv}(y) = x

Examples

========

>>> from sympy import I, oo, erfcinv

>>> from sympy.abc import x

Several special values are known:

>>> erfcinv(1)

>>> erfcinv(0)

Differentiation with respect to x is supported:

>>> from sympy import diff

>>> diff(erfcinv(x), x)

-sqrt(pi)*exp(erfcinv(x)**2)/2

See Also

========

erf: Gaussian error function.

erfc: Complementary error function.

erfi: Imaginary error function.

erf2: Two-argument error function.

erfinv: Inverse error function.

erf2inv: Inverse two-argument error function.

References

==========

.. [1] http://en.wikipedia.org/wiki/Error_function#Inverse_functions

.. [2] http://functions.wolfram.com/GammaBetaErf/InverseErfc/

"""

def fdiff(self, argindex =1):

if argindex == 1:

return -sqrt(S.Pi)*exp(self.func(self.args[0])**2)*S.Half

else:

raise ArgumentIndexError(self, argindex)

def inverse(self, argindex=1):

"""

Returns the inverse of this function.

"""

return erfc

@classmethod

def eval(cls, z):

if z is S.NaN:

return S.NaN

elif z is S.Zero:

return S.Infinity

elif z is S.One:

return S.Zero

elif z == 2:

return S.NegativeInfinity

def _eval_rewrite_as_erfinv(self, z):

return erfinv(1-z)

class erf2inv(Function):

r"""

Two-argument Inverse error function. The erf2inv function is defined as:

.. math ::

\mathrm{erf2}(x, w) = y \quad \Rightarrow \quad \mathrm{erf2inv}(x, y) = w

Examples

========

>>> from sympy import I, oo, erf2inv, erfinv, erfcinv

>>> from sympy.abc import x, y

Several special values are known:

>>> erf2inv(0, 0)

>>> erf2inv(1, 0)

>>> erf2inv(0, 1)

>>> erf2inv(0, y)

erfinv(y)

>>> erf2inv(oo, y)

erfcinv(-y)

Differentiation with respect to x and y is supported:

>>> from sympy import diff

>>> diff(erf2inv(x, y), x)

exp(-x**2 + erf2inv(x, y)**2)

>>> diff(erf2inv(x, y), y)

sqrt(pi)*exp(erf2inv(x, y)**2)/2

See Also

========

erf: Gaussian error function.

erfc: Complementary error function.

erfi: Imaginary error function.

erf2: Two-argument error function.

erfinv: Inverse error function.

erfcinv: Inverse complementary error function.

References

==========

.. [1] http://functions.wolfram.com/GammaBetaErf/InverseErf2/

"""

def fdiff(self, argindex):

x, y = self.args

if argindex == 1:

return exp(self.func(x,y)**2-x**2)

elif argindex == 2:

return sqrt(S.Pi)*S.Half*exp(self.func(x,y)**2)

else:

raise ArgumentIndexError(self, argindex)

@classmethod

def eval(cls, x, y):

if x is S.NaN or y is S.NaN:

return S.NaN

elif x is S.Zero and y is S.Zero:

return S.Zero

elif x is S.Zero and y is S.One:

return S.Infinity

elif x is S.One and y is S.Zero:

return S.One

elif x is S.Zero:

return erfinv(y)

elif x is S.Infinity:

return erfcinv(-y)

elif y is S.Zero:

return x

elif y is S.Infinity:

return erfinv(x)

###############################################################################

#################### EXPONENTIAL INTEGRALS ####################################

###############################################################################

class Ei(Function):

r"""

The classical exponential integral.

For use in SymPy, this function is defined as

.. math:: \operatorname{Ei}(x) = \sum_{n=1}^\infty \frac{x^n}{n\, n!}

+ \log(x) + \gamma,

where `\gamma` is the Euler-Mascheroni constant.

If `x` is a polar number, this defines an analytic function on the

Riemann surface of the logarithm. Otherwise this defines an analytic

function in the cut plane `\mathbb{C} \setminus (-\infty, 0]`.

**Background**

The name *exponential integral* comes from the following statement:

.. math:: \operatorname{Ei}(x) = \int_{-\infty}^x \frac{e^t}{t} \mathrm{d}t

If the integral is interpreted as a Cauchy principal value, this statement

holds for `x > 0` and `\operatorname{Ei}(x)` as defined above.

Note that we carefully avoided defining `\operatorname{Ei}(x)` for

negative real `x`. This is because above integral formula does not hold for

any polar lift of such `x`, indeed all branches of

`\operatorname{Ei}(x)` above the negative reals are imaginary.

However, the following statement holds for all `x \in \mathbb{R}^*`:

.. math:: \int_{-\infty}^x \frac{e^t}{t} \mathrm{d}t =

\frac{\operatorname{Ei}\left(|x|e^{i \arg(x)}\right) +

\operatorname{Ei}\left(|x|e^{- i \arg(x)}\right)}{2},

where the integral is again understood to be a principal value if

`x > 0`, and `|x|e^{i \arg(x)}`,

`|x|e^{- i \arg(x)}` denote two conjugate polar lifts of `x`.

Examples

========

>>> from sympy import Ei, polar_lift, exp_polar, I, pi

>>> from sympy.abc import x

The exponential integral in SymPy is strictly undefined for negative values

of the argument. For convenience, exponential integrals with negative

arguments are immediately converted into an expression that agrees with

the classical integral definition:

>>> Ei(-1)

-I*pi + Ei(exp_polar(I*pi))

This yields a real value:

>>> Ei(-1).n(chop=True)

-0.219383934395520

On the other hand the analytic continuation is not real:

>>> Ei(polar_lift(-1)).n(chop=True)

-0.21938393439552 + 3.14159265358979*I

The exponential integral has a logarithmic branch point at the origin:

>>> Ei(x*exp_polar(2*I*pi))

Ei(x) + 2*I*pi

Differentiation is supported:

>>> Ei(x).diff(x)

exp(x)/x

The exponential integral is related to many other special functions.

For example:

>>> from sympy import uppergamma, expint, Shi

>>> Ei(x).rewrite(expint)

-expint(1, x*exp_polar(I*pi)) - I*pi

>>> Ei(x).rewrite(Shi)

Chi(x) + Shi(x)

See Also

========

expint: Generalised exponential integral.

E1: Special case of the generalised exponential integral.

li: Logarithmic integral.

Li: Offset logarithmic integral.

Si: Sine integral.

Ci: Cosine integral.

Shi: Hyperbolic sine integral.

Chi: Hyperbolic cosine integral.

sympy.functions.special.gamma_functions.uppergamma: Upper incomplete gamma function.

References

==========

.. [1] http://dlmf.nist.gov/6.6

.. [2] http://en.wikipedia.org/wiki/Exponential_integral

.. [3] Abramowitz & Stegun, section 5: http://people.math.sfu.ca/~cbm/aands/page_228.htm

"""

@classmethod

def eval(cls, z):

if z is S.Zero:

return S.NegativeInfinity

elif z is S.Infinity:

return S.Infinity

elif z is S.NegativeInfinity:

return S.Zero

if not z.is_polar and z.is_negative:

# Note: is this a good idea?

return Ei(polar_lift(z)) - pi*I

nz, n = z.extract_branch_factor()

if n:

return Ei(nz) + 2*I*pi*n

def fdiff(self, argindex=1):

from sympy import unpolarify

arg = unpolarify(self.args[0])

if argindex == 1:

return exp(arg)/arg

else:

raise ArgumentIndexError(self, argindex)

def _eval_evalf(self, prec):

if (self.args[0]/polar_lift(-1)).is_positive:

return Function._eval_evalf(self, prec) + (I*pi)._eval_evalf(prec)

return Function._eval_evalf(self, prec)

def _eval_rewrite_as_uppergamma(self, z):

from sympy import uppergamma

# XXX this does not currently work usefully because uppergamma

# immediately turns into expint

return -uppergamma(0, polar_lift(-1)*z) - I*pi

def _eval_rewrite_as_expint(self, z):

return -expint(1, polar_lift(-1)*z) - I*pi

def _eval_rewrite_as_li(self, z):

if isinstance(z, log):

return li(z.args[0])

# TODO:

# Actually it only holds that:

# Ei(z) = li(exp(z))

# for -pi < imag(z) <= pi

return li(exp(z))

def _eval_rewrite_as_Si(self, z):

return Shi(z) + Chi(z)

_eval_rewrite_as_Ci = _eval_rewrite_as_Si

_eval_rewrite_as_Chi = _eval_rewrite_as_Si

_eval_rewrite_as_Shi = _eval_rewrite_as_Si

def _eval_rewrite_as_tractable(self, z):

return exp(z) * _eis(z)

def _eval_nseries(self, x, n, logx):

x0 = self.args[0].limit(x, 0)

if x0 is S.Zero:

f = self._eval_rewrite_as_Si(*self.args)

return f._eval_nseries(x, n, logx)

return super(Ei, self)._eval_nseries(x, n, logx)

class expint(Function):

r"""

Generalized exponential integral.

This function is defined as

.. math:: \operatorname{E}_\nu(z) = z^{\nu - 1} \Gamma(1 - \nu, z),

where `\Gamma(1 - \nu, z)` is the upper incomplete gamma function

(``uppergamma``).

Hence for :math:`z` with positive real part we have

.. math:: \operatorname{E}_\nu(z)

= \int_1^\infty \frac{e^{-zt}}{z^\nu} \mathrm{d}t,

which explains the name.

The representation as an incomplete gamma function provides an analytic

continuation for :math:`\operatorname{E}_\nu(z)`. If :math:`\nu` is a

non-positive integer the exponential integral is thus an unbranched

function of :math:`z`, otherwise there is a branch point at the origin.

Refer to the incomplete gamma function documentation for details of the

branching behavior.

Examples

========

>>> from sympy import expint, S

>>> from sympy.abc import nu, z

Differentiation is supported. Differentiation with respect to z explains

further the name: for integral orders, the exponential integral is an

iterated integral of the exponential function.

>>> expint(nu, z).diff(z)

-expint(nu - 1, z)

Differentiation with respect to nu has no classical expression:

>>> expint(nu, z).diff(nu)

-z**(nu - 1)*meijerg(((), (1, 1)), ((0, 0, -nu + 1), ()), z)

At non-postive integer orders, the exponential integral reduces to the

exponential function:

>>> expint(0, z)

exp(-z)/z

>>> expint(-1, z)

exp(-z)/z + exp(-z)/z**2

At half-integers it reduces to error functions:

>>> expint(S(1)/2, z)

sqrt(pi)*erfc(sqrt(z))/sqrt(z)

At positive integer orders it can be rewritten in terms of exponentials

and expint(1, z). Use expand_func() to do this:

>>> from sympy import expand_func

>>> expand_func(expint(5, z))

z**4*expint(1, z)/24 + (-z**3 + z**2 - 2*z + 6)*exp(-z)/24

The generalised exponential integral is essentially equivalent to the

incomplete gamma function:

>>> from sympy import uppergamma

>>> expint(nu, z).rewrite(uppergamma)

z**(nu - 1)*uppergamma(-nu + 1, z)

As such it is branched at the origin:

>>> from sympy import exp_polar, pi, I

>>> expint(4, z*exp_polar(2*pi*I))

I*pi*z**3/3 + expint(4, z)

>>> expint(nu, z*exp_polar(2*pi*I))

z**(nu - 1)*(exp(2*I*pi*nu) - 1)*gamma(-nu + 1) + expint(nu, z)

See Also

========

Ei: Another related function called exponential integral.

E1: The classical case, returns expint(1, z).

li: Logarithmic integral.

Li: Offset logarithmic integral.

Si: Sine integral.

Ci: Cosine integral.

Shi: Hyperbolic sine integral.

Chi: Hyperbolic cosine integral.

sympy.functions.special.gamma_functions.uppergamma

References

==========

.. [1] http://dlmf.nist.gov/8.19

.. [2] http://functions.wolfram.com/GammaBetaErf/ExpIntegralE/

.. [3] http://en.wikipedia.org/wiki/Exponential_integral

"""

@classmethod

def eval(cls, nu, z):

from sympy import (unpolarify, expand_mul, uppergamma, exp, gamma,

factorial)

nu2 = unpolarify(nu)

if nu != nu2:

return expint(nu2, z)

if nu.is_Integer and nu <= 0 or (not nu.is_Integer and (2*nu).is_Integer):

return unpolarify(expand_mul(z**(nu - 1)*uppergamma(1 - nu, z)))

# Extract branching information. This can be deduced from what is

# explained in lowergamma.eval().

z, n = z.extract_branch_factor()

if n == 0:

return

if nu.is_integer:

if (nu > 0) != True:

return

return expint(nu, z) \

- 2*pi*I*n*(-1)**(nu - 1)/factorial(nu - 1)*unpolarify(z)**(nu - 1)

else:

return (exp(2*I*pi*nu*n) - 1)*z**(nu - 1)*gamma(1 - nu) + expint(nu, z)

def fdiff(self, argindex):

from sympy import meijerg

nu, z = self.args

if argindex == 1:

return -z**(nu - 1)*meijerg([], [1, 1], [0, 0, 1 - nu], [], z)

elif argindex == 2:

return -expint(nu - 1, z)

else:

raise ArgumentIndexError(self, argindex)

def _eval_rewrite_as_uppergamma(self, nu, z):

from sympy import uppergamma

return z**(nu - 1)*uppergamma(1 - nu, z)

def _eval_rewrite_as_Ei(self, nu, z):

from sympy import exp_polar, unpolarify, exp, factorial

if nu == 1:

return -Ei(z*exp_polar(-I*pi)) - I*pi

elif nu.is_Integer and nu > 1:

# DLMF, 8.19.7

x = -unpolarify(z)

return x**(nu - 1)/factorial(nu - 1)*E1(z).rewrite(Ei) + \

exp(x)/factorial(nu - 1) * \

Add(*[factorial(nu - k - 2)*x**k for k in range(nu - 1)])

else:

return self

def _eval_expand_func(self, **hints):

return self.rewrite(Ei).rewrite(expint, **hints)

def _eval_rewrite_as_Si(self, nu, z):

if nu != 1:

return self

return Shi(z) - Chi(z)

_eval_rewrite_as_Ci = _eval_rewrite_as_Si

_eval_rewrite_as_Chi = _eval_rewrite_as_Si

_eval_rewrite_as_Shi = _eval_rewrite_as_Si

def _eval_nseries(self, x, n, logx):

if not self.args[0].has(x):

nu = self.args[0]

if nu == 1:

f = self._eval_rewrite_as_Si(*self.args)

return f._eval_nseries(x, n, logx)

elif nu.is_Integer and nu > 1:

f = self._eval_rewrite_as_Ei(*self.args)

return f._eval_nseries(x, n, logx)

return super(expint, self)._eval_nseries(x, n, logx)

def _sage_(self):

import sage.all as sage

return sage.exp_integral_e(self.args[0]._sage_(), self.args[1]._sage_())

def E1(z):

"""

Classical case of the generalized exponential integral.

This is equivalent to ``expint(1, z)``.

See Also

========

Ei: Exponential integral.

expint: Generalised exponential integral.

li: Logarithmic integral.

Li: Offset logarithmic integral.

Si: Sine integral.

Ci: Cosine integral.

Shi: Hyperbolic sine integral.

Chi: Hyperbolic cosine integral.

"""

return expint(1, z)

class li(Function):

r"""

The classical logarithmic integral.

For the use in SymPy, this function is defined as

.. math:: \operatorname{li}(x) = \int_0^x \frac{1}{\log(t)} \mathrm{d}t \,.

Examples

========

>>> from sympy import I, oo, li

>>> from sympy.abc import z

Several special values are known:

>>> li(0)

>>> li(1)

-oo

>>> li(oo)

Differentiation with respect to z is supported:

>>> from sympy import diff

>>> diff(li(z), z)

1/log(z)

Defining the `li` function via an integral:

The logarithmic integral can also be defined in terms of Ei:

>>> from sympy import Ei

>>> li(z).rewrite(Ei)

Ei(log(z))

>>> diff(li(z).rewrite(Ei), z)

1/log(z)

We can numerically evaluate the logarithmic integral to arbitrary precision

on the whole complex plane (except the singular points):

>>> li(2).evalf(30)

1.04516378011749278484458888919

>>> li(2*I).evalf(30)

1.0652795784357498247001125598 + 3.08346052231061726610939702133*I

We can even compute Soldner's constant by the help of mpmath:

>>> from mpmath import findroot

>>> findroot(li, 2)

1.45136923488338

Further transformations include rewriting `li` in terms of

the trigonometric integrals `Si`, `Ci`, `Shi` and `Chi`:

>>> from sympy import Si, Ci, Shi, Chi

>>> li(z).rewrite(Si)

-log(I*log(z)) - log(1/log(z))/2 + log(log(z))/2 + Ci(I*log(z)) + Shi(log(z))

>>> li(z).rewrite(Ci)

-log(I*log(z)) - log(1/log(z))/2 + log(log(z))/2 + Ci(I*log(z)) + Shi(log(z))

>>> li(z).rewrite(Shi)

-log(1/log(z))/2 + log(log(z))/2 + Chi(log(z)) - Shi(log(z))

>>> li(z).rewrite(Chi)

-log(1/log(z))/2 + log(log(z))/2 + Chi(log(z)) - Shi(log(z))

See Also

========

Li: Offset logarithmic integral.

Ei: Exponential integral.

expint: Generalised exponential integral.

E1: Special case of the generalised exponential integral.

Si: Sine integral.

Ci: Cosine integral.

Shi: Hyperbolic sine integral.

Chi: Hyperbolic cosine integral.

References

==========

.. [1] http://en.wikipedia.org/wiki/Logarithmic_integral

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

.. [3] http://dlmf.nist.gov/6

.. [4] http://mathworld.wolfram.com/SoldnersConstant.html

"""

@classmethod

def eval(cls, z):

if z is S.Zero:

return S.Zero

elif z is S.One:

return S.NegativeInfinity

elif z is S.Infinity:

return S.Infinity

def fdiff(self, argindex=1):

arg = self.args[0]

if argindex == 1:

return S.One / log(arg)

else:

raise ArgumentIndexError(self, argindex)

def _eval_conjugate(self):

z = self.args[0]

# Exclude values on the branch cut (-oo, 0)

if not (z.is_real and z.is_negative):

return self.func(z.conjugate())

def _eval_rewrite_as_Li(self, z):

return Li(z) + li(2)

def _eval_rewrite_as_Ei(self, z):

return Ei(log(z))

def _eval_rewrite_as_uppergamma(self, z):

from sympy import uppergamma

return (-uppergamma(0, -log(z)) +

S.Half*(log(log(z)) - log(S.One/log(z))) - log(-log(z)))

def _eval_rewrite_as_Si(self, z):

return (Ci(I*log(z)) - I*Si(I*log(z)) -

S.Half*(log(S.One/log(z)) - log(log(z))) - log(I*log(z)))

_eval_rewrite_as_Ci = _eval_rewrite_as_Si

def _eval_rewrite_as_Shi(self, z):

return (Chi(log(z)) - Shi(log(z)) - S.Half*(log(S.One/log(z)) - log(log(z))))

_eval_rewrite_as_Chi = _eval_rewrite_as_Shi

def _eval_rewrite_as_hyper(self, z):

return (log(z)*hyper((1, 1), (2, 2), log(z)) +

S.Half*(log(log(z)) - log(S.One/log(z))) + S.EulerGamma)

def _eval_rewrite_as_meijerg(self, z):

return (-log(-log(z)) - S.Half*(log(S.One/log(z)) - log(log(z)))

- meijerg(((), (1,)), ((0, 0), ()), -log(z)))

def _eval_rewrite_as_tractable(self, z):

return z * _eis(log(z))

class Li(Function):

r"""

The offset logarithmic integral.

For the use in SymPy, this function is defined as

.. math:: \operatorname{Li}(x) = \operatorname{li}(x) - \operatorname{li}(2)

Examples

========

>>> from sympy import I, oo, Li

>>> from sympy.abc import z

The following special value is known:

>>> Li(2)

Differentiation with respect to z is supported:

>>> from sympy import diff

>>> diff(Li(z), z)

1/log(z)

The shifted logarithmic integral can be written in terms of `li(z)`:

>>> from sympy import li

>>> Li(z).rewrite(li)

li(z) - li(2)

We can numerically evaluate the logarithmic integral to arbitrary precision

on the whole complex plane (except the singular points):

>>> Li(2).evalf(30)

>>> Li(4).evalf(30)

1.92242131492155809316615998938

See Also

========

li: Logarithmic integral.

Ei: Exponential integral.

expint: Generalised exponential integral.

E1: Special case of the generalised exponential integral.

Si: Sine integral.

Ci: Cosine integral.

Shi: Hyperbolic sine integral.

Chi: Hyperbolic cosine integral.

References

==========

.. [1] http://en.wikipedia.org/wiki/Logarithmic_integral

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

.. [3] http://dlmf.nist.gov/6

"""

@classmethod

def eval(cls, z):

if z is S.Infinity:

return S.Infinity

elif z is 2*S.One:

return S.Zero

def fdiff(self, argindex=1):

arg = self.args[0]

if argindex == 1:

return S.One / log(arg)

else:

raise ArgumentIndexError(self, argindex)

def _eval_evalf(self, prec):

return self.rewrite(li).evalf(prec)

def _eval_rewrite_as_li(self, z):

return li(z) - li(2)

def _eval_rewrite_as_tractable(self, z):

return self.rewrite(li).rewrite("tractable", deep=True)

###############################################################################

#################### TRIGONOMETRIC INTEGRALS ##################################

###############################################################################

class TrigonometricIntegral(Function):

""" Base class for trigonometric integrals. """

@classmethod

def eval(cls, z):

if z == 0:

return cls._atzero

elif z is S.Infinity:

return cls._atinf()

elif z is S.NegativeInfinity:

return cls._atneginf()

nz = z.extract_multiplicatively(polar_lift(I))

if nz is None and cls._trigfunc(0) == 0:

nz = z.extract_multiplicatively(I)

if nz is not None:

return cls._Ifactor(nz, 1)

nz = z.extract_multiplicatively(polar_lift(-I))

if nz is not None:

return cls._Ifactor(nz, -1)

nz = z.extract_multiplicatively(polar_lift(-1))

if nz is None and cls._trigfunc(0) == 0:

nz = z.extract_multiplicatively(-1)

if nz is not None:

return cls._minusfactor(nz)

nz, n = z.extract_branch_factor()

if n == 0 and nz == z:

return

return 2*pi*I*n*cls._trigfunc(0) + cls(nz)

def fdiff(self, argindex=1):

from sympy import unpolarify

arg = unpolarify(self.args[0])

if argindex == 1:

return self._trigfunc(arg)/arg

def _eval_rewrite_as_Ei(self, z):

return self._eval_rewrite_as_expint(z).rewrite(Ei)

def _eval_rewrite_as_uppergamma(self, z):

from sympy import uppergamma

return self._eval_rewrite_as_expint(z).rewrite(uppergamma)

def _eval_nseries(self, x, n, logx):

# NOTE this is fairly inefficient

from sympy import log, EulerGamma, Pow

n += 1

if self.args[0].subs(x, 0) != 0:

return super(TrigonometricIntegral, self)._eval_nseries(x, n, logx)

baseseries = self._trigfunc(x)._eval_nseries(x, n, logx)

if self._trigfunc(0) != 0:

baseseries -= 1

baseseries = baseseries.replace(Pow, lambda t, n: t**n/n, simultaneous=False)

if self._trigfunc(0) != 0:

baseseries += EulerGamma + log(x)

return baseseries.subs(x, self.args[0])._eval_nseries(x, n, logx)

class Si(TrigonometricIntegral):

r"""

Sine integral.

This function is defined by

.. math:: \operatorname{Si}(z) = \int_0^z \frac{\sin{t}}{t} \mathrm{d}t.

It is an entire function.

Examples

========

>>> from sympy import Si

>>> from sympy.abc import z

The sine integral is an antiderivative of sin(z)/z:

>>> Si(z).diff(z)

sin(z)/z

It is unbranched:

>>> from sympy import exp_polar, I, pi

>>> Si(z*exp_polar(2*I*pi))

Si(z)

Sine integral behaves much like ordinary sine under multiplication by ``I``:

>>> Si(I*z)

I*Shi(z)

>>> Si(-z)

-Si(z)

It can also be expressed in terms of exponential integrals, but beware

that the latter is branched:

>>> from sympy import expint

>>> Si(z).rewrite(expint)

-I*(-expint(1, z*exp_polar(-I*pi/2))/2 +

expint(1, z*exp_polar(I*pi/2))/2) + pi/2

It can be rewritten in the form of sinc function (By definition)

>>> from sympy import sinc

>>> Si(z).rewrite(sinc)

Integral(sinc(t), (t, 0, z))

See Also

========

Ci: Cosine integral.

Shi: Hyperbolic sine integral.

Chi: Hyperbolic cosine integral.

Ei: Exponential integral.

expint: Generalised exponential integral.

sinc: unnormalized sinc function

E1: Special case of the generalised exponential integral.

li: Logarithmic integral.

Li: Offset logarithmic integral.

References

==========

.. [1] http://en.wikipedia.org/wiki/Trigonometric_integral

"""

_trigfunc = sin

_atzero = S(0)

@classmethod

def _atinf(cls):

return pi*S.Half

@classmethod

def _atneginf(cls):

return -pi*S.Half

@classmethod

def _minusfactor(cls, z):

return -Si(z)

@classmethod

def _Ifactor(cls, z, sign):

return I*Shi(z)*sign

def _eval_rewrite_as_expint(self, z):

# XXX should we polarify z?

return pi/2 + (E1(polar_lift(I)*z) - E1(polar_lift(-I)*z))/2/I

def _eval_rewrite_as_sinc(self, z):

from sympy import Integral

t = Symbol('t', Dummy=True)

return Integral(sinc(t), (t, 0, z))

def _sage_(self):

import sage.all as sage

return sage.sin_integral(self.args[0]._sage_())

class Ci(TrigonometricIntegral):

r"""

Cosine integral.

This function is defined for positive `x` by

.. math:: \operatorname{Ci}(x) = \gamma + \log{x}

+ \int_0^x \frac{\cos{t} - 1}{t} \mathrm{d}t

= -\int_x^\infty \frac{\cos{t}}{t} \mathrm{d}t,

where `\gamma` is the Euler-Mascheroni constant.

We have

.. math:: \operatorname{Ci}(z) =

-\frac{\operatorname{E}_1\left(e^{i\pi/2} z\right)

+ \operatorname{E}_1\left(e^{-i \pi/2} z\right)}{2}

which holds for all polar `z` and thus provides an analytic

continuation to the Riemann surface of the logarithm.

The formula also holds as stated

for `z \in \mathbb{C}` with `\Re(z) > 0`.

By lifting to the principal branch we obtain an analytic function on the

cut complex plane.

Examples

========

>>> from sympy import Ci

>>> from sympy.abc import z

The cosine integral is a primitive of `\cos(z)/z`:

>>> Ci(z).diff(z)

cos(z)/z

It has a logarithmic branch point at the origin:

>>> from sympy import exp_polar, I, pi

>>> Ci(z*exp_polar(2*I*pi))

Ci(z) + 2*I*pi

The cosine integral behaves somewhat like ordinary `\cos` under multiplication by `i`:

>>> from sympy import polar_lift

>>> Ci(polar_lift(I)*z)

Chi(z) + I*pi/2

>>> Ci(polar_lift(-1)*z)

Ci(z) + I*pi

It can also be expressed in terms of exponential integrals:

>>> from sympy import expint

>>> Ci(z).rewrite(expint)

-expint(1, z*exp_polar(-I*pi/2))/2 - expint(1, z*exp_polar(I*pi/2))/2

See Also

========

Si: Sine integral.

Shi: Hyperbolic sine integral.

Chi: Hyperbolic cosine integral.

Ei: Exponential integral.

expint: Generalised exponential integral.

E1: Special case of the generalised exponential integral.

li: Logarithmic integral.

Li: Offset logarithmic integral.

References

==========

.. [1] http://en.wikipedia.org/wiki/Trigonometric_integral

"""

_trigfunc = cos

_atzero = S.ComplexInfinity

@classmethod

def _atinf(cls):

return S.Zero

@classmethod

def _atneginf(cls):

return I*pi

@classmethod

def _minusfactor(cls, z):

return Ci(z) + I*pi

@classmethod

def _Ifactor(cls, z, sign):

return Chi(z) + I*pi/2*sign

def _eval_rewrite_as_expint(self, z):

return -(E1(polar_lift(I)*z) + E1(polar_lift(-I)*z))/2

def _sage_(self):

import sage.all as sage

return sage.cos_integral(self.args[0]._sage_())

class Shi(TrigonometricIntegral):

r"""

Sinh integral.

This function is defined by

.. math:: \operatorname{Shi}(z) = \int_0^z \frac{\sinh{t}}{t} \mathrm{d}t.

It is an entire function.

Examples

========

>>> from sympy import Shi

>>> from sympy.abc import z

The Sinh integral is a primitive of `\sinh(z)/z`:

>>> Shi(z).diff(z)

sinh(z)/z

It is unbranched:

>>> from sympy import exp_polar, I, pi

>>> Shi(z*exp_polar(2*I*pi))

Shi(z)

The `\sinh` integral behaves much like ordinary `\sinh` under multiplication by `i`:

>>> Shi(I*z)

I*Si(z)

>>> Shi(-z)

-Shi(z)

It can also be expressed in terms of exponential integrals, but beware

that the latter is branched:

>>> from sympy import expint

>>> Shi(z).rewrite(expint)

expint(1, z)/2 - expint(1, z*exp_polar(I*pi))/2 - I*pi/2

See Also

========

Si: Sine integral.

Ci: Cosine integral.

Chi: Hyperbolic cosine integral.

Ei: Exponential integral.

expint: Generalised exponential integral.

E1: Special case of the generalised exponential integral.

li: Logarithmic integral.

Li: Offset logarithmic integral.

References

==========

.. [1] http://en.wikipedia.org/wiki/Trigonometric_integral

"""

_trigfunc = sinh

_atzero = S(0)

@classmethod

def _atinf(cls):

return S.Infinity

@classmethod

def _atneginf(cls):

return S.NegativeInfinity

@classmethod

def _minusfactor(cls, z):

return -Shi(z)

@classmethod

def _Ifactor(cls, z, sign):

return I*Si(z)*sign

def _eval_rewrite_as_expint(self, z):

from sympy import exp_polar

# XXX should we polarify z?

return (E1(z) - E1(exp_polar(I*pi)*z))/2 - I*pi/2

def _sage_(self):

import sage.all as sage

return sage.sinh_integral(self.args[0]._sage_())

class Chi(TrigonometricIntegral):

r"""

Cosh integral.

This function is defined for positive :math:`x` by

.. math:: \operatorname{Chi}(x) = \gamma + \log{x}

+ \int_0^x \frac{\cosh{t} - 1}{t} \mathrm{d}t,

where :math:`\gamma` is the Euler-Mascheroni constant.

We have

.. math:: \operatorname{Chi}(z) = \operatorname{Ci}\left(e^{i \pi/2}z\right)

- i\frac{\pi}{2},

which holds for all polar :math:`z` and thus provides an analytic

continuation to the Riemann surface of the logarithm.

By lifting to the principal branch we obtain an analytic function on the

cut complex plane.

Examples

========

>>> from sympy import Chi

>>> from sympy.abc import z

The `\cosh` integral is a primitive of `\cosh(z)/z`:

>>> Chi(z).diff(z)

cosh(z)/z

It has a logarithmic branch point at the origin:

>>> from sympy import exp_polar, I, pi

>>> Chi(z*exp_polar(2*I*pi))

Chi(z) + 2*I*pi

The `\cosh` integral behaves somewhat like ordinary `\cosh` under multiplication by `i`:

>>> from sympy import polar_lift

>>> Chi(polar_lift(I)*z)

Ci(z) + I*pi/2

>>> Chi(polar_lift(-1)*z)

Chi(z) + I*pi

It can also be expressed in terms of exponential integrals:

>>> from sympy import expint

>>> Chi(z).rewrite(expint)

-expint(1, z)/2 - expint(1, z*exp_polar(I*pi))/2 - I*pi/2

See Also

========

Si: Sine integral.

Ci: Cosine integral.

Shi: Hyperbolic sine integral.

Ei: Exponential integral.

expint: Generalised exponential integral.

E1: Special case of the generalised exponential integral.

li: Logarithmic integral.

Li: Offset logarithmic integral.

References

==========

.. [1] http://en.wikipedia.org/wiki/Trigonometric_integral

"""

_trigfunc = cosh

_atzero = S.ComplexInfinity

@classmethod

def _atinf(cls):

return S.Infinity

@classmethod

def _atneginf(cls):

return S.Infinity

@classmethod

def _minusfactor(cls, z):

return Chi(z) + I*pi

@classmethod

def _Ifactor(cls, z, sign):

return Ci(z) + I*pi/2*sign

def _eval_rewrite_as_expint(self, z):

from sympy import exp_polar

return -I*pi/2 - (E1(z) + E1(exp_polar(I*pi)*z))/2

def _latex(self, printer, exp=None):

if len(self.args) != 1:

raise ValueError("Arg length should be 1")

if exp:

return r'\operatorname{Chi}^{%s}{\left (%s \right )}' \

% (printer._print(exp), printer._print(self.args[0]))

else:

return r'\operatorname{Chi}{\left (%s \right )}' \

% printer._print(self.args[0])

@staticmethod

def _latex_no_arg(printer):

return r'\operatorname{Chi}'

def _sage_(self):

import sage.all as sage

return sage.cosh_integral(self.args[0]._sage_())

###############################################################################

#################### FRESNEL INTEGRALS ########################################

###############################################################################

class FresnelIntegral(Function):

""" Base class for the Fresnel integrals."""

unbranched = True

@classmethod

def eval(cls, z):

# Value at zero

if z is S.Zero:

return S(0)

# Try to pull out factors of -1 and I

prefact = S.One

newarg = z

changed = False

nz = newarg.extract_multiplicatively(-1)

if nz is not None:

prefact = -prefact

newarg = nz

changed = True

nz = newarg.extract_multiplicatively(I)

if nz is not None:

prefact = cls._sign*I*prefact

newarg = nz

changed = True

if changed:

return prefact*cls(newarg)

# Values at positive infinities signs

# if any were extracted automatically

if z is S.Infinity:

return S.Half

elif z is I*S.Infinity:

return cls._sign*I*S.Half

def fdiff(self, argindex=1):

if argindex == 1:

return self._trigfunc(S.Half*pi*self.args[0]**2)

else:

raise ArgumentIndexError(self, argindex)

def _eval_is_real(self):

return self.args[0].is_real

def _eval_conjugate(self):

return self.func(self.args[0].conjugate())

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:

re, im = self.args[0].expand(deep, **hints).as_real_imag()

else:

re, im = self.args[0].as_real_imag()

return (re, im)

def as_real_imag(self, deep=True, **hints):

# Fresnel S

# http://functions.wolfram.com/06.32.19.0003.01

# http://functions.wolfram.com/06.32.19.0006.01

# Fresnel C

# http://functions.wolfram.com/06.33.19.0003.01

# http://functions.wolfram.com/06.33.19.0006.01

x, y = self._as_real_imag(deep=deep, **hints)

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)

class fresnels(FresnelIntegral):

r"""

Fresnel integral S.

This function is defined by

.. math:: \operatorname{S}(z) = \int_0^z \sin{\frac{\pi}{2} t^2} \mathrm{d}t.

It is an entire function.

Examples

========

>>> from sympy import I, oo, fresnels

>>> from sympy.abc import z

Several special values are known:

>>> fresnels(0)

>>> fresnels(oo)

1/2

>>> fresnels(-oo)

-1/2

>>> fresnels(I*oo)

-I/2

>>> fresnels(-I*oo)

I/2

In general one can pull out factors of -1 and `i` from the argument:

>>> fresnels(-z)

-fresnels(z)

>>> fresnels(I*z)

-I*fresnels(z)

The Fresnel S integral obeys the mirror symmetry

`\overline{S(z)} = S(\bar{z})`:

>>> from sympy import conjugate

>>> conjugate(fresnels(z))

fresnels(conjugate(z))

Differentiation with respect to `z` is supported:

>>> from sympy import diff

>>> diff(fresnels(z), z)

sin(pi*z**2/2)

Defining the Fresnel functions via an integral

>>> from sympy import integrate, pi, sin, gamma, expand_func

>>> integrate(sin(pi*z**2/2), z)

3*fresnels(z)*gamma(3/4)/(4*gamma(7/4))

>>> expand_func(integrate(sin(pi*z**2/2), z))

fresnels(z)

We can numerically evaluate the Fresnel integral to arbitrary precision

on the whole complex plane:

>>> fresnels(2).evalf(30)

0.343415678363698242195300815958

>>> fresnels(-2*I).evalf(30)

0.343415678363698242195300815958*I

See Also

========

fresnelc: Fresnel cosine integral.

References

==========

.. [1] http://en.wikipedia.org/wiki/Fresnel_integral

.. [2] http://dlmf.nist.gov/7

.. [3] http://mathworld.wolfram.com/FresnelIntegrals.html

.. [4] http://functions.wolfram.com/GammaBetaErf/FresnelS

.. [5] The converging factors for the fresnel integrals

by John W. Wrench Jr. and Vicki Alley

"""

_trigfunc = sin

_sign = -S.One

@staticmethod

@cacheit

def taylor_term(n, x, *previous_terms):

if n < 0:

return S.Zero

else:

x = sympify(x)

if len(previous_terms) > 1:

p = previous_terms[-1]

return (-pi**2*x**4*(4*n - 1)/(8*n*(2*n + 1)*(4*n + 3))) * p

else:

return x**3 * (-x**4)**n * (S(2)**(-2*n - 1)*pi**(2*n + 1)) / ((4*n + 3)*factorial(2*n + 1))

def _eval_rewrite_as_erf(self, z):

return (S.One + I)/4 * (erf((S.One + I)/2*sqrt(pi)*z) - I*erf((S.One - I)/2*sqrt(pi)*z))

def _eval_rewrite_as_hyper(self, z):

return pi*z**3/6 * hyper([S(3)/4], [S(3)/2, S(7)/4], -pi**2*z**4/16)

def _eval_rewrite_as_meijerg(self, z):

return (pi*z**(S(9)/4) / (sqrt(2)*(z**2)**(S(3)/4)*(-z)**(S(3)/4))

* meijerg([], [1], [S(3)/4], [S(1)/4, 0], -pi**2*z**4/16))

def _eval_aseries(self, n, args0, x, logx):

from sympy import Order

point = args0[0]

# Expansion at oo

if point is S.Infinity:

z = self.args[0]

# expansion of S(x) = S1(x*sqrt(pi/2)), see reference[5] page 1-8

p = [(-1)**k * factorial(4*k + 1) /

(2**(2*k + 2) * z**(4*k + 3) * 2**(2*k)*factorial(2*k))

for k in range(0, n)]

q = [1/(2*z)] + [(-1)**k * factorial(4*k - 1) /

(2**(2*k + 1) * z**(4*k + 1) * 2**(2*k - 1)*factorial(2*k - 1))

for k in range(1, n)]

p = [-sqrt(2/pi)*t for t in p] + [Order(1/z**n, x)]

q = [-sqrt(2/pi)*t for t in q] + [Order(1/z**n, x)]

return S.Half + (sin(z**2)*Add(*p) + cos(z**2)*Add(*q)).subs(x, sqrt(2/pi)*x)

# All other points are not handled

return super(fresnels, self)._eval_aseries(n, args0, x, logx)

class fresnelc(FresnelIntegral):

r"""

Fresnel integral C.

This function is defined by

.. math:: \operatorname{C}(z) = \int_0^z \cos{\frac{\pi}{2} t^2} \mathrm{d}t.

It is an entire function.

Examples

========

>>> from sympy import I, oo, fresnelc

>>> from sympy.abc import z

Several special values are known:

>>> fresnelc(0)

>>> fresnelc(oo)

1/2

>>> fresnelc(-oo)

-1/2

>>> fresnelc(I*oo)

I/2

>>> fresnelc(-I*oo)

-I/2

In general one can pull out factors of -1 and `i` from the argument:

>>> fresnelc(-z)

-fresnelc(z)

>>> fresnelc(I*z)

I*fresnelc(z)

The Fresnel C integral obeys the mirror symmetry

`\overline{C(z)} = C(\bar{z})`:

>>> from sympy import conjugate

>>> conjugate(fresnelc(z))

fresnelc(conjugate(z))

Differentiation with respect to `z` is supported:

>>> from sympy import diff

>>> diff(fresnelc(z), z)

cos(pi*z**2/2)

Defining the Fresnel functions via an integral

>>> from sympy import integrate, pi, cos, gamma, expand_func

>>> integrate(cos(pi*z**2/2), z)

fresnelc(z)*gamma(1/4)/(4*gamma(5/4))

>>> expand_func(integrate(cos(pi*z**2/2), z))

fresnelc(z)

We can numerically evaluate the Fresnel integral to arbitrary precision

on the whole complex plane:

>>> fresnelc(2).evalf(30)

0.488253406075340754500223503357

>>> fresnelc(-2*I).evalf(30)

-0.488253406075340754500223503357*I

See Also

========

fresnels: Fresnel sine integral.

References

==========

.. [1] http://en.wikipedia.org/wiki/Fresnel_integral

.. [2] http://dlmf.nist.gov/7

.. [3] http://mathworld.wolfram.com/FresnelIntegrals.html

.. [4] http://functions.wolfram.com/GammaBetaErf/FresnelC

.. [5] The converging factors for the fresnel integrals

by John W. Wrench Jr. and Vicki Alley

"""

_trigfunc = cos

_sign = S.One

@staticmethod

@cacheit

def taylor_term(n, x, *previous_terms):

if n < 0:

return S.Zero

else:

x = sympify(x)

if len(previous_terms) > 1:

p = previous_terms[-1]

return (-pi**2*x**4*(4*n - 3)/(8*n*(2*n - 1)*(4*n + 1))) * p

else:

return x * (-x**4)**n * (S(2)**(-2*n)*pi**(2*n)) / ((4*n + 1)*factorial(2*n))

def _eval_rewrite_as_erf(self, z):

return (S.One - I)/4 * (erf((S.One + I)/2*sqrt(pi)*z) + I*erf((S.One - I)/2*sqrt(pi)*z))

def _eval_rewrite_as_hyper(self, z):

return z * hyper([S.One/4], [S.One/2, S(5)/4], -pi**2*z**4/16)

def _eval_rewrite_as_meijerg(self, z):

return (pi*z**(S(3)/4) / (sqrt(2)*root(z**2, 4)*root(-z, 4))

* meijerg([], [1], [S(1)/4], [S(3)/4, 0], -pi**2*z**4/16))

def _eval_aseries(self, n, args0, x, logx):

from sympy import Order

point = args0[0]

# Expansion at oo

if point is S.Infinity:

z = self.args[0]

# expansion of C(x) = C1(x*sqrt(pi/2)), see reference[5] page 1-8

p = [(-1)**k * factorial(4*k + 1) /

(2**(2*k + 2) * z**(4*k + 3) * 2**(2*k)*factorial(2*k))

for k in range(0, n)]

q = [1/(2*z)] + [(-1)**k * factorial(4*k - 1) /

(2**(2*k + 1) * z**(4*k + 1) * 2**(2*k - 1)*factorial(2*k - 1))

for k in range(1, n)]

p = [-sqrt(2/pi)*t for t in p] + [Order(1/z**n, x)]

q = [ sqrt(2/pi)*t for t in q] + [Order(1/z**n, x)]

return S.Half + (cos(z**2)*Add(*p) + sin(z**2)*Add(*q)).subs(x, sqrt(2/pi)*x)

# All other points are not handled

return super(fresnelc, self)._eval_aseries(n, args0, x, logx)

###############################################################################

#################### HELPER FUNCTIONS #########################################

###############################################################################

class _erfs(Function):

"""

Helper function to make the `\\mathrm{erf}(z)` function

tractable for the Gruntz algorithm.

"""

def _eval_aseries(self, n, args0, x, logx):

from sympy import Order

point = args0[0]

# Expansion at oo

if point is S.Infinity:

z = self.args[0]

l = [ 1/sqrt(S.Pi) * factorial(2*k)*(-S(

4))**(-k)/factorial(k) * (1/z)**(2*k + 1) for k in range(0, n) ]

o = Order(1/z**(2*n + 1), x)

# It is very inefficient to first add the order and then do the nseries

return (Add(*l))._eval_nseries(x, n, logx) + o

# Expansion at I*oo

t = point.extract_multiplicatively(S.ImaginaryUnit)

if t is S.Infinity:

z = self.args[0]

# TODO: is the series really correct?

l = [ 1/sqrt(S.Pi) * factorial(2*k)*(-S(

4))**(-k)/factorial(k) * (1/z)**(2*k + 1) for k in range(0, n) ]

o = Order(1/z**(2*n + 1), x)

# It is very inefficient to first add the order and then do the nseries

return (Add(*l))._eval_nseries(x, n, logx) + o

# All other points are not handled

return super(_erfs, self)._eval_aseries(n, args0, x, logx)

def fdiff(self, argindex=1):

if argindex == 1:

z = self.args[0]

return -2/sqrt(S.Pi) + 2*z*_erfs(z)

else:

raise ArgumentIndexError(self, argindex)

def _eval_rewrite_as_intractable(self, z):

return (S.One - erf(z))*exp(z**2)

class _eis(Function):

"""

Helper function to make the `\\mathrm{Ei}(z)` and `\\mathrm{li}(z)` functions

tractable for the Gruntz algorithm.

"""

def _eval_aseries(self, n, args0, x, logx):

from sympy import Order

if args0[0] != S.Infinity:

return super(_erfs, self)._eval_aseries(n, args0, x, logx)

z = self.args[0]

l = [ factorial(k) * (1/z)**(k + 1) for k in range(0, n) ]

o = Order(1/z**(n + 1), x)

# It is very inefficient to first add the order and then do the nseries

return (Add(*l))._eval_nseries(x, n, logx) + o

def fdiff(self, argindex=1):

if argindex == 1:

z = self.args[0]

return S.One / z - _eis(z)

else:

raise ArgumentIndexError(self, argindex)

def _eval_rewrite_as_intractable(self, z):

return exp(-z)*Ei(z)

def _eval_nseries(self, x, n, logx):

x0 = self.args[0].limit(x, 0)

if x0 is S.Zero:

f = self._eval_rewrite_as_intractable(*self.args)

return f._eval_nseries(x, n, logx)

return super(_eis, self)._eval_nseries(x, n, logx)

Coverage for sympy/functions/special/error_functions.py : 6%

796 statements 49 run 747 missing 0 excluded