You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

1292 lines
42 KiB

from itertools import product
from typing import Tuple as tTuple
from sympy.core.add import Add
from sympy.core.cache import cacheit
from sympy.core.expr import Expr
from sympy.core.function import (Function, ArgumentIndexError, expand_log,
expand_mul, FunctionClass, PoleError, expand_multinomial, expand_complex)
from sympy.core.logic import fuzzy_and, fuzzy_not, fuzzy_or
from sympy.core.mul import Mul
from sympy.core.numbers import Integer, Rational, pi, I, ImaginaryUnit
from sympy.core.parameters import global_parameters
from sympy.core.power import Pow
from sympy.core.singleton import S
from sympy.core.symbol import Wild, Dummy
from sympy.core.sympify import sympify
from sympy.functions.combinatorial.factorials import factorial
from sympy.functions.elementary.complexes import arg, unpolarify, im, re, Abs
from sympy.functions.elementary.miscellaneous import sqrt
from sympy.ntheory import multiplicity, perfect_power
from sympy.ntheory.factor_ import factorint
# NOTE IMPORTANT
# The series expansion code in this file is an important part of the gruntz
# algorithm for determining limits. _eval_nseries has to return a generalized
# power series with coefficients in C(log(x), log).
# In more detail, the result of _eval_nseries(self, x, n) must be
# c_0*x**e_0 + ... (finitely many terms)
# where e_i are numbers (not necessarily integers) and c_i involve only
# numbers, the function log, and log(x). [This also means it must not contain
# log(x(1+p)), this *has* to be expanded to log(x)+log(1+p) if x.is_positive and
# p.is_positive.]
class ExpBase(Function):
unbranched = True
_singularities = (S.ComplexInfinity,)
@property
def kind(self):
return self.exp.kind
def inverse(self, argindex=1):
"""
Returns the inverse function of ``exp(x)``.
"""
return log
def as_numer_denom(self):
"""
Returns this with a positive exponent as a 2-tuple (a fraction).
Examples
========
>>> from sympy import exp
>>> from sympy.abc import x
>>> exp(-x).as_numer_denom()
(1, exp(x))
>>> exp(x).as_numer_denom()
(exp(x), 1)
"""
# this should be the same as Pow.as_numer_denom wrt
# exponent handling
exp = self.exp
neg_exp = exp.is_negative
if not neg_exp and not (-exp).is_negative:
neg_exp = exp.could_extract_minus_sign()
if neg_exp:
return S.One, self.func(-exp)
return self, S.One
@property
def exp(self):
"""
Returns the exponent of the function.
"""
return self.args[0]
def as_base_exp(self):
"""
Returns the 2-tuple (base, exponent).
"""
return self.func(1), Mul(*self.args)
def _eval_adjoint(self):
return self.func(self.exp.adjoint())
def _eval_conjugate(self):
return self.func(self.exp.conjugate())
def _eval_transpose(self):
return self.func(self.exp.transpose())
def _eval_is_finite(self):
arg = self.exp
if arg.is_infinite:
if arg.is_extended_negative:
return True
if arg.is_extended_positive:
return False
if arg.is_finite:
return True
def _eval_is_rational(self):
s = self.func(*self.args)
if s.func == self.func:
z = s.exp.is_zero
if z:
return True
elif s.exp.is_rational and fuzzy_not(z):
return False
else:
return s.is_rational
def _eval_is_zero(self):
return self.exp is S.NegativeInfinity
def _eval_power(self, other):
"""exp(arg)**e -> exp(arg*e) if assumptions allow it.
"""
b, e = self.as_base_exp()
return Pow._eval_power(Pow(b, e, evaluate=False), other)
def _eval_expand_power_exp(self, **hints):
from sympy.concrete.products import Product
from sympy.concrete.summations import Sum
arg = self.args[0]
if arg.is_Add and arg.is_commutative:
return Mul.fromiter(self.func(x) for x in arg.args)
elif isinstance(arg, Sum) and arg.is_commutative:
return Product(self.func(arg.function), *arg.limits)
return self.func(arg)
class exp_polar(ExpBase):
r"""
Represent a *polar number* (see g-function Sphinx documentation).
Explanation
===========
``exp_polar`` represents the function
`Exp: \mathbb{C} \rightarrow \mathcal{S}`, sending the complex number
`z = a + bi` to the polar number `r = exp(a), \theta = b`. It is one of
the main functions to construct polar numbers.
Examples
========
>>> from sympy import exp_polar, pi, I, exp
The main difference is that polar numbers do not "wrap around" at `2 \pi`:
>>> exp(2*pi*I)
1
>>> exp_polar(2*pi*I)
exp_polar(2*I*pi)
apart from that they behave mostly like classical complex numbers:
>>> exp_polar(2)*exp_polar(3)
exp_polar(5)
See Also
========
sympy.simplify.powsimp.powsimp
polar_lift
periodic_argument
principal_branch
"""
is_polar = True
is_comparable = False # cannot be evalf'd
def _eval_Abs(self): # Abs is never a polar number
return exp(re(self.args[0]))
def _eval_evalf(self, prec):
""" Careful! any evalf of polar numbers is flaky """
i = im(self.args[0])
try:
bad = (i <= -pi or i > pi)
except TypeError:
bad = True
if bad:
return self # cannot evalf for this argument
res = exp(self.args[0])._eval_evalf(prec)
if i > 0 and im(res) < 0:
# i ~ pi, but exp(I*i) evaluated to argument slightly bigger than pi
return re(res)
return res
def _eval_power(self, other):
return self.func(self.args[0]*other)
def _eval_is_extended_real(self):
if self.args[0].is_extended_real:
return True
def as_base_exp(self):
# XXX exp_polar(0) is special!
if self.args[0] == 0:
return self, S.One
return ExpBase.as_base_exp(self)
class ExpMeta(FunctionClass):
def __instancecheck__(cls, instance):
if exp in instance.__class__.__mro__:
return True
return isinstance(instance, Pow) and instance.base is S.Exp1
class exp(ExpBase, metaclass=ExpMeta):
"""
The exponential function, :math:`e^x`.
Examples
========
>>> from sympy import exp, I, pi
>>> from sympy.abc import x
>>> exp(x)
exp(x)
>>> exp(x).diff(x)
exp(x)
>>> exp(I*pi)
-1
Parameters
==========
arg : Expr
See Also
========
log
"""
def fdiff(self, argindex=1):
"""
Returns the first derivative of this function.
"""
if argindex == 1:
return self
else:
raise ArgumentIndexError(self, argindex)
def _eval_refine(self, assumptions):
from sympy.assumptions import ask, Q
arg = self.args[0]
if arg.is_Mul:
Ioo = I*S.Infinity
if arg in [Ioo, -Ioo]:
return S.NaN
coeff = arg.as_coefficient(pi*I)
if coeff:
if ask(Q.integer(2*coeff)):
if ask(Q.even(coeff)):
return S.One
elif ask(Q.odd(coeff)):
return S.NegativeOne
elif ask(Q.even(coeff + S.Half)):
return -I
elif ask(Q.odd(coeff + S.Half)):
return I
@classmethod
def eval(cls, arg):
from sympy.calculus import AccumBounds
from sympy.matrices.matrices import MatrixBase
from sympy.sets.setexpr import SetExpr
from sympy.simplify.simplify import logcombine
if isinstance(arg, MatrixBase):
return arg.exp()
elif global_parameters.exp_is_pow:
return Pow(S.Exp1, arg)
elif arg.is_Number:
if arg is S.NaN:
return S.NaN
elif arg.is_zero:
return S.One
elif arg is S.One:
return S.Exp1
elif arg is S.Infinity:
return S.Infinity
elif arg is S.NegativeInfinity:
return S.Zero
elif arg is S.ComplexInfinity:
return S.NaN
elif isinstance(arg, log):
return arg.args[0]
elif isinstance(arg, AccumBounds):
return AccumBounds(exp(arg.min), exp(arg.max))
elif isinstance(arg, SetExpr):
return arg._eval_func(cls)
elif arg.is_Mul:
coeff = arg.as_coefficient(pi*I)
if coeff:
if (2*coeff).is_integer:
if coeff.is_even:
return S.One
elif coeff.is_odd:
return S.NegativeOne
elif (coeff + S.Half).is_even:
return -I
elif (coeff + S.Half).is_odd:
return I
elif coeff.is_Rational:
ncoeff = coeff % 2 # restrict to [0, 2pi)
if ncoeff > 1: # restrict to (-pi, pi]
ncoeff -= 2
if ncoeff != coeff:
return cls(ncoeff*pi*I)
# Warning: code in risch.py will be very sensitive to changes
# in this (see DifferentialExtension).
# look for a single log factor
coeff, terms = arg.as_coeff_Mul()
# but it can't be multiplied by oo
if coeff in [S.NegativeInfinity, S.Infinity]:
if terms.is_number:
if coeff is S.NegativeInfinity:
terms = -terms
if re(terms).is_zero and terms is not S.Zero:
return S.NaN
if re(terms).is_positive and im(terms) is not S.Zero:
return S.ComplexInfinity
if re(terms).is_negative:
return S.Zero
return None
coeffs, log_term = [coeff], None
for term in Mul.make_args(terms):
term_ = logcombine(term)
if isinstance(term_, log):
if log_term is None:
log_term = term_.args[0]
else:
return None
elif term.is_comparable:
coeffs.append(term)
else:
return None
return log_term**Mul(*coeffs) if log_term else None
elif arg.is_Add:
out = []
add = []
argchanged = False
for a in arg.args:
if a is S.One:
add.append(a)
continue
newa = cls(a)
if isinstance(newa, cls):
if newa.args[0] != a:
add.append(newa.args[0])
argchanged = True
else:
add.append(a)
else:
out.append(newa)
if out or argchanged:
return Mul(*out)*cls(Add(*add), evaluate=False)
if arg.is_zero:
return S.One
@property
def base(self):
"""
Returns the base of the exponential function.
"""
return S.Exp1
@staticmethod
@cacheit
def taylor_term(n, x, *previous_terms):
"""
Calculates the next term in the Taylor series expansion.
"""
if n < 0:
return S.Zero
if n == 0:
return S.One
x = sympify(x)
if previous_terms:
p = previous_terms[-1]
if p is not None:
return p * x / n
return x**n/factorial(n)
def as_real_imag(self, deep=True, **hints):
"""
Returns this function as a 2-tuple representing a complex number.
Examples
========
>>> from sympy import exp, I
>>> from sympy.abc import x
>>> exp(x).as_real_imag()
(exp(re(x))*cos(im(x)), exp(re(x))*sin(im(x)))
>>> exp(1).as_real_imag()
(E, 0)
>>> exp(I).as_real_imag()
(cos(1), sin(1))
>>> exp(1+I).as_real_imag()
(E*cos(1), E*sin(1))
See Also
========
sympy.functions.elementary.complexes.re
sympy.functions.elementary.complexes.im
"""
from sympy.functions.elementary.trigonometric import cos, sin
re, im = self.args[0].as_real_imag()
if deep:
re = re.expand(deep, **hints)
im = im.expand(deep, **hints)
cos, sin = cos(im), sin(im)
return (exp(re)*cos, exp(re)*sin)
def _eval_subs(self, old, new):
# keep processing of power-like args centralized in Pow
if old.is_Pow: # handle (exp(3*log(x))).subs(x**2, z) -> z**(3/2)
old = exp(old.exp*log(old.base))
elif old is S.Exp1 and new.is_Function:
old = exp
if isinstance(old, exp) or old is S.Exp1:
f = lambda a: Pow(*a.as_base_exp(), evaluate=False) if (
a.is_Pow or isinstance(a, exp)) else a
return Pow._eval_subs(f(self), f(old), new)
if old is exp and not new.is_Function:
return new**self.exp._subs(old, new)
return Function._eval_subs(self, old, new)
def _eval_is_extended_real(self):
if self.args[0].is_extended_real:
return True
elif self.args[0].is_imaginary:
arg2 = -S(2) * I * self.args[0] / pi
return arg2.is_even
def _eval_is_complex(self):
def complex_extended_negative(arg):
yield arg.is_complex
yield arg.is_extended_negative
return fuzzy_or(complex_extended_negative(self.args[0]))
def _eval_is_algebraic(self):
if (self.exp / pi / I).is_rational:
return True
if fuzzy_not(self.exp.is_zero):
if self.exp.is_algebraic:
return False
elif (self.exp / pi).is_rational:
return False
def _eval_is_extended_positive(self):
if self.exp.is_extended_real:
return self.args[0] is not S.NegativeInfinity
elif self.exp.is_imaginary:
arg2 = -I * self.args[0] / pi
return arg2.is_even
def _eval_nseries(self, x, n, logx, cdir=0):
# NOTE Please see the comment at the beginning of this file, labelled
# IMPORTANT.
from sympy.functions.elementary.complexes import sign
from sympy.functions.elementary.integers import ceiling
from sympy.series.limits import limit
from sympy.series.order import Order
from sympy.simplify.powsimp import powsimp
arg = self.exp
arg_series = arg._eval_nseries(x, n=n, logx=logx)
if arg_series.is_Order:
return 1 + arg_series
arg0 = limit(arg_series.removeO(), x, 0)
if arg0 is S.NegativeInfinity:
return Order(x**n, x)
if arg0 is S.Infinity:
return self
# checking for indecisiveness/ sign terms in arg0
if any(isinstance(arg, (sign, ImaginaryUnit)) for arg in arg0.args):
return self
t = Dummy("t")
nterms = n
try:
cf = Order(arg.as_leading_term(x, logx=logx), x).getn()
except (NotImplementedError, PoleError):
cf = 0
if cf and cf > 0:
nterms = ceiling(n/cf)
exp_series = exp(t)._taylor(t, nterms)
r = exp(arg0)*exp_series.subs(t, arg_series - arg0)
rep = {logx: log(x)} if logx is not None else {}
if r.subs(rep) == self:
return r
if cf and cf > 1:
r += Order((arg_series - arg0)**n, x)/x**((cf-1)*n)
else:
r += Order((arg_series - arg0)**n, x)
r = r.expand()
r = powsimp(r, deep=True, combine='exp')
# powsimp may introduce unexpanded (-1)**Rational; see PR #17201
simplerat = lambda x: x.is_Rational and x.q in [3, 4, 6]
w = Wild('w', properties=[simplerat])
r = r.replace(S.NegativeOne**w, expand_complex(S.NegativeOne**w))
return r
def _taylor(self, x, n):
l = []
g = None
for i in range(n):
g = self.taylor_term(i, self.args[0], g)
g = g.nseries(x, n=n)
l.append(g.removeO())
return Add(*l)
def _eval_as_leading_term(self, x, logx=None, cdir=0):
from sympy.calculus.util import AccumBounds
arg = self.args[0].cancel().as_leading_term(x, logx=logx)
arg0 = arg.subs(x, 0)
if arg is S.NaN:
return S.NaN
if isinstance(arg0, AccumBounds):
# This check addresses a corner case involving AccumBounds.
# if isinstance(arg, AccumBounds) is True, then arg0 can either be 0,
# AccumBounds(-oo, 0) or AccumBounds(-oo, oo).
# Check out function: test_issue_18473() in test_exponential.py and
# test_limits.py for more information.
if re(cdir) < S.Zero:
return exp(-arg0)
return exp(arg0)
if arg0 is S.NaN:
arg0 = arg.limit(x, 0)
if arg0.is_infinite is False:
return exp(arg0)
raise PoleError("Cannot expand %s around 0" % (self))
def _eval_rewrite_as_sin(self, arg, **kwargs):
from sympy.functions.elementary.trigonometric import sin
return sin(I*arg + pi/2) - I*sin(I*arg)
def _eval_rewrite_as_cos(self, arg, **kwargs):
from sympy.functions.elementary.trigonometric import cos
return cos(I*arg) + I*cos(I*arg + pi/2)
def _eval_rewrite_as_tanh(self, arg, **kwargs):
from sympy.functions.elementary.hyperbolic import tanh
return (1 + tanh(arg/2))/(1 - tanh(arg/2))
def _eval_rewrite_as_sqrt(self, arg, **kwargs):
from sympy.functions.elementary.trigonometric import sin, cos
if arg.is_Mul:
coeff = arg.coeff(pi*I)
if coeff and coeff.is_number:
cosine, sine = cos(pi*coeff), sin(pi*coeff)
if not isinstance(cosine, cos) and not isinstance (sine, sin):
return cosine + I*sine
def _eval_rewrite_as_Pow(self, arg, **kwargs):
if arg.is_Mul:
logs = [a for a in arg.args if isinstance(a, log) and len(a.args) == 1]
if logs:
return Pow(logs[0].args[0], arg.coeff(logs[0]))
def match_real_imag(expr):
r"""
Try to match expr with $a + Ib$ for real $a$ and $b$.
``match_real_imag`` returns a tuple containing the real and imaginary
parts of expr or ``(None, None)`` if direct matching is not possible. Contrary
to :func:`~.re()`, :func:`~.im()``, and ``as_real_imag()``, this helper will not force things
by returning expressions themselves containing ``re()`` or ``im()`` and it
does not expand its argument either.
"""
r_, i_ = expr.as_independent(I, as_Add=True)
if i_ == 0 and r_.is_real:
return (r_, i_)
i_ = i_.as_coefficient(I)
if i_ and i_.is_real and r_.is_real:
return (r_, i_)
else:
return (None, None) # simpler to check for than None
class log(Function):
r"""
The natural logarithm function `\ln(x)` or `\log(x)`.
Explanation
===========
Logarithms are taken with the natural base, `e`. To get
a logarithm of a different base ``b``, use ``log(x, b)``,
which is essentially short-hand for ``log(x)/log(b)``.
``log`` represents the principal branch of the natural
logarithm. As such it has a branch cut along the negative
real axis and returns values having a complex argument in
`(-\pi, \pi]`.
Examples
========
>>> from sympy import log, sqrt, S, I
>>> log(8, 2)
3
>>> log(S(8)/3, 2)
-log(3)/log(2) + 3
>>> log(-1 + I*sqrt(3))
log(2) + 2*I*pi/3
See Also
========
exp
"""
args: tTuple[Expr]
_singularities = (S.Zero, S.ComplexInfinity)
def fdiff(self, argindex=1):
"""
Returns the first derivative of the function.
"""
if argindex == 1:
return 1/self.args[0]
else:
raise ArgumentIndexError(self, argindex)
def inverse(self, argindex=1):
r"""
Returns `e^x`, the inverse function of `\log(x)`.
"""
return exp
@classmethod
def eval(cls, arg, base=None):
from sympy.calculus import AccumBounds
from sympy.sets.setexpr import SetExpr
arg = sympify(arg)
if base is not None:
base = sympify(base)
if base == 1:
if arg == 1:
return S.NaN
else:
return S.ComplexInfinity
try:
# handle extraction of powers of the base now
# or else expand_log in Mul would have to handle this
n = multiplicity(base, arg)
if n:
return n + log(arg / base**n) / log(base)
else:
return log(arg)/log(base)
except ValueError:
pass
if base is not S.Exp1:
return cls(arg)/cls(base)
else:
return cls(arg)
if arg.is_Number:
if arg.is_zero:
return S.ComplexInfinity
elif arg is S.One:
return S.Zero
elif arg is S.Infinity:
return S.Infinity
elif arg is S.NegativeInfinity:
return S.Infinity
elif arg is S.NaN:
return S.NaN
elif arg.is_Rational and arg.p == 1:
return -cls(arg.q)
if arg.is_Pow and arg.base is S.Exp1 and arg.exp.is_extended_real:
return arg.exp
if isinstance(arg, exp) and arg.exp.is_extended_real:
return arg.exp
elif isinstance(arg, exp) and arg.exp.is_number:
r_, i_ = match_real_imag(arg.exp)
if i_ and i_.is_comparable:
i_ %= 2*pi
if i_ > pi:
i_ -= 2*pi
return r_ + expand_mul(i_ * I, deep=False)
elif isinstance(arg, exp_polar):
return unpolarify(arg.exp)
elif isinstance(arg, AccumBounds):
if arg.min.is_positive:
return AccumBounds(log(arg.min), log(arg.max))
elif arg.min.is_zero:
return AccumBounds(S.NegativeInfinity, log(arg.max))
else:
return S.NaN
elif isinstance(arg, SetExpr):
return arg._eval_func(cls)
if arg.is_number:
if arg.is_negative:
return pi * I + cls(-arg)
elif arg is S.ComplexInfinity:
return S.ComplexInfinity
elif arg is S.Exp1:
return S.One
if arg.is_zero:
return S.ComplexInfinity
# don't autoexpand Pow or Mul (see the issue 3351):
if not arg.is_Add:
coeff = arg.as_coefficient(I)
if coeff is not None:
if coeff is S.Infinity:
return S.Infinity
elif coeff is S.NegativeInfinity:
return S.Infinity
elif coeff.is_Rational:
if coeff.is_nonnegative:
return pi * I * S.Half + cls(coeff)
else:
return -pi * I * S.Half + cls(-coeff)
if arg.is_number and arg.is_algebraic:
# Match arg = coeff*(r_ + i_*I) with coeff>0, r_ and i_ real.
coeff, arg_ = arg.as_independent(I, as_Add=False)
if coeff.is_negative:
coeff *= -1
arg_ *= -1
arg_ = expand_mul(arg_, deep=False)
r_, i_ = arg_.as_independent(I, as_Add=True)
i_ = i_.as_coefficient(I)
if coeff.is_real and i_ and i_.is_real and r_.is_real:
if r_.is_zero:
if i_.is_positive:
return pi * I * S.Half + cls(coeff * i_)
elif i_.is_negative:
return -pi * I * S.Half + cls(coeff * -i_)
else:
from sympy.simplify import ratsimp
# Check for arguments involving rational multiples of pi
t = (i_/r_).cancel()
t1 = (-t).cancel()
atan_table = _log_atan_table()
if t in atan_table:
modulus = ratsimp(coeff * Abs(arg_))
if r_.is_positive:
return cls(modulus) + I * atan_table[t]
else:
return cls(modulus) + I * (atan_table[t] - pi)
elif t1 in atan_table:
modulus = ratsimp(coeff * Abs(arg_))
if r_.is_positive:
return cls(modulus) + I * (-atan_table[t1])
else:
return cls(modulus) + I * (pi - atan_table[t1])
def as_base_exp(self):
"""
Returns this function in the form (base, exponent).
"""
return self, S.One
@staticmethod
@cacheit
def taylor_term(n, x, *previous_terms): # of log(1+x)
r"""
Returns the next term in the Taylor series expansion of `\log(1+x)`.
"""
from sympy.simplify.powsimp import powsimp
if n < 0:
return S.Zero
x = sympify(x)
if n == 0:
return x
if previous_terms:
p = previous_terms[-1]
if p is not None:
return powsimp((-n) * p * x / (n + 1), deep=True, combine='exp')
return (1 - 2*(n % 2)) * x**(n + 1)/(n + 1)
def _eval_expand_log(self, deep=True, **hints):
from sympy.concrete import Sum, Product
force = hints.get('force', False)
factor = hints.get('factor', False)
if (len(self.args) == 2):
return expand_log(self.func(*self.args), deep=deep, force=force)
arg = self.args[0]
if arg.is_Integer:
# remove perfect powers
p = perfect_power(arg)
logarg = None
coeff = 1
if p is not False:
arg, coeff = p
logarg = self.func(arg)
# expand as product of its prime factors if factor=True
if factor:
p = factorint(arg)
if arg not in p.keys():
logarg = sum(n*log(val) for val, n in p.items())
if logarg is not None:
return coeff*logarg
elif arg.is_Rational:
return log(arg.p) - log(arg.q)
elif arg.is_Mul:
expr = []
nonpos = []
for x in arg.args:
if force or x.is_positive or x.is_polar:
a = self.func(x)
if isinstance(a, log):
expr.append(self.func(x)._eval_expand_log(**hints))
else:
expr.append(a)
elif x.is_negative:
a = self.func(-x)
expr.append(a)
nonpos.append(S.NegativeOne)
else:
nonpos.append(x)
return Add(*expr) + log(Mul(*nonpos))
elif arg.is_Pow or isinstance(arg, exp):
if force or (arg.exp.is_extended_real and (arg.base.is_positive or ((arg.exp+1)
.is_positive and (arg.exp-1).is_nonpositive))) or arg.base.is_polar:
b = arg.base
e = arg.exp
a = self.func(b)
if isinstance(a, log):
return unpolarify(e) * a._eval_expand_log(**hints)
else:
return unpolarify(e) * a
elif isinstance(arg, Product):
if force or arg.function.is_positive:
return Sum(log(arg.function), *arg.limits)
return self.func(arg)
def _eval_simplify(self, **kwargs):
from sympy.simplify.simplify import expand_log, simplify, inversecombine
if len(self.args) == 2: # it's unevaluated
return simplify(self.func(*self.args), **kwargs)
expr = self.func(simplify(self.args[0], **kwargs))
if kwargs['inverse']:
expr = inversecombine(expr)
expr = expand_log(expr, deep=True)
return min([expr, self], key=kwargs['measure'])
def as_real_imag(self, deep=True, **hints):
"""
Returns this function as a complex coordinate.
Examples
========
>>> from sympy import I, log
>>> from sympy.abc import x
>>> log(x).as_real_imag()
(log(Abs(x)), arg(x))
>>> log(I).as_real_imag()
(0, pi/2)
>>> log(1 + I).as_real_imag()
(log(sqrt(2)), pi/4)
>>> log(I*x).as_real_imag()
(log(Abs(x)), arg(I*x))
"""
sarg = self.args[0]
if deep:
sarg = self.args[0].expand(deep, **hints)
sarg_abs = Abs(sarg)
if sarg_abs == sarg:
return self, S.Zero
sarg_arg = arg(sarg)
if hints.get('log', False): # Expand the log
hints['complex'] = False
return (log(sarg_abs).expand(deep, **hints), sarg_arg)
else:
return log(sarg_abs), sarg_arg
def _eval_is_rational(self):
s = self.func(*self.args)
if s.func == self.func:
if (self.args[0] - 1).is_zero:
return True
if s.args[0].is_rational and fuzzy_not((self.args[0] - 1).is_zero):
return False
else:
return s.is_rational
def _eval_is_algebraic(self):
s = self.func(*self.args)
if s.func == self.func:
if (self.args[0] - 1).is_zero:
return True
elif fuzzy_not((self.args[0] - 1).is_zero):
if self.args[0].is_algebraic:
return False
else:
return s.is_algebraic
def _eval_is_extended_real(self):
return self.args[0].is_extended_positive
def _eval_is_complex(self):
z = self.args[0]
return fuzzy_and([z.is_complex, fuzzy_not(z.is_zero)])
def _eval_is_finite(self):
arg = self.args[0]
if arg.is_zero:
return False
return arg.is_finite
def _eval_is_extended_positive(self):
return (self.args[0] - 1).is_extended_positive
def _eval_is_zero(self):
return (self.args[0] - 1).is_zero
def _eval_is_extended_nonnegative(self):
return (self.args[0] - 1).is_extended_nonnegative
def _eval_nseries(self, x, n, logx, cdir=0):
# NOTE Please see the comment at the beginning of this file, labelled
# IMPORTANT.
from sympy.series.order import Order
from sympy.simplify.simplify import logcombine
from sympy.core.symbol import Dummy
if self.args[0] == x:
return log(x) if logx is None else logx
arg = self.args[0]
t = Dummy('t', positive=True)
if cdir == 0:
cdir = 1
z = arg.subs(x, cdir*t)
k, l = Wild("k"), Wild("l")
r = z.match(k*t**l)
if r is not None:
k, l = r[k], r[l]
if l != 0 and not l.has(t) and not k.has(t):
r = l*log(x) if logx is None else l*logx
r += log(k) - l*log(cdir) # XXX true regardless of assumptions?
return r
def coeff_exp(term, x):
coeff, exp = S.One, S.Zero
for factor in Mul.make_args(term):
if factor.has(x):
base, exp = factor.as_base_exp()
if base != x:
try:
return term.leadterm(x)
except ValueError:
return term, S.Zero
else:
coeff *= factor
return coeff, exp
# TODO new and probably slow
try:
a, b = z.leadterm(t, logx=logx, cdir=1)
except (ValueError, NotImplementedError, PoleError):
s = z._eval_nseries(t, n=n, logx=logx, cdir=1)
while s.is_Order:
n += 1
s = z._eval_nseries(t, n=n, logx=logx, cdir=1)
try:
a, b = s.removeO().leadterm(t, cdir=1)
except ValueError:
a, b = s.removeO().as_leading_term(t, cdir=1), S.Zero
p = (z/(a*t**b) - 1)._eval_nseries(t, n=n, logx=logx, cdir=1)
if p.has(exp):
p = logcombine(p)
if isinstance(p, Order):
n = p.getn()
_, d = coeff_exp(p, t)
logx = log(x) if logx is None else logx
if not d.is_positive:
res = log(a) - b*log(cdir) + b*logx
_res = res
logflags = {"deep": True, "log": True, "mul": False, "power_exp": False,
"power_base": False, "multinomial": False, "basic": False, "force": True,
"factor": False}
expr = self.expand(**logflags)
if (not a.could_extract_minus_sign() and
logx.could_extract_minus_sign()):
_res = _res.subs(-logx, -log(x)).expand(**logflags)
else:
_res = _res.subs(logx, log(x)).expand(**logflags)
if _res == expr:
return res
return res + Order(x**n, x)
def mul(d1, d2):
res = {}
for e1, e2 in product(d1, d2):
ex = e1 + e2
if ex < n:
res[ex] = res.get(ex, S.Zero) + d1[e1]*d2[e2]
return res
pterms = {}
for term in Add.make_args(p.removeO()):
co1, e1 = coeff_exp(term, t)
pterms[e1] = pterms.get(e1, S.Zero) + co1
k = S.One
terms = {}
pk = pterms
while k*d < n:
coeff = -S.NegativeOne**k/k
for ex in pk:
_ = terms.get(ex, S.Zero) + coeff*pk[ex]
terms[ex] = _.nsimplify()
pk = mul(pk, pterms)
k += S.One
res = log(a) - b*log(cdir) + b*logx
for ex in terms:
res += terms[ex]*t**(ex)
if a.is_negative and im(z) != 0:
from sympy.functions.special.delta_functions import Heaviside
for i, term in enumerate(z.lseries(t)):
if not term.is_real or i == 5:
break
if i < 5:
coeff, _ = term.as_coeff_exponent(t)
res += -2*I*pi*Heaviside(-im(coeff), 0)
res = res.subs(t, x/cdir)
return res + Order(x**n, x)
def _eval_as_leading_term(self, x, logx=None, cdir=0):
# NOTE
# Refer https://github.com/sympy/sympy/pull/23592 for more information
# on each of the following steps involved in this method.
arg0 = self.args[0].together()
# STEP 1
t = Dummy('t', positive=True)
if cdir == 0:
cdir = 1
z = arg0.subs(x, cdir*t)
# STEP 2
try:
c, e = z.leadterm(t, logx=logx, cdir=1)
except ValueError:
arg = arg0.as_leading_term(x, logx=logx, cdir=cdir)
return log(arg)
if c.has(t):
c = c.subs(t, x/cdir)
if e != 0:
raise PoleError("Cannot expand %s around 0" % (self))
return log(c)
# STEP 3
if c == S.One and e == S.Zero:
return (arg0 - S.One).as_leading_term(x, logx=logx)
# STEP 4
res = log(c) - e*log(cdir)
logx = log(x) if logx is None else logx
res += e*logx
# STEP 5
if c.is_negative and im(z) != 0:
from sympy.functions.special.delta_functions import Heaviside
for i, term in enumerate(z.lseries(t)):
if not term.is_real or i == 5:
break
if i < 5:
coeff, _ = term.as_coeff_exponent(t)
res += -2*I*pi*Heaviside(-im(coeff), 0)
return res
class LambertW(Function):
r"""
The Lambert W function $W(z)$ is defined as the inverse
function of $w \exp(w)$ [1]_.
Explanation
===========
In other words, the value of $W(z)$ is such that $z = W(z) \exp(W(z))$
for any complex number $z$. The Lambert W function is a multivalued
function with infinitely many branches $W_k(z)$, indexed by
$k \in \mathbb{Z}$. Each branch gives a different solution $w$
of the equation $z = w \exp(w)$.
The Lambert W function has two partially real branches: the
principal branch ($k = 0$) is real for real $z > -1/e$, and the
$k = -1$ branch is real for $-1/e < z < 0$. All branches except
$k = 0$ have a logarithmic singularity at $z = 0$.
Examples
========
>>> from sympy import LambertW
>>> LambertW(1.2)
0.635564016364870
>>> LambertW(1.2, -1).n()
-1.34747534407696 - 4.41624341514535*I
>>> LambertW(-1).is_real
False
References
==========
.. [1] https://en.wikipedia.org/wiki/Lambert_W_function
"""
_singularities = (-Pow(S.Exp1, -1, evaluate=False), S.ComplexInfinity)
@classmethod
def eval(cls, x, k=None):
if k == S.Zero:
return cls(x)
elif k is None:
k = S.Zero
if k.is_zero:
if x.is_zero:
return S.Zero
if x is S.Exp1:
return S.One
if x == -1/S.Exp1:
return S.NegativeOne
if x == -log(2)/2:
return -log(2)
if x == 2*log(2):
return log(2)
if x == -pi/2:
return I*pi/2
if x == exp(1 + S.Exp1):
return S.Exp1
if x is S.Infinity:
return S.Infinity
if x.is_zero:
return S.Zero
if fuzzy_not(k.is_zero):
if x.is_zero:
return S.NegativeInfinity
if k is S.NegativeOne:
if x == -pi/2:
return -I*pi/2
elif x == -1/S.Exp1:
return S.NegativeOne
elif x == -2*exp(-2):
return -Integer(2)
def fdiff(self, argindex=1):
"""
Return the first derivative of this function.
"""
x = self.args[0]
if len(self.args) == 1:
if argindex == 1:
return LambertW(x)/(x*(1 + LambertW(x)))
else:
k = self.args[1]
if argindex == 1:
return LambertW(x, k)/(x*(1 + LambertW(x, k)))
raise ArgumentIndexError(self, argindex)
def _eval_is_extended_real(self):
x = self.args[0]
if len(self.args) == 1:
k = S.Zero
else:
k = self.args[1]
if k.is_zero:
if (x + 1/S.Exp1).is_positive:
return True
elif (x + 1/S.Exp1).is_nonpositive:
return False
elif (k + 1).is_zero:
if x.is_negative and (x + 1/S.Exp1).is_positive:
return True
elif x.is_nonpositive or (x + 1/S.Exp1).is_nonnegative:
return False
elif fuzzy_not(k.is_zero) and fuzzy_not((k + 1).is_zero):
if x.is_extended_real:
return False
def _eval_is_finite(self):
return self.args[0].is_finite
def _eval_is_algebraic(self):
s = self.func(*self.args)
if s.func == self.func:
if fuzzy_not(self.args[0].is_zero) and self.args[0].is_algebraic:
return False
else:
return s.is_algebraic
def _eval_as_leading_term(self, x, logx=None, cdir=0):
if len(self.args) == 1:
arg = self.args[0]
arg0 = arg.subs(x, 0).cancel()
if not arg0.is_zero:
return self.func(arg0)
return arg.as_leading_term(x)
def _eval_nseries(self, x, n, logx, cdir=0):
if len(self.args) == 1:
from sympy.functions.elementary.integers import ceiling
from sympy.series.order import Order
arg = self.args[0].nseries(x, n=n, logx=logx)
lt = arg.as_leading_term(x, logx=logx)
lte = 1
if lt.is_Pow:
lte = lt.exp
if ceiling(n/lte) >= 1:
s = Add(*[(-S.One)**(k - 1)*Integer(k)**(k - 2)/
factorial(k - 1)*arg**k for k in range(1, ceiling(n/lte))])
s = expand_multinomial(s)
else:
s = S.Zero
return s + Order(x**n, x)
return super()._eval_nseries(x, n, logx)
def _eval_is_zero(self):
x = self.args[0]
if len(self.args) == 1:
return x.is_zero
else:
return fuzzy_and([x.is_zero, self.args[1].is_zero])
@cacheit
def _log_atan_table():
return {
# first quadrant only
sqrt(3): pi / 3,
1: pi / 4,
sqrt(5 - 2 * sqrt(5)): pi / 5,
sqrt(2) * sqrt(5 - sqrt(5)) / (1 + sqrt(5)): pi / 5,
sqrt(5 + 2 * sqrt(5)): pi * Rational(2, 5),
sqrt(2) * sqrt(sqrt(5) + 5) / (-1 + sqrt(5)): pi * Rational(2, 5),
sqrt(3) / 3: pi / 6,
sqrt(2) - 1: pi / 8,
sqrt(2 - sqrt(2)) / sqrt(sqrt(2) + 2): pi / 8,
sqrt(2) + 1: pi * Rational(3, 8),
sqrt(sqrt(2) + 2) / sqrt(2 - sqrt(2)): pi * Rational(3, 8),
sqrt(1 - 2 * sqrt(5) / 5): pi / 10,
(-sqrt(2) + sqrt(10)) / (2 * sqrt(sqrt(5) + 5)): pi / 10,
sqrt(1 + 2 * sqrt(5) / 5): pi * Rational(3, 10),
(sqrt(2) + sqrt(10)) / (2 * sqrt(5 - sqrt(5))): pi * Rational(3, 10),
2 - sqrt(3): pi / 12,
(-1 + sqrt(3)) / (1 + sqrt(3)): pi / 12,
2 + sqrt(3): pi * Rational(5, 12),
(1 + sqrt(3)) / (-1 + sqrt(3)): pi * Rational(5, 12)
}