Coverage for /usr/lib/python3/dist-packages/sympy/functions/elementary/exponential.py: 18%

830 statements  

« prev     ^ index     » next       coverage.py v7.9.1, created at 2025-06-14 15:55 +0200

1from itertools import product 

2from typing import Tuple as tTuple 

3 

4from sympy.core.add import Add 

5from sympy.core.cache import cacheit 

6from sympy.core.expr import Expr 

7from sympy.core.function import (Function, ArgumentIndexError, expand_log, 

8 expand_mul, FunctionClass, PoleError, expand_multinomial, expand_complex) 

9from sympy.core.logic import fuzzy_and, fuzzy_not, fuzzy_or 

10from sympy.core.mul import Mul 

11from sympy.core.numbers import Integer, Rational, pi, I, ImaginaryUnit 

12from sympy.core.parameters import global_parameters 

13from sympy.core.power import Pow 

14from sympy.core.singleton import S 

15from sympy.core.symbol import Wild, Dummy 

16from sympy.core.sympify import sympify 

17from sympy.functions.combinatorial.factorials import factorial 

18from sympy.functions.elementary.complexes import arg, unpolarify, im, re, Abs 

19from sympy.functions.elementary.miscellaneous import sqrt 

20from sympy.ntheory import multiplicity, perfect_power 

21from sympy.ntheory.factor_ import factorint 

22 

23# NOTE IMPORTANT 

24# The series expansion code in this file is an important part of the gruntz 

25# algorithm for determining limits. _eval_nseries has to return a generalized 

26# power series with coefficients in C(log(x), log). 

27# In more detail, the result of _eval_nseries(self, x, n) must be 

28# c_0*x**e_0 + ... (finitely many terms) 

29# where e_i are numbers (not necessarily integers) and c_i involve only 

30# numbers, the function log, and log(x). [This also means it must not contain 

31# log(x(1+p)), this *has* to be expanded to log(x)+log(1+p) if x.is_positive and 

32# p.is_positive.] 

33 

34 

35class ExpBase(Function): 

36 

37 unbranched = True 

38 _singularities = (S.ComplexInfinity,) 

39 

40 @property 

41 def kind(self): 

42 return self.exp.kind 

43 

44 def inverse(self, argindex=1): 

45 """ 

46 Returns the inverse function of ``exp(x)``. 

47 """ 

48 return log 

49 

50 def as_numer_denom(self): 

51 """ 

52 Returns this with a positive exponent as a 2-tuple (a fraction). 

53 

54 Examples 

55 ======== 

56 

57 >>> from sympy import exp 

58 >>> from sympy.abc import x 

59 >>> exp(-x).as_numer_denom() 

60 (1, exp(x)) 

61 >>> exp(x).as_numer_denom() 

62 (exp(x), 1) 

63 """ 

64 # this should be the same as Pow.as_numer_denom wrt 

65 # exponent handling 

66 exp = self.exp 

67 neg_exp = exp.is_negative 

68 if not neg_exp and not (-exp).is_negative: 

69 neg_exp = exp.could_extract_minus_sign() 

70 if neg_exp: 

71 return S.One, self.func(-exp) 

72 return self, S.One 

73 

74 @property 

75 def exp(self): 

76 """ 

77 Returns the exponent of the function. 

78 """ 

79 return self.args[0] 

80 

81 def as_base_exp(self): 

82 """ 

83 Returns the 2-tuple (base, exponent). 

84 """ 

85 return self.func(1), Mul(*self.args) 

86 

87 def _eval_adjoint(self): 

88 return self.func(self.exp.adjoint()) 

89 

90 def _eval_conjugate(self): 

91 return self.func(self.exp.conjugate()) 

92 

93 def _eval_transpose(self): 

94 return self.func(self.exp.transpose()) 

95 

96 def _eval_is_finite(self): 

97 arg = self.exp 

98 if arg.is_infinite: 

99 if arg.is_extended_negative: 

100 return True 

101 if arg.is_extended_positive: 

102 return False 

103 if arg.is_finite: 

104 return True 

105 

106 def _eval_is_rational(self): 

107 s = self.func(*self.args) 

108 if s.func == self.func: 

109 z = s.exp.is_zero 

110 if z: 

111 return True 

112 elif s.exp.is_rational and fuzzy_not(z): 

113 return False 

114 else: 

115 return s.is_rational 

116 

117 def _eval_is_zero(self): 

118 return self.exp is S.NegativeInfinity 

119 

120 def _eval_power(self, other): 

121 """exp(arg)**e -> exp(arg*e) if assumptions allow it. 

122 """ 

123 b, e = self.as_base_exp() 

124 return Pow._eval_power(Pow(b, e, evaluate=False), other) 

125 

126 def _eval_expand_power_exp(self, **hints): 

127 from sympy.concrete.products import Product 

128 from sympy.concrete.summations import Sum 

129 arg = self.args[0] 

130 if arg.is_Add and arg.is_commutative: 

131 return Mul.fromiter(self.func(x) for x in arg.args) 

132 elif isinstance(arg, Sum) and arg.is_commutative: 

133 return Product(self.func(arg.function), *arg.limits) 

134 return self.func(arg) 

135 

136 

137class exp_polar(ExpBase): 

138 r""" 

139 Represent a *polar number* (see g-function Sphinx documentation). 

140 

141 Explanation 

142 =========== 

143 

144 ``exp_polar`` represents the function 

145 `Exp: \mathbb{C} \rightarrow \mathcal{S}`, sending the complex number 

146 `z = a + bi` to the polar number `r = exp(a), \theta = b`. It is one of 

147 the main functions to construct polar numbers. 

148 

149 Examples 

150 ======== 

151 

152 >>> from sympy import exp_polar, pi, I, exp 

153 

154 The main difference is that polar numbers do not "wrap around" at `2 \pi`: 

155 

156 >>> exp(2*pi*I) 

157 1 

158 >>> exp_polar(2*pi*I) 

159 exp_polar(2*I*pi) 

160 

161 apart from that they behave mostly like classical complex numbers: 

162 

163 >>> exp_polar(2)*exp_polar(3) 

164 exp_polar(5) 

165 

166 See Also 

167 ======== 

168 

169 sympy.simplify.powsimp.powsimp 

170 polar_lift 

171 periodic_argument 

172 principal_branch 

173 """ 

174 

175 is_polar = True 

176 is_comparable = False # cannot be evalf'd 

177 

178 def _eval_Abs(self): # Abs is never a polar number 

179 return exp(re(self.args[0])) 

180 

181 def _eval_evalf(self, prec): 

182 """ Careful! any evalf of polar numbers is flaky """ 

183 i = im(self.args[0]) 

184 try: 

185 bad = (i <= -pi or i > pi) 

186 except TypeError: 

187 bad = True 

188 if bad: 

189 return self # cannot evalf for this argument 

190 res = exp(self.args[0])._eval_evalf(prec) 

191 if i > 0 and im(res) < 0: 

192 # i ~ pi, but exp(I*i) evaluated to argument slightly bigger than pi 

193 return re(res) 

194 return res 

195 

196 def _eval_power(self, other): 

197 return self.func(self.args[0]*other) 

198 

199 def _eval_is_extended_real(self): 

200 if self.args[0].is_extended_real: 

201 return True 

202 

203 def as_base_exp(self): 

204 # XXX exp_polar(0) is special! 

205 if self.args[0] == 0: 

206 return self, S.One 

207 return ExpBase.as_base_exp(self) 

208 

209 

210class ExpMeta(FunctionClass): 

211 def __instancecheck__(cls, instance): 

212 if exp in instance.__class__.__mro__: 

213 return True 

214 return isinstance(instance, Pow) and instance.base is S.Exp1 

215 

216 

217class exp(ExpBase, metaclass=ExpMeta): 

218 """ 

219 The exponential function, :math:`e^x`. 

220 

221 Examples 

222 ======== 

223 

224 >>> from sympy import exp, I, pi 

225 >>> from sympy.abc import x 

226 >>> exp(x) 

227 exp(x) 

228 >>> exp(x).diff(x) 

229 exp(x) 

230 >>> exp(I*pi) 

231 -1 

232 

233 Parameters 

234 ========== 

235 

236 arg : Expr 

237 

238 See Also 

239 ======== 

240 

241 log 

242 """ 

243 

244 def fdiff(self, argindex=1): 

245 """ 

246 Returns the first derivative of this function. 

247 """ 

248 if argindex == 1: 

249 return self 

250 else: 

251 raise ArgumentIndexError(self, argindex) 

252 

253 def _eval_refine(self, assumptions): 

254 from sympy.assumptions import ask, Q 

255 arg = self.args[0] 

256 if arg.is_Mul: 

257 Ioo = I*S.Infinity 

258 if arg in [Ioo, -Ioo]: 

259 return S.NaN 

260 

261 coeff = arg.as_coefficient(pi*I) 

262 if coeff: 

263 if ask(Q.integer(2*coeff)): 

264 if ask(Q.even(coeff)): 

265 return S.One 

266 elif ask(Q.odd(coeff)): 

267 return S.NegativeOne 

268 elif ask(Q.even(coeff + S.Half)): 

269 return -I 

270 elif ask(Q.odd(coeff + S.Half)): 

271 return I 

272 

273 @classmethod 

274 def eval(cls, arg): 

275 from sympy.calculus import AccumBounds 

276 from sympy.matrices.matrices import MatrixBase 

277 from sympy.sets.setexpr import SetExpr 

278 from sympy.simplify.simplify import logcombine 

279 if isinstance(arg, MatrixBase): 

280 return arg.exp() 

281 elif global_parameters.exp_is_pow: 

282 return Pow(S.Exp1, arg) 

283 elif arg.is_Number: 

284 if arg is S.NaN: 

285 return S.NaN 

286 elif arg.is_zero: 

287 return S.One 

288 elif arg is S.One: 

289 return S.Exp1 

290 elif arg is S.Infinity: 

291 return S.Infinity 

292 elif arg is S.NegativeInfinity: 

293 return S.Zero 

294 elif arg is S.ComplexInfinity: 

295 return S.NaN 

296 elif isinstance(arg, log): 

297 return arg.args[0] 

298 elif isinstance(arg, AccumBounds): 

299 return AccumBounds(exp(arg.min), exp(arg.max)) 

300 elif isinstance(arg, SetExpr): 

301 return arg._eval_func(cls) 

302 elif arg.is_Mul: 

303 coeff = arg.as_coefficient(pi*I) 

304 if coeff: 

305 if (2*coeff).is_integer: 

306 if coeff.is_even: 

307 return S.One 

308 elif coeff.is_odd: 

309 return S.NegativeOne 

310 elif (coeff + S.Half).is_even: 

311 return -I 

312 elif (coeff + S.Half).is_odd: 

313 return I 

314 elif coeff.is_Rational: 

315 ncoeff = coeff % 2 # restrict to [0, 2pi) 

316 if ncoeff > 1: # restrict to (-pi, pi] 

317 ncoeff -= 2 

318 if ncoeff != coeff: 

319 return cls(ncoeff*pi*I) 

320 

321 # Warning: code in risch.py will be very sensitive to changes 

322 # in this (see DifferentialExtension). 

323 

324 # look for a single log factor 

325 

326 coeff, terms = arg.as_coeff_Mul() 

327 

328 # but it can't be multiplied by oo 

329 if coeff in [S.NegativeInfinity, S.Infinity]: 

330 if terms.is_number: 

331 if coeff is S.NegativeInfinity: 

332 terms = -terms 

333 if re(terms).is_zero and terms is not S.Zero: 

334 return S.NaN 

335 if re(terms).is_positive and im(terms) is not S.Zero: 

336 return S.ComplexInfinity 

337 if re(terms).is_negative: 

338 return S.Zero 

339 return None 

340 

341 coeffs, log_term = [coeff], None 

342 for term in Mul.make_args(terms): 

343 term_ = logcombine(term) 

344 if isinstance(term_, log): 

345 if log_term is None: 

346 log_term = term_.args[0] 

347 else: 

348 return None 

349 elif term.is_comparable: 

350 coeffs.append(term) 

351 else: 

352 return None 

353 

354 return log_term**Mul(*coeffs) if log_term else None 

355 

356 elif arg.is_Add: 

357 out = [] 

358 add = [] 

359 argchanged = False 

360 for a in arg.args: 

361 if a is S.One: 

362 add.append(a) 

363 continue 

364 newa = cls(a) 

365 if isinstance(newa, cls): 

366 if newa.args[0] != a: 

367 add.append(newa.args[0]) 

368 argchanged = True 

369 else: 

370 add.append(a) 

371 else: 

372 out.append(newa) 

373 if out or argchanged: 

374 return Mul(*out)*cls(Add(*add), evaluate=False) 

375 

376 if arg.is_zero: 

377 return S.One 

378 

379 @property 

380 def base(self): 

381 """ 

382 Returns the base of the exponential function. 

383 """ 

384 return S.Exp1 

385 

386 @staticmethod 

387 @cacheit 

388 def taylor_term(n, x, *previous_terms): 

389 """ 

390 Calculates the next term in the Taylor series expansion. 

391 """ 

392 if n < 0: 

393 return S.Zero 

394 if n == 0: 

395 return S.One 

396 x = sympify(x) 

397 if previous_terms: 

398 p = previous_terms[-1] 

399 if p is not None: 

400 return p * x / n 

401 return x**n/factorial(n) 

402 

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

404 """ 

405 Returns this function as a 2-tuple representing a complex number. 

406 

407 Examples 

408 ======== 

409 

410 >>> from sympy import exp, I 

411 >>> from sympy.abc import x 

412 >>> exp(x).as_real_imag() 

413 (exp(re(x))*cos(im(x)), exp(re(x))*sin(im(x))) 

414 >>> exp(1).as_real_imag() 

415 (E, 0) 

416 >>> exp(I).as_real_imag() 

417 (cos(1), sin(1)) 

418 >>> exp(1+I).as_real_imag() 

419 (E*cos(1), E*sin(1)) 

420 

421 See Also 

422 ======== 

423 

424 sympy.functions.elementary.complexes.re 

425 sympy.functions.elementary.complexes.im 

426 """ 

427 from sympy.functions.elementary.trigonometric import cos, sin 

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

429 if deep: 

430 re = re.expand(deep, **hints) 

431 im = im.expand(deep, **hints) 

432 cos, sin = cos(im), sin(im) 

433 return (exp(re)*cos, exp(re)*sin) 

434 

435 def _eval_subs(self, old, new): 

436 # keep processing of power-like args centralized in Pow 

437 if old.is_Pow: # handle (exp(3*log(x))).subs(x**2, z) -> z**(3/2) 

438 old = exp(old.exp*log(old.base)) 

439 elif old is S.Exp1 and new.is_Function: 

440 old = exp 

441 if isinstance(old, exp) or old is S.Exp1: 

442 f = lambda a: Pow(*a.as_base_exp(), evaluate=False) if ( 

443 a.is_Pow or isinstance(a, exp)) else a 

444 return Pow._eval_subs(f(self), f(old), new) 

445 

446 if old is exp and not new.is_Function: 

447 return new**self.exp._subs(old, new) 

448 return Function._eval_subs(self, old, new) 

449 

450 def _eval_is_extended_real(self): 

451 if self.args[0].is_extended_real: 

452 return True 

453 elif self.args[0].is_imaginary: 

454 arg2 = -S(2) * I * self.args[0] / pi 

455 return arg2.is_even 

456 

457 def _eval_is_complex(self): 

458 def complex_extended_negative(arg): 

459 yield arg.is_complex 

460 yield arg.is_extended_negative 

461 return fuzzy_or(complex_extended_negative(self.args[0])) 

462 

463 def _eval_is_algebraic(self): 

464 if (self.exp / pi / I).is_rational: 

465 return True 

466 if fuzzy_not(self.exp.is_zero): 

467 if self.exp.is_algebraic: 

468 return False 

469 elif (self.exp / pi).is_rational: 

470 return False 

471 

472 def _eval_is_extended_positive(self): 

473 if self.exp.is_extended_real: 

474 return self.args[0] is not S.NegativeInfinity 

475 elif self.exp.is_imaginary: 

476 arg2 = -I * self.args[0] / pi 

477 return arg2.is_even 

478 

479 def _eval_nseries(self, x, n, logx, cdir=0): 

480 # NOTE Please see the comment at the beginning of this file, labelled 

481 # IMPORTANT. 

482 from sympy.functions.elementary.complexes import sign 

483 from sympy.functions.elementary.integers import ceiling 

484 from sympy.series.limits import limit 

485 from sympy.series.order import Order 

486 from sympy.simplify.powsimp import powsimp 

487 arg = self.exp 

488 arg_series = arg._eval_nseries(x, n=n, logx=logx) 

489 if arg_series.is_Order: 

490 return 1 + arg_series 

491 arg0 = limit(arg_series.removeO(), x, 0) 

492 if arg0 is S.NegativeInfinity: 

493 return Order(x**n, x) 

494 if arg0 is S.Infinity: 

495 return self 

496 # checking for indecisiveness/ sign terms in arg0 

497 if any(isinstance(arg, (sign, ImaginaryUnit)) for arg in arg0.args): 

498 return self 

499 t = Dummy("t") 

500 nterms = n 

501 try: 

502 cf = Order(arg.as_leading_term(x, logx=logx), x).getn() 

503 except (NotImplementedError, PoleError): 

504 cf = 0 

505 if cf and cf > 0: 

506 nterms = ceiling(n/cf) 

507 exp_series = exp(t)._taylor(t, nterms) 

508 r = exp(arg0)*exp_series.subs(t, arg_series - arg0) 

509 rep = {logx: log(x)} if logx is not None else {} 

510 if r.subs(rep) == self: 

511 return r 

512 if cf and cf > 1: 

513 r += Order((arg_series - arg0)**n, x)/x**((cf-1)*n) 

514 else: 

515 r += Order((arg_series - arg0)**n, x) 

516 r = r.expand() 

517 r = powsimp(r, deep=True, combine='exp') 

518 # powsimp may introduce unexpanded (-1)**Rational; see PR #17201 

519 simplerat = lambda x: x.is_Rational and x.q in [3, 4, 6] 

520 w = Wild('w', properties=[simplerat]) 

521 r = r.replace(S.NegativeOne**w, expand_complex(S.NegativeOne**w)) 

522 return r 

523 

524 def _taylor(self, x, n): 

525 l = [] 

526 g = None 

527 for i in range(n): 

528 g = self.taylor_term(i, self.args[0], g) 

529 g = g.nseries(x, n=n) 

530 l.append(g.removeO()) 

531 return Add(*l) 

532 

533 def _eval_as_leading_term(self, x, logx=None, cdir=0): 

534 from sympy.calculus.util import AccumBounds 

535 arg = self.args[0].cancel().as_leading_term(x, logx=logx) 

536 arg0 = arg.subs(x, 0) 

537 if arg is S.NaN: 

538 return S.NaN 

539 if isinstance(arg0, AccumBounds): 

540 # This check addresses a corner case involving AccumBounds. 

541 # if isinstance(arg, AccumBounds) is True, then arg0 can either be 0, 

542 # AccumBounds(-oo, 0) or AccumBounds(-oo, oo). 

543 # Check out function: test_issue_18473() in test_exponential.py and 

544 # test_limits.py for more information. 

545 if re(cdir) < S.Zero: 

546 return exp(-arg0) 

547 return exp(arg0) 

548 if arg0 is S.NaN: 

549 arg0 = arg.limit(x, 0) 

550 if arg0.is_infinite is False: 

551 return exp(arg0) 

552 raise PoleError("Cannot expand %s around 0" % (self)) 

553 

554 def _eval_rewrite_as_sin(self, arg, **kwargs): 

555 from sympy.functions.elementary.trigonometric import sin 

556 return sin(I*arg + pi/2) - I*sin(I*arg) 

557 

558 def _eval_rewrite_as_cos(self, arg, **kwargs): 

559 from sympy.functions.elementary.trigonometric import cos 

560 return cos(I*arg) + I*cos(I*arg + pi/2) 

561 

562 def _eval_rewrite_as_tanh(self, arg, **kwargs): 

563 from sympy.functions.elementary.hyperbolic import tanh 

564 return (1 + tanh(arg/2))/(1 - tanh(arg/2)) 

565 

566 def _eval_rewrite_as_sqrt(self, arg, **kwargs): 

567 from sympy.functions.elementary.trigonometric import sin, cos 

568 if arg.is_Mul: 

569 coeff = arg.coeff(pi*I) 

570 if coeff and coeff.is_number: 

571 cosine, sine = cos(pi*coeff), sin(pi*coeff) 

572 if not isinstance(cosine, cos) and not isinstance (sine, sin): 

573 return cosine + I*sine 

574 

575 def _eval_rewrite_as_Pow(self, arg, **kwargs): 

576 if arg.is_Mul: 

577 logs = [a for a in arg.args if isinstance(a, log) and len(a.args) == 1] 

578 if logs: 

579 return Pow(logs[0].args[0], arg.coeff(logs[0])) 

580 

581 

582def match_real_imag(expr): 

583 r""" 

584 Try to match expr with $a + Ib$ for real $a$ and $b$. 

585 

586 ``match_real_imag`` returns a tuple containing the real and imaginary 

587 parts of expr or ``(None, None)`` if direct matching is not possible. Contrary 

588 to :func:`~.re()`, :func:`~.im()``, and ``as_real_imag()``, this helper will not force things 

589 by returning expressions themselves containing ``re()`` or ``im()`` and it 

590 does not expand its argument either. 

591 

592 """ 

593 r_, i_ = expr.as_independent(I, as_Add=True) 

594 if i_ == 0 and r_.is_real: 

595 return (r_, i_) 

596 i_ = i_.as_coefficient(I) 

597 if i_ and i_.is_real and r_.is_real: 

598 return (r_, i_) 

599 else: 

600 return (None, None) # simpler to check for than None 

601 

602 

603class log(Function): 

604 r""" 

605 The natural logarithm function `\ln(x)` or `\log(x)`. 

606 

607 Explanation 

608 =========== 

609 

610 Logarithms are taken with the natural base, `e`. To get 

611 a logarithm of a different base ``b``, use ``log(x, b)``, 

612 which is essentially short-hand for ``log(x)/log(b)``. 

613 

614 ``log`` represents the principal branch of the natural 

615 logarithm. As such it has a branch cut along the negative 

616 real axis and returns values having a complex argument in 

617 `(-\pi, \pi]`. 

618 

619 Examples 

620 ======== 

621 

622 >>> from sympy import log, sqrt, S, I 

623 >>> log(8, 2) 

624 3 

625 >>> log(S(8)/3, 2) 

626 -log(3)/log(2) + 3 

627 >>> log(-1 + I*sqrt(3)) 

628 log(2) + 2*I*pi/3 

629 

630 See Also 

631 ======== 

632 

633 exp 

634 

635 """ 

636 

637 args: tTuple[Expr] 

638 

639 _singularities = (S.Zero, S.ComplexInfinity) 

640 

641 def fdiff(self, argindex=1): 

642 """ 

643 Returns the first derivative of the function. 

644 """ 

645 if argindex == 1: 

646 return 1/self.args[0] 

647 else: 

648 raise ArgumentIndexError(self, argindex) 

649 

650 def inverse(self, argindex=1): 

651 r""" 

652 Returns `e^x`, the inverse function of `\log(x)`. 

653 """ 

654 return exp 

655 

656 @classmethod 

657 def eval(cls, arg, base=None): 

658 from sympy.calculus import AccumBounds 

659 from sympy.sets.setexpr import SetExpr 

660 

661 arg = sympify(arg) 

662 

663 if base is not None: 

664 base = sympify(base) 

665 if base == 1: 

666 if arg == 1: 

667 return S.NaN 

668 else: 

669 return S.ComplexInfinity 

670 try: 

671 # handle extraction of powers of the base now 

672 # or else expand_log in Mul would have to handle this 

673 n = multiplicity(base, arg) 

674 if n: 

675 return n + log(arg / base**n) / log(base) 

676 else: 

677 return log(arg)/log(base) 

678 except ValueError: 

679 pass 

680 if base is not S.Exp1: 

681 return cls(arg)/cls(base) 

682 else: 

683 return cls(arg) 

684 

685 if arg.is_Number: 

686 if arg.is_zero: 

687 return S.ComplexInfinity 

688 elif arg is S.One: 

689 return S.Zero 

690 elif arg is S.Infinity: 

691 return S.Infinity 

692 elif arg is S.NegativeInfinity: 

693 return S.Infinity 

694 elif arg is S.NaN: 

695 return S.NaN 

696 elif arg.is_Rational and arg.p == 1: 

697 return -cls(arg.q) 

698 

699 if arg.is_Pow and arg.base is S.Exp1 and arg.exp.is_extended_real: 

700 return arg.exp 

701 if isinstance(arg, exp) and arg.exp.is_extended_real: 

702 return arg.exp 

703 elif isinstance(arg, exp) and arg.exp.is_number: 

704 r_, i_ = match_real_imag(arg.exp) 

705 if i_ and i_.is_comparable: 

706 i_ %= 2*pi 

707 if i_ > pi: 

708 i_ -= 2*pi 

709 return r_ + expand_mul(i_ * I, deep=False) 

710 elif isinstance(arg, exp_polar): 

711 return unpolarify(arg.exp) 

712 elif isinstance(arg, AccumBounds): 

713 if arg.min.is_positive: 

714 return AccumBounds(log(arg.min), log(arg.max)) 

715 elif arg.min.is_zero: 

716 return AccumBounds(S.NegativeInfinity, log(arg.max)) 

717 else: 

718 return S.NaN 

719 elif isinstance(arg, SetExpr): 

720 return arg._eval_func(cls) 

721 

722 if arg.is_number: 

723 if arg.is_negative: 

724 return pi * I + cls(-arg) 

725 elif arg is S.ComplexInfinity: 

726 return S.ComplexInfinity 

727 elif arg is S.Exp1: 

728 return S.One 

729 

730 if arg.is_zero: 

731 return S.ComplexInfinity 

732 

733 # don't autoexpand Pow or Mul (see the issue 3351): 

734 if not arg.is_Add: 

735 coeff = arg.as_coefficient(I) 

736 

737 if coeff is not None: 

738 if coeff is S.Infinity: 

739 return S.Infinity 

740 elif coeff is S.NegativeInfinity: 

741 return S.Infinity 

742 elif coeff.is_Rational: 

743 if coeff.is_nonnegative: 

744 return pi * I * S.Half + cls(coeff) 

745 else: 

746 return -pi * I * S.Half + cls(-coeff) 

747 

748 if arg.is_number and arg.is_algebraic: 

749 # Match arg = coeff*(r_ + i_*I) with coeff>0, r_ and i_ real. 

750 coeff, arg_ = arg.as_independent(I, as_Add=False) 

751 if coeff.is_negative: 

752 coeff *= -1 

753 arg_ *= -1 

754 arg_ = expand_mul(arg_, deep=False) 

755 r_, i_ = arg_.as_independent(I, as_Add=True) 

756 i_ = i_.as_coefficient(I) 

757 if coeff.is_real and i_ and i_.is_real and r_.is_real: 

758 if r_.is_zero: 

759 if i_.is_positive: 

760 return pi * I * S.Half + cls(coeff * i_) 

761 elif i_.is_negative: 

762 return -pi * I * S.Half + cls(coeff * -i_) 

763 else: 

764 from sympy.simplify import ratsimp 

765 # Check for arguments involving rational multiples of pi 

766 t = (i_/r_).cancel() 

767 t1 = (-t).cancel() 

768 atan_table = _log_atan_table() 

769 if t in atan_table: 

770 modulus = ratsimp(coeff * Abs(arg_)) 

771 if r_.is_positive: 

772 return cls(modulus) + I * atan_table[t] 

773 else: 

774 return cls(modulus) + I * (atan_table[t] - pi) 

775 elif t1 in atan_table: 

776 modulus = ratsimp(coeff * Abs(arg_)) 

777 if r_.is_positive: 

778 return cls(modulus) + I * (-atan_table[t1]) 

779 else: 

780 return cls(modulus) + I * (pi - atan_table[t1]) 

781 

782 def as_base_exp(self): 

783 """ 

784 Returns this function in the form (base, exponent). 

785 """ 

786 return self, S.One 

787 

788 @staticmethod 

789 @cacheit 

790 def taylor_term(n, x, *previous_terms): # of log(1+x) 

791 r""" 

792 Returns the next term in the Taylor series expansion of `\log(1+x)`. 

793 """ 

794 from sympy.simplify.powsimp import powsimp 

795 if n < 0: 

796 return S.Zero 

797 x = sympify(x) 

798 if n == 0: 

799 return x 

800 if previous_terms: 

801 p = previous_terms[-1] 

802 if p is not None: 

803 return powsimp((-n) * p * x / (n + 1), deep=True, combine='exp') 

804 return (1 - 2*(n % 2)) * x**(n + 1)/(n + 1) 

805 

806 def _eval_expand_log(self, deep=True, **hints): 

807 from sympy.concrete import Sum, Product 

808 force = hints.get('force', False) 

809 factor = hints.get('factor', False) 

810 if (len(self.args) == 2): 

811 return expand_log(self.func(*self.args), deep=deep, force=force) 

812 arg = self.args[0] 

813 if arg.is_Integer: 

814 # remove perfect powers 

815 p = perfect_power(arg) 

816 logarg = None 

817 coeff = 1 

818 if p is not False: 

819 arg, coeff = p 

820 logarg = self.func(arg) 

821 # expand as product of its prime factors if factor=True 

822 if factor: 

823 p = factorint(arg) 

824 if arg not in p.keys(): 

825 logarg = sum(n*log(val) for val, n in p.items()) 

826 if logarg is not None: 

827 return coeff*logarg 

828 elif arg.is_Rational: 

829 return log(arg.p) - log(arg.q) 

830 elif arg.is_Mul: 

831 expr = [] 

832 nonpos = [] 

833 for x in arg.args: 

834 if force or x.is_positive or x.is_polar: 

835 a = self.func(x) 

836 if isinstance(a, log): 

837 expr.append(self.func(x)._eval_expand_log(**hints)) 

838 else: 

839 expr.append(a) 

840 elif x.is_negative: 

841 a = self.func(-x) 

842 expr.append(a) 

843 nonpos.append(S.NegativeOne) 

844 else: 

845 nonpos.append(x) 

846 return Add(*expr) + log(Mul(*nonpos)) 

847 elif arg.is_Pow or isinstance(arg, exp): 

848 if force or (arg.exp.is_extended_real and (arg.base.is_positive or ((arg.exp+1) 

849 .is_positive and (arg.exp-1).is_nonpositive))) or arg.base.is_polar: 

850 b = arg.base 

851 e = arg.exp 

852 a = self.func(b) 

853 if isinstance(a, log): 

854 return unpolarify(e) * a._eval_expand_log(**hints) 

855 else: 

856 return unpolarify(e) * a 

857 elif isinstance(arg, Product): 

858 if force or arg.function.is_positive: 

859 return Sum(log(arg.function), *arg.limits) 

860 

861 return self.func(arg) 

862 

863 def _eval_simplify(self, **kwargs): 

864 from sympy.simplify.simplify import expand_log, simplify, inversecombine 

865 if len(self.args) == 2: # it's unevaluated 

866 return simplify(self.func(*self.args), **kwargs) 

867 

868 expr = self.func(simplify(self.args[0], **kwargs)) 

869 if kwargs['inverse']: 

870 expr = inversecombine(expr) 

871 expr = expand_log(expr, deep=True) 

872 return min([expr, self], key=kwargs['measure']) 

873 

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

875 """ 

876 Returns this function as a complex coordinate. 

877 

878 Examples 

879 ======== 

880 

881 >>> from sympy import I, log 

882 >>> from sympy.abc import x 

883 >>> log(x).as_real_imag() 

884 (log(Abs(x)), arg(x)) 

885 >>> log(I).as_real_imag() 

886 (0, pi/2) 

887 >>> log(1 + I).as_real_imag() 

888 (log(sqrt(2)), pi/4) 

889 >>> log(I*x).as_real_imag() 

890 (log(Abs(x)), arg(I*x)) 

891 

892 """ 

893 sarg = self.args[0] 

894 if deep: 

895 sarg = self.args[0].expand(deep, **hints) 

896 sarg_abs = Abs(sarg) 

897 if sarg_abs == sarg: 

898 return self, S.Zero 

899 sarg_arg = arg(sarg) 

900 if hints.get('log', False): # Expand the log 

901 hints['complex'] = False 

902 return (log(sarg_abs).expand(deep, **hints), sarg_arg) 

903 else: 

904 return log(sarg_abs), sarg_arg 

905 

906 def _eval_is_rational(self): 

907 s = self.func(*self.args) 

908 if s.func == self.func: 

909 if (self.args[0] - 1).is_zero: 

910 return True 

911 if s.args[0].is_rational and fuzzy_not((self.args[0] - 1).is_zero): 

912 return False 

913 else: 

914 return s.is_rational 

915 

916 def _eval_is_algebraic(self): 

917 s = self.func(*self.args) 

918 if s.func == self.func: 

919 if (self.args[0] - 1).is_zero: 

920 return True 

921 elif fuzzy_not((self.args[0] - 1).is_zero): 

922 if self.args[0].is_algebraic: 

923 return False 

924 else: 

925 return s.is_algebraic 

926 

927 def _eval_is_extended_real(self): 

928 return self.args[0].is_extended_positive 

929 

930 def _eval_is_complex(self): 

931 z = self.args[0] 

932 return fuzzy_and([z.is_complex, fuzzy_not(z.is_zero)]) 

933 

934 def _eval_is_finite(self): 

935 arg = self.args[0] 

936 if arg.is_zero: 

937 return False 

938 return arg.is_finite 

939 

940 def _eval_is_extended_positive(self): 

941 return (self.args[0] - 1).is_extended_positive 

942 

943 def _eval_is_zero(self): 

944 return (self.args[0] - 1).is_zero 

945 

946 def _eval_is_extended_nonnegative(self): 

947 return (self.args[0] - 1).is_extended_nonnegative 

948 

949 def _eval_nseries(self, x, n, logx, cdir=0): 

950 # NOTE Please see the comment at the beginning of this file, labelled 

951 # IMPORTANT. 

952 from sympy.series.order import Order 

953 from sympy.simplify.simplify import logcombine 

954 from sympy.core.symbol import Dummy 

955 

956 if self.args[0] == x: 

957 return log(x) if logx is None else logx 

958 arg = self.args[0] 

959 t = Dummy('t', positive=True) 

960 if cdir == 0: 

961 cdir = 1 

962 z = arg.subs(x, cdir*t) 

963 

964 k, l = Wild("k"), Wild("l") 

965 r = z.match(k*t**l) 

966 if r is not None: 

967 k, l = r[k], r[l] 

968 if l != 0 and not l.has(t) and not k.has(t): 

969 r = l*log(x) if logx is None else l*logx 

970 r += log(k) - l*log(cdir) # XXX true regardless of assumptions? 

971 return r 

972 

973 def coeff_exp(term, x): 

974 coeff, exp = S.One, S.Zero 

975 for factor in Mul.make_args(term): 

976 if factor.has(x): 

977 base, exp = factor.as_base_exp() 

978 if base != x: 

979 try: 

980 return term.leadterm(x) 

981 except ValueError: 

982 return term, S.Zero 

983 else: 

984 coeff *= factor 

985 return coeff, exp 

986 

987 # TODO new and probably slow 

988 try: 

989 a, b = z.leadterm(t, logx=logx, cdir=1) 

990 except (ValueError, NotImplementedError, PoleError): 

991 s = z._eval_nseries(t, n=n, logx=logx, cdir=1) 

992 while s.is_Order: 

993 n += 1 

994 s = z._eval_nseries(t, n=n, logx=logx, cdir=1) 

995 try: 

996 a, b = s.removeO().leadterm(t, cdir=1) 

997 except ValueError: 

998 a, b = s.removeO().as_leading_term(t, cdir=1), S.Zero 

999 

1000 p = (z/(a*t**b) - 1)._eval_nseries(t, n=n, logx=logx, cdir=1) 

1001 if p.has(exp): 

1002 p = logcombine(p) 

1003 if isinstance(p, Order): 

1004 n = p.getn() 

1005 _, d = coeff_exp(p, t) 

1006 logx = log(x) if logx is None else logx 

1007 

1008 if not d.is_positive: 

1009 res = log(a) - b*log(cdir) + b*logx 

1010 _res = res 

1011 logflags = {"deep": True, "log": True, "mul": False, "power_exp": False, 

1012 "power_base": False, "multinomial": False, "basic": False, "force": True, 

1013 "factor": False} 

1014 expr = self.expand(**logflags) 

1015 if (not a.could_extract_minus_sign() and 

1016 logx.could_extract_minus_sign()): 

1017 _res = _res.subs(-logx, -log(x)).expand(**logflags) 

1018 else: 

1019 _res = _res.subs(logx, log(x)).expand(**logflags) 

1020 if _res == expr: 

1021 return res 

1022 return res + Order(x**n, x) 

1023 

1024 def mul(d1, d2): 

1025 res = {} 

1026 for e1, e2 in product(d1, d2): 

1027 ex = e1 + e2 

1028 if ex < n: 

1029 res[ex] = res.get(ex, S.Zero) + d1[e1]*d2[e2] 

1030 return res 

1031 

1032 pterms = {} 

1033 

1034 for term in Add.make_args(p.removeO()): 

1035 co1, e1 = coeff_exp(term, t) 

1036 pterms[e1] = pterms.get(e1, S.Zero) + co1 

1037 

1038 k = S.One 

1039 terms = {} 

1040 pk = pterms 

1041 

1042 while k*d < n: 

1043 coeff = -S.NegativeOne**k/k 

1044 for ex in pk: 

1045 _ = terms.get(ex, S.Zero) + coeff*pk[ex] 

1046 terms[ex] = _.nsimplify() 

1047 pk = mul(pk, pterms) 

1048 k += S.One 

1049 

1050 res = log(a) - b*log(cdir) + b*logx 

1051 for ex in terms: 

1052 res += terms[ex]*t**(ex) 

1053 

1054 if a.is_negative and im(z) != 0: 

1055 from sympy.functions.special.delta_functions import Heaviside 

1056 for i, term in enumerate(z.lseries(t)): 

1057 if not term.is_real or i == 5: 

1058 break 

1059 if i < 5: 

1060 coeff, _ = term.as_coeff_exponent(t) 

1061 res += -2*I*pi*Heaviside(-im(coeff), 0) 

1062 

1063 res = res.subs(t, x/cdir) 

1064 return res + Order(x**n, x) 

1065 

1066 def _eval_as_leading_term(self, x, logx=None, cdir=0): 

1067 # NOTE 

1068 # Refer https://github.com/sympy/sympy/pull/23592 for more information 

1069 # on each of the following steps involved in this method. 

1070 arg0 = self.args[0].together() 

1071 

1072 # STEP 1 

1073 t = Dummy('t', positive=True) 

1074 if cdir == 0: 

1075 cdir = 1 

1076 z = arg0.subs(x, cdir*t) 

1077 

1078 # STEP 2 

1079 try: 

1080 c, e = z.leadterm(t, logx=logx, cdir=1) 

1081 except ValueError: 

1082 arg = arg0.as_leading_term(x, logx=logx, cdir=cdir) 

1083 return log(arg) 

1084 if c.has(t): 

1085 c = c.subs(t, x/cdir) 

1086 if e != 0: 

1087 raise PoleError("Cannot expand %s around 0" % (self)) 

1088 return log(c) 

1089 

1090 # STEP 3 

1091 if c == S.One and e == S.Zero: 

1092 return (arg0 - S.One).as_leading_term(x, logx=logx) 

1093 

1094 # STEP 4 

1095 res = log(c) - e*log(cdir) 

1096 logx = log(x) if logx is None else logx 

1097 res += e*logx 

1098 

1099 # STEP 5 

1100 if c.is_negative and im(z) != 0: 

1101 from sympy.functions.special.delta_functions import Heaviside 

1102 for i, term in enumerate(z.lseries(t)): 

1103 if not term.is_real or i == 5: 

1104 break 

1105 if i < 5: 

1106 coeff, _ = term.as_coeff_exponent(t) 

1107 res += -2*I*pi*Heaviside(-im(coeff), 0) 

1108 return res 

1109 

1110 

1111class LambertW(Function): 

1112 r""" 

1113 The Lambert W function $W(z)$ is defined as the inverse 

1114 function of $w \exp(w)$ [1]_. 

1115 

1116 Explanation 

1117 =========== 

1118 

1119 In other words, the value of $W(z)$ is such that $z = W(z) \exp(W(z))$ 

1120 for any complex number $z$. The Lambert W function is a multivalued 

1121 function with infinitely many branches $W_k(z)$, indexed by 

1122 $k \in \mathbb{Z}$. Each branch gives a different solution $w$ 

1123 of the equation $z = w \exp(w)$. 

1124 

1125 The Lambert W function has two partially real branches: the 

1126 principal branch ($k = 0$) is real for real $z > -1/e$, and the 

1127 $k = -1$ branch is real for $-1/e < z < 0$. All branches except 

1128 $k = 0$ have a logarithmic singularity at $z = 0$. 

1129 

1130 Examples 

1131 ======== 

1132 

1133 >>> from sympy import LambertW 

1134 >>> LambertW(1.2) 

1135 0.635564016364870 

1136 >>> LambertW(1.2, -1).n() 

1137 -1.34747534407696 - 4.41624341514535*I 

1138 >>> LambertW(-1).is_real 

1139 False 

1140 

1141 References 

1142 ========== 

1143 

1144 .. [1] https://en.wikipedia.org/wiki/Lambert_W_function 

1145 """ 

1146 _singularities = (-Pow(S.Exp1, -1, evaluate=False), S.ComplexInfinity) 

1147 

1148 @classmethod 

1149 def eval(cls, x, k=None): 

1150 if k == S.Zero: 

1151 return cls(x) 

1152 elif k is None: 

1153 k = S.Zero 

1154 

1155 if k.is_zero: 

1156 if x.is_zero: 

1157 return S.Zero 

1158 if x is S.Exp1: 

1159 return S.One 

1160 if x == -1/S.Exp1: 

1161 return S.NegativeOne 

1162 if x == -log(2)/2: 

1163 return -log(2) 

1164 if x == 2*log(2): 

1165 return log(2) 

1166 if x == -pi/2: 

1167 return I*pi/2 

1168 if x == exp(1 + S.Exp1): 

1169 return S.Exp1 

1170 if x is S.Infinity: 

1171 return S.Infinity 

1172 if x.is_zero: 

1173 return S.Zero 

1174 

1175 if fuzzy_not(k.is_zero): 

1176 if x.is_zero: 

1177 return S.NegativeInfinity 

1178 if k is S.NegativeOne: 

1179 if x == -pi/2: 

1180 return -I*pi/2 

1181 elif x == -1/S.Exp1: 

1182 return S.NegativeOne 

1183 elif x == -2*exp(-2): 

1184 return -Integer(2) 

1185 

1186 def fdiff(self, argindex=1): 

1187 """ 

1188 Return the first derivative of this function. 

1189 """ 

1190 x = self.args[0] 

1191 

1192 if len(self.args) == 1: 

1193 if argindex == 1: 

1194 return LambertW(x)/(x*(1 + LambertW(x))) 

1195 else: 

1196 k = self.args[1] 

1197 if argindex == 1: 

1198 return LambertW(x, k)/(x*(1 + LambertW(x, k))) 

1199 

1200 raise ArgumentIndexError(self, argindex) 

1201 

1202 def _eval_is_extended_real(self): 

1203 x = self.args[0] 

1204 if len(self.args) == 1: 

1205 k = S.Zero 

1206 else: 

1207 k = self.args[1] 

1208 if k.is_zero: 

1209 if (x + 1/S.Exp1).is_positive: 

1210 return True 

1211 elif (x + 1/S.Exp1).is_nonpositive: 

1212 return False 

1213 elif (k + 1).is_zero: 

1214 if x.is_negative and (x + 1/S.Exp1).is_positive: 

1215 return True 

1216 elif x.is_nonpositive or (x + 1/S.Exp1).is_nonnegative: 

1217 return False 

1218 elif fuzzy_not(k.is_zero) and fuzzy_not((k + 1).is_zero): 

1219 if x.is_extended_real: 

1220 return False 

1221 

1222 def _eval_is_finite(self): 

1223 return self.args[0].is_finite 

1224 

1225 def _eval_is_algebraic(self): 

1226 s = self.func(*self.args) 

1227 if s.func == self.func: 

1228 if fuzzy_not(self.args[0].is_zero) and self.args[0].is_algebraic: 

1229 return False 

1230 else: 

1231 return s.is_algebraic 

1232 

1233 def _eval_as_leading_term(self, x, logx=None, cdir=0): 

1234 if len(self.args) == 1: 

1235 arg = self.args[0] 

1236 arg0 = arg.subs(x, 0).cancel() 

1237 if not arg0.is_zero: 

1238 return self.func(arg0) 

1239 return arg.as_leading_term(x) 

1240 

1241 def _eval_nseries(self, x, n, logx, cdir=0): 

1242 if len(self.args) == 1: 

1243 from sympy.functions.elementary.integers import ceiling 

1244 from sympy.series.order import Order 

1245 arg = self.args[0].nseries(x, n=n, logx=logx) 

1246 lt = arg.as_leading_term(x, logx=logx) 

1247 lte = 1 

1248 if lt.is_Pow: 

1249 lte = lt.exp 

1250 if ceiling(n/lte) >= 1: 

1251 s = Add(*[(-S.One)**(k - 1)*Integer(k)**(k - 2)/ 

1252 factorial(k - 1)*arg**k for k in range(1, ceiling(n/lte))]) 

1253 s = expand_multinomial(s) 

1254 else: 

1255 s = S.Zero 

1256 

1257 return s + Order(x**n, x) 

1258 return super()._eval_nseries(x, n, logx) 

1259 

1260 def _eval_is_zero(self): 

1261 x = self.args[0] 

1262 if len(self.args) == 1: 

1263 return x.is_zero 

1264 else: 

1265 return fuzzy_and([x.is_zero, self.args[1].is_zero]) 

1266 

1267 

1268@cacheit 

1269def _log_atan_table(): 

1270 return { 

1271 # first quadrant only 

1272 sqrt(3): pi / 3, 

1273 1: pi / 4, 

1274 sqrt(5 - 2 * sqrt(5)): pi / 5, 

1275 sqrt(2) * sqrt(5 - sqrt(5)) / (1 + sqrt(5)): pi / 5, 

1276 sqrt(5 + 2 * sqrt(5)): pi * Rational(2, 5), 

1277 sqrt(2) * sqrt(sqrt(5) + 5) / (-1 + sqrt(5)): pi * Rational(2, 5), 

1278 sqrt(3) / 3: pi / 6, 

1279 sqrt(2) - 1: pi / 8, 

1280 sqrt(2 - sqrt(2)) / sqrt(sqrt(2) + 2): pi / 8, 

1281 sqrt(2) + 1: pi * Rational(3, 8), 

1282 sqrt(sqrt(2) + 2) / sqrt(2 - sqrt(2)): pi * Rational(3, 8), 

1283 sqrt(1 - 2 * sqrt(5) / 5): pi / 10, 

1284 (-sqrt(2) + sqrt(10)) / (2 * sqrt(sqrt(5) + 5)): pi / 10, 

1285 sqrt(1 + 2 * sqrt(5) / 5): pi * Rational(3, 10), 

1286 (sqrt(2) + sqrt(10)) / (2 * sqrt(5 - sqrt(5))): pi * Rational(3, 10), 

1287 2 - sqrt(3): pi / 12, 

1288 (-1 + sqrt(3)) / (1 + sqrt(3)): pi / 12, 

1289 2 + sqrt(3): pi * Rational(5, 12), 

1290 (1 + sqrt(3)) / (-1 + sqrt(3)): pi * Rational(5, 12) 

1291 }