Coverage for /usr/lib/python3/dist-packages/sympy/functions/special/error

1""" This module contains various functions that are special cases

2 of incomplete gamma functions. It should probably be renamed. """

4from sympy.core import EulerGamma # Must be imported from core, not core.numbers

5from sympy.core.add import Add

6from sympy.core.cache import cacheit

7from sympy.core.function import Function, ArgumentIndexError, expand_mul

8from sympy.core.numbers import I, pi, Rational

9from sympy.core.relational import is_eq

10from sympy.core.power import Pow

11from sympy.core.singleton import S

12from sympy.core.symbol import Symbol

13from sympy.core.sympify import sympify

14from sympy.functions.combinatorial.factorials import factorial, factorial2, RisingFactorial

15from sympy.functions.elementary.complexes import polar_lift, re, unpolarify

16from sympy.functions.elementary.integers import ceiling, floor

17from sympy.functions.elementary.miscellaneous import sqrt, root

18from sympy.functions.elementary.exponential import exp, log, exp_polar

19from sympy.functions.elementary.hyperbolic import cosh, sinh

20from sympy.functions.elementary.trigonometric import cos, sin, sinc

21from sympy.functions.special.hyper import hyper, meijerg

23# TODO series expansions

24# TODO see the "Note:" in Ei

26# Helper function

27def real_to_real_as_real_imag(self, deep=True, **hints):

28 if self.args[0].is_extended_real:

29 if deep:

30 hints['complex'] = False

31 return (self.expand(deep, **hints), S.Zero)

32 else:

33 return (self, S.Zero)

34 if deep:

35 x, y = self.args[0].expand(deep, **hints).as_real_imag()

36 else:

37 x, y = self.args[0].as_real_imag()

38 re = (self.func(x + I*y) + self.func(x - I*y))/2

39 im = (self.func(x + I*y) - self.func(x - I*y))/(2*I)

40 return (re, im)

43###############################################################################

44################################ ERROR FUNCTION ###############################

45###############################################################################

48class erf(Function):

49 r"""

50 The Gauss error function.

52 Explanation

53 ===========

55 This function is defined as:

57 .. math ::

58 \mathrm{erf}(x) = \frac{2}{\sqrt{\pi}} \int_0^x e^{-t^2} \mathrm{d}t.

60 Examples

61 ========

63 >>> from sympy import I, oo, erf

64 >>> from sympy.abc import z

66 Several special values are known:

68 >>> erf(0)

69 0

70 >>> erf(oo)

71 1

72 >>> erf(-oo)

73 -1

74 >>> erf(I*oo)

75 oo*I

76 >>> erf(-I*oo)

77 -oo*I

79 In general one can pull out factors of -1 and $I$ from the argument:

81 >>> erf(-z)

82 -erf(z)

84 The error function obeys the mirror symmetry:

86 >>> from sympy import conjugate

87 >>> conjugate(erf(z))

88 erf(conjugate(z))

90 Differentiation with respect to $z$ is supported:

92 >>> from sympy import diff

93 >>> diff(erf(z), z)

94 2*exp(-z**2)/sqrt(pi)

96 We can numerically evaluate the error function to arbitrary precision

97 on the whole complex plane:

99 >>> erf(4).evalf(30)

100 0.999999984582742099719981147840

101

102 >>> erf(-4*I).evalf(30)

103 -1296959.73071763923152794095062*I

104

105 See Also

106 ========

107

108 erfc: Complementary error function.

109 erfi: Imaginary error function.

110 erf2: Two-argument error function.

111 erfinv: Inverse error function.

112 erfcinv: Inverse Complementary error function.

113 erf2inv: Inverse two-argument error function.

114

115 References

116 ==========

117

118 .. [1] https://en.wikipedia.org/wiki/Error_function

119 .. [2] https://dlmf.nist.gov/7

120 .. [3] https://mathworld.wolfram.com/Erf.html

121 .. [4] https://functions.wolfram.com/GammaBetaErf/Erf

122

123 """

124

125 unbranched = True

126

127 def fdiff(self, argindex=1):

128 if argindex == 1:

129 return 2*exp(-self.args[0]**2)/sqrt(pi)

130 else:

131 raise ArgumentIndexError(self, argindex)

132

133

134 def inverse(self, argindex=1):

135 """

136 Returns the inverse of this function.

137

138 """

139 return erfinv

140

141 @classmethod

142 def eval(cls, arg):

143 if arg.is_Number:

144 if arg is S.NaN:

145 return S.NaN

146 elif arg is S.Infinity:

147 return S.One

148 elif arg is S.NegativeInfinity:

149 return S.NegativeOne

150 elif arg.is_zero:

151 return S.Zero

152

153 if isinstance(arg, erfinv):

154 return arg.args[0]

155

156 if isinstance(arg, erfcinv):

157 return S.One - arg.args[0]

158

159 if arg.is_zero:

160 return S.Zero

161

162 # Only happens with unevaluated erf2inv

163 if isinstance(arg, erf2inv) and arg.args[0].is_zero:

164 return arg.args[1]

165

166 # Try to pull out factors of I

167 t = arg.extract_multiplicatively(I)

168 if t in (S.Infinity, S.NegativeInfinity):

169 return arg

170

171 # Try to pull out factors of -1

172 if arg.could_extract_minus_sign():

173 return -cls(-arg)

174

175 @staticmethod

176 @cacheit

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

178 if n < 0 or n % 2 == 0:

179 return S.Zero

180 else:

181 x = sympify(x)

182 k = floor((n - 1)/S(2))

183 if len(previous_terms) > 2:

184 return -previous_terms[-2] * x**2 * (n - 2)/(n*k)

185 else:

186 return 2*S.NegativeOne**k * x**n/(n*factorial(k)*sqrt(pi))

187

188 def _eval_conjugate(self):

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

190

191 def _eval_is_real(self):

192 return self.args[0].is_extended_real

193

194 def _eval_is_finite(self):

195 if self.args[0].is_finite:

196 return True

197 else:

198 return self.args[0].is_extended_real

199

200 def _eval_is_zero(self):

201 return self.args[0].is_zero

202

203 def _eval_rewrite_as_uppergamma(self, z, **kwargs):

204 from sympy.functions.special.gamma_functions import uppergamma

205 return sqrt(z**2)/z*(S.One - uppergamma(S.Half, z**2)/sqrt(pi))

206

207 def _eval_rewrite_as_fresnels(self, z, **kwargs):

208 arg = (S.One - I)*z/sqrt(pi)

209 return (S.One + I)*(fresnelc(arg) - I*fresnels(arg))

210

211 def _eval_rewrite_as_fresnelc(self, z, **kwargs):

212 arg = (S.One - I)*z/sqrt(pi)

213 return (S.One + I)*(fresnelc(arg) - I*fresnels(arg))

214

215 def _eval_rewrite_as_meijerg(self, z, **kwargs):

216 return z/sqrt(pi)*meijerg([S.Half], [], [0], [Rational(-1, 2)], z**2)

217

218 def _eval_rewrite_as_hyper(self, z, **kwargs):

219 return 2*z/sqrt(pi)*hyper([S.Half], [3*S.Half], -z**2)

220

221 def _eval_rewrite_as_expint(self, z, **kwargs):

222 return sqrt(z**2)/z - z*expint(S.Half, z**2)/sqrt(pi)

223

224 def _eval_rewrite_as_tractable(self, z, limitvar=None, **kwargs):

225 from sympy.series.limits import limit

226 if limitvar:

227 lim = limit(z, limitvar, S.Infinity)

228 if lim is S.NegativeInfinity:

229 return S.NegativeOne + _erfs(-z)*exp(-z**2)

230 return S.One - _erfs(z)*exp(-z**2)

231

232 def _eval_rewrite_as_erfc(self, z, **kwargs):

233 return S.One - erfc(z)

234

235 def _eval_rewrite_as_erfi(self, z, **kwargs):

236 return -I*erfi(I*z)

237

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

239 arg = self.args[0].as_leading_term(x, logx=logx, cdir=cdir)

240 arg0 = arg.subs(x, 0)

241

242 if arg0 is S.ComplexInfinity:

243 arg0 = arg.limit(x, 0, dir='-' if cdir == -1 else '+')

244 if x in arg.free_symbols and arg0.is_zero:

245 return 2*arg/sqrt(pi)

246 else:

247 return self.func(arg0)

248

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

250 from sympy.series.order import Order

251 point = args0[0]

252

253 if point in [S.Infinity, S.NegativeInfinity]:

254 z = self.args[0]

255

256 try:

257 _, ex = z.leadterm(x)

258 except (ValueError, NotImplementedError):

259 return self

260

261 ex = -ex # as x->1/x for aseries

262 if ex.is_positive:

263 newn = ceiling(n/ex)

264 s = [S.NegativeOne**k * factorial2(2*k - 1) / (z**(2*k + 1) * 2**k)

265 for k in range(newn)] + [Order(1/z**newn, x)]

266 return S.One - (exp(-z**2)/sqrt(pi)) * Add(*s)

267

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

269

270 as_real_imag = real_to_real_as_real_imag

271

272

273class erfc(Function):

274 r"""

275 Complementary Error Function.

276

277 Explanation

278 ===========

279

280 The function is defined as:

281

282 .. math ::

283 \mathrm{erfc}(x) = \frac{2}{\sqrt{\pi}} \int_x^\infty e^{-t^2} \mathrm{d}t

284

285 Examples

286 ========

287

288 >>> from sympy import I, oo, erfc

289 >>> from sympy.abc import z

290

291 Several special values are known:

292

293 >>> erfc(0)

294 1

295 >>> erfc(oo)

296 0

297 >>> erfc(-oo)

298 2

299 >>> erfc(I*oo)

300 -oo*I

301 >>> erfc(-I*oo)

302 oo*I

303

304 The error function obeys the mirror symmetry:

305

306 >>> from sympy import conjugate

307 >>> conjugate(erfc(z))

308 erfc(conjugate(z))

309

310 Differentiation with respect to $z$ is supported:

311

312 >>> from sympy import diff

313 >>> diff(erfc(z), z)

314 -2*exp(-z**2)/sqrt(pi)

315

316 It also follows

317

318 >>> erfc(-z)

319 2 - erfc(z)

320

321 We can numerically evaluate the complementary error function to arbitrary

322 precision on the whole complex plane:

323

324 >>> erfc(4).evalf(30)

325 0.0000000154172579002800188521596734869

326

327 >>> erfc(4*I).evalf(30)

328 1.0 - 1296959.73071763923152794095062*I

329

330 See Also

331 ========

332

333 erf: Gaussian error function.

334 erfi: Imaginary error function.

335 erf2: Two-argument error function.

336 erfinv: Inverse error function.

337 erfcinv: Inverse Complementary error function.

338 erf2inv: Inverse two-argument error function.

339

340 References

341 ==========

342

343 .. [1] https://en.wikipedia.org/wiki/Error_function

344 .. [2] https://dlmf.nist.gov/7

345 .. [3] https://mathworld.wolfram.com/Erfc.html

346 .. [4] https://functions.wolfram.com/GammaBetaErf/Erfc

347

348 """

349

350 unbranched = True

351

352 def fdiff(self, argindex=1):

353 if argindex == 1:

354 return -2*exp(-self.args[0]**2)/sqrt(pi)

355 else:

356 raise ArgumentIndexError(self, argindex)

357

358 def inverse(self, argindex=1):

359 """

360 Returns the inverse of this function.

361

362 """

363 return erfcinv

364

365 @classmethod

366 def eval(cls, arg):

367 if arg.is_Number:

368 if arg is S.NaN:

369 return S.NaN

370 elif arg is S.Infinity:

371 return S.Zero

372 elif arg.is_zero:

373 return S.One

374

375 if isinstance(arg, erfinv):

376 return S.One - arg.args[0]

377

378 if isinstance(arg, erfcinv):

379 return arg.args[0]

380

381 if arg.is_zero:

382 return S.One

383

384 # Try to pull out factors of I

385 t = arg.extract_multiplicatively(I)

386 if t in (S.Infinity, S.NegativeInfinity):

387 return -arg

388

389 # Try to pull out factors of -1

390 if arg.could_extract_minus_sign():

391 return 2 - cls(-arg)

392

393 @staticmethod

394 @cacheit

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

396 if n == 0:

397 return S.One

398 elif n < 0 or n % 2 == 0:

399 return S.Zero

400 else:

401 x = sympify(x)

402 k = floor((n - 1)/S(2))

403 if len(previous_terms) > 2:

404 return -previous_terms[-2] * x**2 * (n - 2)/(n*k)

405 else:

406 return -2*S.NegativeOne**k * x**n/(n*factorial(k)*sqrt(pi))

407

408 def _eval_conjugate(self):

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

410

411 def _eval_is_real(self):

412 return self.args[0].is_extended_real

413

414 def _eval_rewrite_as_tractable(self, z, limitvar=None, **kwargs):

415 return self.rewrite(erf).rewrite("tractable", deep=True, limitvar=limitvar)

416

417 def _eval_rewrite_as_erf(self, z, **kwargs):

418 return S.One - erf(z)

419

420 def _eval_rewrite_as_erfi(self, z, **kwargs):

421 return S.One + I*erfi(I*z)

422

423 def _eval_rewrite_as_fresnels(self, z, **kwargs):

424 arg = (S.One - I)*z/sqrt(pi)

425 return S.One - (S.One + I)*(fresnelc(arg) - I*fresnels(arg))

426

427 def _eval_rewrite_as_fresnelc(self, z, **kwargs):

428 arg = (S.One-I)*z/sqrt(pi)

429 return S.One - (S.One + I)*(fresnelc(arg) - I*fresnels(arg))

430

431 def _eval_rewrite_as_meijerg(self, z, **kwargs):

432 return S.One - z/sqrt(pi)*meijerg([S.Half], [], [0], [Rational(-1, 2)], z**2)

433

434 def _eval_rewrite_as_hyper(self, z, **kwargs):

435 return S.One - 2*z/sqrt(pi)*hyper([S.Half], [3*S.Half], -z**2)

436

437 def _eval_rewrite_as_uppergamma(self, z, **kwargs):

438 from sympy.functions.special.gamma_functions import uppergamma

439 return S.One - sqrt(z**2)/z*(S.One - uppergamma(S.Half, z**2)/sqrt(pi))

440

441 def _eval_rewrite_as_expint(self, z, **kwargs):

442 return S.One - sqrt(z**2)/z + z*expint(S.Half, z**2)/sqrt(pi)

443

444 def _eval_expand_func(self, **hints):

445 return self.rewrite(erf)

446

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

448 arg = self.args[0].as_leading_term(x, logx=logx, cdir=cdir)

449 arg0 = arg.subs(x, 0)

450

451 if arg0 is S.ComplexInfinity:

452 arg0 = arg.limit(x, 0, dir='-' if cdir == -1 else '+')

453 if arg0.is_zero:

454 return S.One

455 else:

456 return self.func(arg0)

457

458 as_real_imag = real_to_real_as_real_imag

459

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

461 return S.One - erf(*self.args)._eval_aseries(n, args0, x, logx)

462

463

464class erfi(Function):

465 r"""

466 Imaginary error function.

467

468 Explanation

469 ===========

470

471 The function erfi is defined as:

472

473 .. math ::

474 \mathrm{erfi}(x) = \frac{2}{\sqrt{\pi}} \int_0^x e^{t^2} \mathrm{d}t

475

476 Examples

477 ========

478

479 >>> from sympy import I, oo, erfi

480 >>> from sympy.abc import z

481

482 Several special values are known:

483

484 >>> erfi(0)

485 0

486 >>> erfi(oo)

487 oo

488 >>> erfi(-oo)

489 -oo

490 >>> erfi(I*oo)

491 I

492 >>> erfi(-I*oo)

493 -I

494

495 In general one can pull out factors of -1 and $I$ from the argument:

496

497 >>> erfi(-z)

498 -erfi(z)

499

500 >>> from sympy import conjugate

501 >>> conjugate(erfi(z))

502 erfi(conjugate(z))

503

504 Differentiation with respect to $z$ is supported:

505

506 >>> from sympy import diff

507 >>> diff(erfi(z), z)

508 2*exp(z**2)/sqrt(pi)

509

510 We can numerically evaluate the imaginary error function to arbitrary

511 precision on the whole complex plane:

512

513 >>> erfi(2).evalf(30)

514 18.5648024145755525987042919132

515

516 >>> erfi(-2*I).evalf(30)

517 -0.995322265018952734162069256367*I

518

519 See Also

520 ========

521

522 erf: Gaussian error function.

523 erfc: Complementary error function.

524 erf2: Two-argument error function.

525 erfinv: Inverse error function.

526 erfcinv: Inverse Complementary error function.

527 erf2inv: Inverse two-argument error function.

528

529 References

530 ==========

531

532 .. [1] https://en.wikipedia.org/wiki/Error_function

533 .. [2] https://mathworld.wolfram.com/Erfi.html

534 .. [3] https://functions.wolfram.com/GammaBetaErf/Erfi

535

536 """

537

538 unbranched = True

539

540 def fdiff(self, argindex=1):

541 if argindex == 1:

542 return 2*exp(self.args[0]**2)/sqrt(pi)

543 else:

544 raise ArgumentIndexError(self, argindex)

545

546 @classmethod

547 def eval(cls, z):

548 if z.is_Number:

549 if z is S.NaN:

550 return S.NaN

551 elif z.is_zero:

552 return S.Zero

553 elif z is S.Infinity:

554 return S.Infinity

555

556 if z.is_zero:

557 return S.Zero

558

559 # Try to pull out factors of -1

560 if z.could_extract_minus_sign():

561 return -cls(-z)

562

563 # Try to pull out factors of I

564 nz = z.extract_multiplicatively(I)

565 if nz is not None:

566 if nz is S.Infinity:

567 return I

568 if isinstance(nz, erfinv):

569 return I*nz.args[0]

570 if isinstance(nz, erfcinv):

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

572 # Only happens with unevaluated erf2inv

573 if isinstance(nz, erf2inv) and nz.args[0].is_zero:

574 return I*nz.args[1]

575

576 @staticmethod

577 @cacheit

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

579 if n < 0 or n % 2 == 0:

580 return S.Zero

581 else:

582 x = sympify(x)

583 k = floor((n - 1)/S(2))

584 if len(previous_terms) > 2:

585 return previous_terms[-2] * x**2 * (n - 2)/(n*k)

586 else:

587 return 2 * x**n/(n*factorial(k)*sqrt(pi))

588

589 def _eval_conjugate(self):

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

591

592 def _eval_is_extended_real(self):

593 return self.args[0].is_extended_real

594

595 def _eval_is_zero(self):

596 return self.args[0].is_zero

597

598 def _eval_rewrite_as_tractable(self, z, limitvar=None, **kwargs):

599 return self.rewrite(erf).rewrite("tractable", deep=True, limitvar=limitvar)

600

601 def _eval_rewrite_as_erf(self, z, **kwargs):

602 return -I*erf(I*z)

603

604 def _eval_rewrite_as_erfc(self, z, **kwargs):

605 return I*erfc(I*z) - I

606

607 def _eval_rewrite_as_fresnels(self, z, **kwargs):

608 arg = (S.One + I)*z/sqrt(pi)

609 return (S.One - I)*(fresnelc(arg) - I*fresnels(arg))

610

611 def _eval_rewrite_as_fresnelc(self, z, **kwargs):

612 arg = (S.One + I)*z/sqrt(pi)

613 return (S.One - I)*(fresnelc(arg) - I*fresnels(arg))

614

615 def _eval_rewrite_as_meijerg(self, z, **kwargs):

616 return z/sqrt(pi)*meijerg([S.Half], [], [0], [Rational(-1, 2)], -z**2)

617

618 def _eval_rewrite_as_hyper(self, z, **kwargs):

619 return 2*z/sqrt(pi)*hyper([S.Half], [3*S.Half], z**2)

620

621 def _eval_rewrite_as_uppergamma(self, z, **kwargs):

622 from sympy.functions.special.gamma_functions import uppergamma

623 return sqrt(-z**2)/z*(uppergamma(S.Half, -z**2)/sqrt(pi) - S.One)

624

625 def _eval_rewrite_as_expint(self, z, **kwargs):

626 return sqrt(-z**2)/z - z*expint(S.Half, -z**2)/sqrt(pi)

627

628 def _eval_expand_func(self, **hints):

629 return self.rewrite(erf)

630

631 as_real_imag = real_to_real_as_real_imag

632

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

634 arg = self.args[0].as_leading_term(x, logx=logx, cdir=cdir)

635 arg0 = arg.subs(x, 0)

636

637 if x in arg.free_symbols and arg0.is_zero:

638 return 2*arg/sqrt(pi)

639 elif arg0.is_finite:

640 return self.func(arg0)

641 return self.func(arg)

642

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

644 from sympy.series.order import Order

645 point = args0[0]

646

647 if point is S.Infinity:

648 z = self.args[0]

649 s = [factorial2(2*k - 1) / (2**k * z**(2*k + 1))

650 for k in range(n)] + [Order(1/z**n, x)]

651 return -I + (exp(z**2)/sqrt(pi)) * Add(*s)

652

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

654

655

656class erf2(Function):

657 r"""

658 Two-argument error function.

659

660 Explanation

661 ===========

662

663 This function is defined as:

664

665 .. math ::

666 \mathrm{erf2}(x, y) = \frac{2}{\sqrt{\pi}} \int_x^y e^{-t^2} \mathrm{d}t

667

668 Examples

669 ========

670

671 >>> from sympy import oo, erf2

672 >>> from sympy.abc import x, y

673

674 Several special values are known:

675

676 >>> erf2(0, 0)

677 0

678 >>> erf2(x, x)

679 0

680 >>> erf2(x, oo)

681 1 - erf(x)

682 >>> erf2(x, -oo)

683 -erf(x) - 1

684 >>> erf2(oo, y)

685 erf(y) - 1

686 >>> erf2(-oo, y)

687 erf(y) + 1

688

689 In general one can pull out factors of -1:

690

691 >>> erf2(-x, -y)

692 -erf2(x, y)

693

694 The error function obeys the mirror symmetry:

695

696 >>> from sympy import conjugate

697 >>> conjugate(erf2(x, y))

698 erf2(conjugate(x), conjugate(y))

699

700 Differentiation with respect to $x$, $y$ is supported:

701

702 >>> from sympy import diff

703 >>> diff(erf2(x, y), x)

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

705 >>> diff(erf2(x, y), y)

706 2*exp(-y**2)/sqrt(pi)

707

708 See Also

709 ========

710

711 erf: Gaussian error function.

712 erfc: Complementary error function.

713 erfi: Imaginary error function.

714 erfinv: Inverse error function.

715 erfcinv: Inverse Complementary error function.

716 erf2inv: Inverse two-argument error function.

717

718 References

719 ==========

720

721 .. [1] https://functions.wolfram.com/GammaBetaErf/Erf2/

722

723 """

724

725

726 def fdiff(self, argindex):

727 x, y = self.args

728 if argindex == 1:

729 return -2*exp(-x**2)/sqrt(pi)

730 elif argindex == 2:

731 return 2*exp(-y**2)/sqrt(pi)

732 else:

733 raise ArgumentIndexError(self, argindex)

734

735 @classmethod

736 def eval(cls, x, y):

737 chk = (S.Infinity, S.NegativeInfinity, S.Zero)

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

739 return S.NaN

740 elif x == y:

741 return S.Zero

742 elif x in chk or y in chk:

743 return erf(y) - erf(x)

744

745 if isinstance(y, erf2inv) and y.args[0] == x:

746 return y.args[1]

747

748 if x.is_zero or y.is_zero or x.is_extended_real and x.is_infinite or \

749 y.is_extended_real and y.is_infinite:

750 return erf(y) - erf(x)

751

752 #Try to pull out -1 factor

753 sign_x = x.could_extract_minus_sign()

754 sign_y = y.could_extract_minus_sign()

755 if (sign_x and sign_y):

756 return -cls(-x, -y)

757 elif (sign_x or sign_y):

758 return erf(y)-erf(x)

759

760 def _eval_conjugate(self):

761 return self.func(self.args[0].conjugate(), self.args[1].conjugate())

762

763 def _eval_is_extended_real(self):

764 return self.args[0].is_extended_real and self.args[1].is_extended_real

765

766 def _eval_rewrite_as_erf(self, x, y, **kwargs):

767 return erf(y) - erf(x)

768

769 def _eval_rewrite_as_erfc(self, x, y, **kwargs):

770 return erfc(x) - erfc(y)

771

772 def _eval_rewrite_as_erfi(self, x, y, **kwargs):

773 return I*(erfi(I*x)-erfi(I*y))

774

775 def _eval_rewrite_as_fresnels(self, x, y, **kwargs):

776 return erf(y).rewrite(fresnels) - erf(x).rewrite(fresnels)

777

778 def _eval_rewrite_as_fresnelc(self, x, y, **kwargs):

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

780

781 def _eval_rewrite_as_meijerg(self, x, y, **kwargs):

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

783

784 def _eval_rewrite_as_hyper(self, x, y, **kwargs):

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

786

787 def _eval_rewrite_as_uppergamma(self, x, y, **kwargs):

788 from sympy.functions.special.gamma_functions import uppergamma

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

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

791

792 def _eval_rewrite_as_expint(self, x, y, **kwargs):

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

794

795 def _eval_expand_func(self, **hints):

796 return self.rewrite(erf)

797

798 def _eval_is_zero(self):

799 return is_eq(*self.args)

800

801class erfinv(Function):

802 r"""

803 Inverse Error Function. The erfinv function is defined as:

804

805 .. math ::

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

807

808 Examples

809 ========

810

811 >>> from sympy import erfinv

812 >>> from sympy.abc import x

813

814 Several special values are known:

815

816 >>> erfinv(0)

817 0

818 >>> erfinv(1)

819 oo

820

821 Differentiation with respect to $x$ is supported:

822

823 >>> from sympy import diff

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

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

826

827 We can numerically evaluate the inverse error function to arbitrary

828 precision on [-1, 1]:

829

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

831 0.179143454621291692285822705344

832

833 See Also

834 ========

835

836 erf: Gaussian error function.

837 erfc: Complementary error function.

838 erfi: Imaginary error function.

839 erf2: Two-argument error function.

840 erfcinv: Inverse Complementary error function.

841 erf2inv: Inverse two-argument error function.

842

843 References

844 ==========

845

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

847 .. [2] https://functions.wolfram.com/GammaBetaErf/InverseErf/

848

849 """

850

851

852 def fdiff(self, argindex =1):

853 if argindex == 1:

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

855 else :

856 raise ArgumentIndexError(self, argindex)

857

858 def inverse(self, argindex=1):

859 """

860 Returns the inverse of this function.

861

862 """

863 return erf

864

865 @classmethod

866 def eval(cls, z):

867 if z is S.NaN:

868 return S.NaN

869 elif z is S.NegativeOne:

870 return S.NegativeInfinity

871 elif z.is_zero:

872 return S.Zero

873 elif z is S.One:

874 return S.Infinity

875

876 if isinstance(z, erf) and z.args[0].is_extended_real:

877 return z.args[0]

878

879 if z.is_zero:

880 return S.Zero

881

882 # Try to pull out factors of -1

883 nz = z.extract_multiplicatively(-1)

884 if nz is not None and (isinstance(nz, erf) and (nz.args[0]).is_extended_real):

885 return -nz.args[0]

886

887 def _eval_rewrite_as_erfcinv(self, z, **kwargs):

888 return erfcinv(1-z)

889

890 def _eval_is_zero(self):

891 return self.args[0].is_zero

892

893

894class erfcinv (Function):

895 r"""

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

897

898 .. math ::

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

900

901 Examples

902 ========

903

904 >>> from sympy import erfcinv

905 >>> from sympy.abc import x

906

907 Several special values are known:

908

909 >>> erfcinv(1)

910 0

911 >>> erfcinv(0)

912 oo

913

914 Differentiation with respect to $x$ is supported:

915

916 >>> from sympy import diff

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

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

919

920 See Also

921 ========

922

923 erf: Gaussian error function.

924 erfc: Complementary error function.

925 erfi: Imaginary error function.

926 erf2: Two-argument error function.

927 erfinv: Inverse error function.

928 erf2inv: Inverse two-argument error function.

929

930 References

931 ==========

932

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

934 .. [2] https://functions.wolfram.com/GammaBetaErf/InverseErfc/

935

936 """

937

938

939 def fdiff(self, argindex =1):

940 if argindex == 1:

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

942 else:

943 raise ArgumentIndexError(self, argindex)

944

945 def inverse(self, argindex=1):

946 """

947 Returns the inverse of this function.

948

949 """

950 return erfc

951

952 @classmethod

953 def eval(cls, z):

954 if z is S.NaN:

955 return S.NaN

956 elif z.is_zero:

957 return S.Infinity

958 elif z is S.One:

959 return S.Zero

960 elif z == 2:

961 return S.NegativeInfinity

962

963 if z.is_zero:

964 return S.Infinity

965

966 def _eval_rewrite_as_erfinv(self, z, **kwargs):

967 return erfinv(1-z)

968

969 def _eval_is_zero(self):

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

971

972 def _eval_is_infinite(self):

973 return self.args[0].is_zero

974

975

976class erf2inv(Function):

977 r"""

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

979

980 .. math ::

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

982

983 Examples

984 ========

985

986 >>> from sympy import erf2inv, oo

987 >>> from sympy.abc import x, y

988

989 Several special values are known:

990

991 >>> erf2inv(0, 0)

992 0

993 >>> erf2inv(1, 0)

994 1

995 >>> erf2inv(0, 1)

996 oo

997 >>> erf2inv(0, y)

998 erfinv(y)

999 >>> erf2inv(oo, y)

1000 erfcinv(-y)

1001

1002 Differentiation with respect to $x$ and $y$ is supported:

1003

1004 >>> from sympy import diff

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

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

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

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

1009

1010 See Also

1011 ========

1012

1013 erf: Gaussian error function.

1014 erfc: Complementary error function.

1015 erfi: Imaginary error function.

1016 erf2: Two-argument error function.

1017 erfinv: Inverse error function.

1018 erfcinv: Inverse complementary error function.

1019

1020 References

1021 ==========

1022

1023 .. [1] https://functions.wolfram.com/GammaBetaErf/InverseErf2/

1024

1025 """

1026

1027

1028 def fdiff(self, argindex):

1029 x, y = self.args

1030 if argindex == 1:

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

1032 elif argindex == 2:

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

1034 else:

1035 raise ArgumentIndexError(self, argindex)

1036

1037 @classmethod

1038 def eval(cls, x, y):

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

1040 return S.NaN

1041 elif x.is_zero and y.is_zero:

1042 return S.Zero

1043 elif x.is_zero and y is S.One:

1044 return S.Infinity

1045 elif x is S.One and y.is_zero:

1046 return S.One

1047 elif x.is_zero:

1048 return erfinv(y)

1049 elif x is S.Infinity:

1050 return erfcinv(-y)

1051 elif y.is_zero:

1052 return x

1053 elif y is S.Infinity:

1054 return erfinv(x)

1055

1056 if x.is_zero:

1057 if y.is_zero:

1058 return S.Zero

1059 else:

1060 return erfinv(y)

1061 if y.is_zero:

1062 return x

1063

1064 def _eval_is_zero(self):

1065 x, y = self.args

1066 if x.is_zero and y.is_zero:

1067 return True

1068

1069###############################################################################

1070#################### EXPONENTIAL INTEGRALS ####################################

1071###############################################################################

1072

1073class Ei(Function):

1074 r"""

1075 The classical exponential integral.

1076

1077 Explanation

1078 ===========

1079

1080 For use in SymPy, this function is defined as

1081

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

1083 + \log(x) + \gamma,

1084

1085 where $\gamma$ is the Euler-Mascheroni constant.

1086

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

1088 Riemann surface of the logarithm. Otherwise this defines an analytic

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

1090

1091 **Background**

1092

1093 The name exponential integral comes from the following statement:

1094

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

1096

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

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

1099

1100 Examples

1101 ========

1102

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

1104 >>> from sympy.abc import x

1105

1106 >>> Ei(-1)

1107 Ei(-1)

1108

1109 This yields a real value:

1110

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

1112 -0.219383934395520

1113

1114 On the other hand the analytic continuation is not real:

1115

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

1117 -0.21938393439552 + 3.14159265358979*I

1118

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

1120

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

1122 Ei(x) + 2*I*pi

1123

1124 Differentiation is supported:

1125

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

1127 exp(x)/x

1128

1129 The exponential integral is related to many other special functions.

1130 For example:

1131

1132 >>> from sympy import expint, Shi

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

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

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

1136 Chi(x) + Shi(x)

1137

1138 See Also

1139 ========

1140

1141 expint: Generalised exponential integral.

1142 E1: Special case of the generalised exponential integral.

1143 li: Logarithmic integral.

1144 Li: Offset logarithmic integral.

1145 Si: Sine integral.

1146 Ci: Cosine integral.

1147 Shi: Hyperbolic sine integral.

1148 Chi: Hyperbolic cosine integral.

1149 uppergamma: Upper incomplete gamma function.

1150

1151 References

1152 ==========

1153

1154 .. [1] https://dlmf.nist.gov/6.6

1155 .. [2] https://en.wikipedia.org/wiki/Exponential_integral

1156 .. [3] Abramowitz & Stegun, section 5: https://web.archive.org/web/20201128173312/http://people.math.sfu.ca/~cbm/aands/page_228.htm

1157

1158 """

1159

1160

1161 @classmethod

1162 def eval(cls, z):

1163 if z.is_zero:

1164 return S.NegativeInfinity

1165 elif z is S.Infinity:

1166 return S.Infinity

1167 elif z is S.NegativeInfinity:

1168 return S.Zero

1169

1170 if z.is_zero:

1171 return S.NegativeInfinity

1172

1173 nz, n = z.extract_branch_factor()

1174 if n:

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

1176

1177 def fdiff(self, argindex=1):

1178 arg = unpolarify(self.args[0])

1179 if argindex == 1:

1180 return exp(arg)/arg

1181 else:

1182 raise ArgumentIndexError(self, argindex)

1183

1184 def _eval_evalf(self, prec):

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

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

1187 return Function._eval_evalf(self, prec)

1188

1189 def _eval_rewrite_as_uppergamma(self, z, **kwargs):

1190 from sympy.functions.special.gamma_functions import uppergamma

1191 # XXX this does not currently work usefully because uppergamma

1192 # immediately turns into expint

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

1194

1195 def _eval_rewrite_as_expint(self, z, **kwargs):

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

1197

1198 def _eval_rewrite_as_li(self, z, **kwargs):

1199 if isinstance(z, log):

1200 return li(z.args[0])

1201 # TODO:

1202 # Actually it only holds that:

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

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

1205 return li(exp(z))

1206

1207 def _eval_rewrite_as_Si(self, z, **kwargs):

1208 if z.is_negative:

1209 return Shi(z) + Chi(z) - I*pi

1210 else:

1211 return Shi(z) + Chi(z)

1212 _eval_rewrite_as_Ci = _eval_rewrite_as_Si

1213 _eval_rewrite_as_Chi = _eval_rewrite_as_Si

1214 _eval_rewrite_as_Shi = _eval_rewrite_as_Si

1215

1216 def _eval_rewrite_as_tractable(self, z, limitvar=None, **kwargs):

1217 return exp(z) * _eis(z)

1218

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

1220 from sympy import re

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

1222 arg = self.args[0].as_leading_term(x, cdir=cdir)

1223 cdir = arg.dir(x, cdir)

1224 if x0.is_zero:

1225 c, e = arg.as_coeff_exponent(x)

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

1227 return log(c) + e*logx + EulerGamma - (

1228 I*pi if re(cdir).is_negative else S.Zero)

1229 return super()._eval_as_leading_term(x, logx=logx, cdir=cdir)

1230

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

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

1233 if x0.is_zero:

1234 f = self._eval_rewrite_as_Si(*self.args)

1235 return f._eval_nseries(x, n, logx)

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

1237

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

1239 from sympy.series.order import Order

1240 point = args0[0]

1241

1242 if point is S.Infinity:

1243 z = self.args[0]

1244 s = [factorial(k) / (z)**k for k in range(n)] + \

1245 [Order(1/z**n, x)]

1246 return (exp(z)/z) * Add(*s)

1247

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

1249

1250

1251class expint(Function):

1252 r"""

1253 Generalized exponential integral.

1254

1255 Explanation

1256 ===========

1257

1258 This function is defined as

1259

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

1261

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

1263 (``uppergamma``).

1264

1265 Hence for $z$ with positive real part we have

1266

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

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

1269

1270 which explains the name.

1271

1272 The representation as an incomplete gamma function provides an analytic

1273 continuation for $\operatorname{E}_\nu(z)$. If $\nu$ is a

1274 non-positive integer, the exponential integral is thus an unbranched

1275 function of $z$, otherwise there is a branch point at the origin.

1276 Refer to the incomplete gamma function documentation for details of the

1277 branching behavior.

1278

1279 Examples

1280 ========

1281

1282 >>> from sympy import expint, S

1283 >>> from sympy.abc import nu, z

1284

1285 Differentiation is supported. Differentiation with respect to $z$ further

1286 explains the name: for integral orders, the exponential integral is an

1287 iterated integral of the exponential function.

1288

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

1290 -expint(nu - 1, z)

1291

1292 Differentiation with respect to $\nu$ has no classical expression:

1293

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

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

1296

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

1298 exponential function:

1299

1300 >>> expint(0, z)

1301 exp(-z)/z

1302 >>> expint(-1, z)

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

1304

1305 At half-integers it reduces to error functions:

1306

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

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

1309

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

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

1312

1313 >>> from sympy import expand_func

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

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

1316

1317 The generalised exponential integral is essentially equivalent to the

1318 incomplete gamma function:

1319

1320 >>> from sympy import uppergamma

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

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

1323

1324 As such it is branched at the origin:

1325

1326 >>> from sympy import exp_polar, pi, I

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

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

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

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

1331

1332 See Also

1333 ========

1334

1335 Ei: Another related function called exponential integral.

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

1337 li: Logarithmic integral.

1338 Li: Offset logarithmic integral.

1339 Si: Sine integral.

1340 Ci: Cosine integral.

1341 Shi: Hyperbolic sine integral.

1342 Chi: Hyperbolic cosine integral.

1343 uppergamma

1344

1345 References

1346 ==========

1347

1348 .. [1] https://dlmf.nist.gov/8.19

1349 .. [2] https://functions.wolfram.com/GammaBetaErf/ExpIntegralE/

1350 .. [3] https://en.wikipedia.org/wiki/Exponential_integral

1351

1352 """

1353

1354

1355 @classmethod

1356 def eval(cls, nu, z):

1357 from sympy.functions.special.gamma_functions import (gamma, uppergamma)

1358 nu2 = unpolarify(nu)

1359 if nu != nu2:

1360 return expint(nu2, z)

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

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

1363

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

1365 # explained in lowergamma.eval().

1366 z, n = z.extract_branch_factor()

1367 if n is S.Zero:

1368 return

1369 if nu.is_integer:

1370 if not nu > 0:

1371 return

1372 return expint(nu, z) \

1373 - 2*pi*I*n*S.NegativeOne**(nu - 1)/factorial(nu - 1)*unpolarify(z)**(nu - 1)

1374 else:

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

1376

1377 def fdiff(self, argindex):

1378 nu, z = self.args

1379 if argindex == 1:

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

1381 elif argindex == 2:

1382 return -expint(nu - 1, z)

1383 else:

1384 raise ArgumentIndexError(self, argindex)

1385

1386 def _eval_rewrite_as_uppergamma(self, nu, z, **kwargs):

1387 from sympy.functions.special.gamma_functions import uppergamma

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

1389

1390 def _eval_rewrite_as_Ei(self, nu, z, **kwargs):

1391 if nu == 1:

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

1393 elif nu.is_Integer and nu > 1:

1394 # DLMF, 8.19.7

1395 x = -unpolarify(z)

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

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

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

1399 else:

1400 return self

1401

1402 def _eval_expand_func(self, **hints):

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

1404

1405 def _eval_rewrite_as_Si(self, nu, z, **kwargs):

1406 if nu != 1:

1407 return self

1408 return Shi(z) - Chi(z)

1409 _eval_rewrite_as_Ci = _eval_rewrite_as_Si

1410 _eval_rewrite_as_Chi = _eval_rewrite_as_Si

1411 _eval_rewrite_as_Shi = _eval_rewrite_as_Si

1412

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

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

1415 nu = self.args[0]

1416 if nu == 1:

1417 f = self._eval_rewrite_as_Si(*self.args)

1418 return f._eval_nseries(x, n, logx)

1419 elif nu.is_Integer and nu > 1:

1420 f = self._eval_rewrite_as_Ei(*self.args)

1421 return f._eval_nseries(x, n, logx)

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

1423

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

1425 from sympy.series.order import Order

1426 point = args0[1]

1427 nu = self.args[0]

1428

1429 if point is S.Infinity:

1430 z = self.args[1]

1431 s = [S.NegativeOne**k * RisingFactorial(nu, k) / z**k for k in range(n)] + [Order(1/z**n, x)]

1432 return (exp(-z)/z) * Add(*s)

1433

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

1435

1436

1437def E1(z):

1438 """

1439 Classical case of the generalized exponential integral.

1440

1441 Explanation

1442 ===========

1443

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

1445

1446 Examples

1447 ========

1448

1449 >>> from sympy import E1

1450 >>> E1(0)

1451 expint(1, 0)

1452

1453 >>> E1(5)

1454 expint(1, 5)

1455

1456 See Also

1457 ========

1458

1459 Ei: Exponential integral.

1460 expint: Generalised exponential integral.

1461 li: Logarithmic integral.

1462 Li: Offset logarithmic integral.

1463 Si: Sine integral.

1464 Ci: Cosine integral.

1465 Shi: Hyperbolic sine integral.

1466 Chi: Hyperbolic cosine integral.

1467

1468 """

1469 return expint(1, z)

1470

1471

1472class li(Function):

1473 r"""

1474 The classical logarithmic integral.

1475

1476 Explanation

1477 ===========

1478

1479 For use in SymPy, this function is defined as

1480

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

1482

1483 Examples

1484 ========

1485

1486 >>> from sympy import I, oo, li

1487 >>> from sympy.abc import z

1488

1489 Several special values are known:

1490

1491 >>> li(0)

1492 0

1493 >>> li(1)

1494 -oo

1495 >>> li(oo)

1496 oo

1497

1498 Differentiation with respect to $z$ is supported:

1499

1500 >>> from sympy import diff

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

1502 1/log(z)

1503

1504 Defining the ``li`` function via an integral:

1505 >>> from sympy import integrate

1506 >>> integrate(li(z))

1507 z*li(z) - Ei(2*log(z))

1508

1509 >>> integrate(li(z),z)

1510 z*li(z) - Ei(2*log(z))

1511

1512

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

1514

1515 >>> from sympy import Ei

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

1517 Ei(log(z))

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

1519 1/log(z)

1520

1521 We can numerically evaluate the logarithmic integral to arbitrary precision

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

1523

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

1525 1.04516378011749278484458888919

1526

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

1528 1.0652795784357498247001125598 + 3.08346052231061726610939702133*I

1529

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

1531

1532 >>> from mpmath import findroot

1533 >>> findroot(li, 2)

1534 1.45136923488338

1535

1536 Further transformations include rewriting ``li`` in terms of

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

1538

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

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

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

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

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

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

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

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

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

1548

1549 See Also

1550 ========

1551

1552 Li: Offset logarithmic integral.

1553 Ei: Exponential integral.

1554 expint: Generalised exponential integral.

1555 E1: Special case of the generalised exponential integral.

1556 Si: Sine integral.

1557 Ci: Cosine integral.

1558 Shi: Hyperbolic sine integral.

1559 Chi: Hyperbolic cosine integral.

1560

1561 References

1562 ==========

1563

1564 .. [1] https://en.wikipedia.org/wiki/Logarithmic_integral

1565 .. [2] https://mathworld.wolfram.com/LogarithmicIntegral.html

1566 .. [3] https://dlmf.nist.gov/6

1567 .. [4] https://mathworld.wolfram.com/SoldnersConstant.html

1568

1569 """

1570

1571

1572 @classmethod

1573 def eval(cls, z):

1574 if z.is_zero:

1575 return S.Zero

1576 elif z is S.One:

1577 return S.NegativeInfinity

1578 elif z is S.Infinity:

1579 return S.Infinity

1580 if z.is_zero:

1581 return S.Zero

1582

1583 def fdiff(self, argindex=1):

1584 arg = self.args[0]

1585 if argindex == 1:

1586 return S.One / log(arg)

1587 else:

1588 raise ArgumentIndexError(self, argindex)

1589

1590 def _eval_conjugate(self):

1591 z = self.args[0]

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

1593 if not z.is_extended_negative:

1594 return self.func(z.conjugate())

1595

1596 def _eval_rewrite_as_Li(self, z, **kwargs):

1597 return Li(z) + li(2)

1598

1599 def _eval_rewrite_as_Ei(self, z, **kwargs):

1600 return Ei(log(z))

1601

1602 def _eval_rewrite_as_uppergamma(self, z, **kwargs):

1603 from sympy.functions.special.gamma_functions import uppergamma

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

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

1606

1607 def _eval_rewrite_as_Si(self, z, **kwargs):

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

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

1610

1611 _eval_rewrite_as_Ci = _eval_rewrite_as_Si

1612

1613 def _eval_rewrite_as_Shi(self, z, **kwargs):

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

1615

1616 _eval_rewrite_as_Chi = _eval_rewrite_as_Shi

1617

1618 def _eval_rewrite_as_hyper(self, z, **kwargs):

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

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

1621

1622 def _eval_rewrite_as_meijerg(self, z, **kwargs):

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

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

1625

1626 def _eval_rewrite_as_tractable(self, z, limitvar=None, **kwargs):

1627 return z * _eis(log(z))

1628

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

1630 z = self.args[0]

1631 s = [(log(z))**k / (factorial(k) * k) for k in range(1, n)]

1632 return EulerGamma + log(log(z)) + Add(*s)

1633

1634 def _eval_is_zero(self):

1635 z = self.args[0]

1636 if z.is_zero:

1637 return True

1638

1639class Li(Function):

1640 r"""

1641 The offset logarithmic integral.

1642

1643 Explanation

1644 ===========

1645

1646 For use in SymPy, this function is defined as

1647

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

1649

1650 Examples

1651 ========

1652

1653 >>> from sympy import Li

1654 >>> from sympy.abc import z

1655

1656 The following special value is known:

1657

1658 >>> Li(2)

1659 0

1660

1661 Differentiation with respect to $z$ is supported:

1662

1663 >>> from sympy import diff

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

1665 1/log(z)

1666

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

1668

1669 >>> from sympy import li

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

1671 li(z) - li(2)

1672

1673 We can numerically evaluate the logarithmic integral to arbitrary precision

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

1675

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

1677 0

1678

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

1680 1.92242131492155809316615998938

1681

1682 See Also

1683 ========

1684

1685 li: Logarithmic integral.

1686 Ei: Exponential integral.

1687 expint: Generalised exponential integral.

1688 E1: Special case of the generalised exponential integral.

1689 Si: Sine integral.

1690 Ci: Cosine integral.

1691 Shi: Hyperbolic sine integral.

1692 Chi: Hyperbolic cosine integral.

1693

1694 References

1695 ==========

1696

1697 .. [1] https://en.wikipedia.org/wiki/Logarithmic_integral

1698 .. [2] https://mathworld.wolfram.com/LogarithmicIntegral.html

1699 .. [3] https://dlmf.nist.gov/6

1700

1701 """

1702

1703

1704 @classmethod

1705 def eval(cls, z):

1706 if z is S.Infinity:

1707 return S.Infinity

1708 elif z == S(2):

1709 return S.Zero

1710

1711 def fdiff(self, argindex=1):

1712 arg = self.args[0]

1713 if argindex == 1:

1714 return S.One / log(arg)

1715 else:

1716 raise ArgumentIndexError(self, argindex)

1717

1718 def _eval_evalf(self, prec):

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

1720

1721 def _eval_rewrite_as_li(self, z, **kwargs):

1722 return li(z) - li(2)

1723

1724 def _eval_rewrite_as_tractable(self, z, limitvar=None, **kwargs):

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

1726

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

1728 f = self._eval_rewrite_as_li(*self.args)

1729 return f._eval_nseries(x, n, logx)

1730

1731###############################################################################

1732#################### TRIGONOMETRIC INTEGRALS ##################################

1733###############################################################################

1734

1735class TrigonometricIntegral(Function):

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

1737

1738

1739 @classmethod

1740 def eval(cls, z):

1741 if z is S.Zero:

1742 return cls._atzero

1743 elif z is S.Infinity:

1744 return cls._atinf()

1745 elif z is S.NegativeInfinity:

1746 return cls._atneginf()

1747

1748 if z.is_zero:

1749 return cls._atzero

1750

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

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

1753 nz = z.extract_multiplicatively(I)

1754 if nz is not None:

1755 return cls._Ifactor(nz, 1)

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

1757 if nz is not None:

1758 return cls._Ifactor(nz, -1)

1759

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

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

1762 nz = z.extract_multiplicatively(-1)

1763 if nz is not None:

1764 return cls._minusfactor(nz)

1765

1766 nz, n = z.extract_branch_factor()

1767 if n == 0 and nz == z:

1768 return

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

1770

1771 def fdiff(self, argindex=1):

1772 arg = unpolarify(self.args[0])

1773 if argindex == 1:

1774 return self._trigfunc(arg)/arg

1775 else:

1776 raise ArgumentIndexError(self, argindex)

1777

1778 def _eval_rewrite_as_Ei(self, z, **kwargs):

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

1780

1781 def _eval_rewrite_as_uppergamma(self, z, **kwargs):

1782 from sympy.functions.special.gamma_functions import uppergamma

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

1784

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

1786 # NOTE this is fairly inefficient

1787 n += 1

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

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

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

1791 if self._trigfunc(0) != 0:

1792 baseseries -= 1

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

1794 if self._trigfunc(0) != 0:

1795 baseseries += EulerGamma + log(x)

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

1797

1798

1799class Si(TrigonometricIntegral):

1800 r"""

1801 Sine integral.

1802

1803 Explanation

1804 ===========

1805

1806 This function is defined by

1807

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

1809

1810 It is an entire function.

1811

1812 Examples

1813 ========

1814

1815 >>> from sympy import Si

1816 >>> from sympy.abc import z

1817

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

1819

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

1821 sin(z)/z

1822

1823 It is unbranched:

1824

1825 >>> from sympy import exp_polar, I, pi

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

1827 Si(z)

1828

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

1830

1831 >>> Si(I*z)

1832 I*Shi(z)

1833 >>> Si(-z)

1834 -Si(z)

1835

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

1837 that the latter is branched:

1838

1839 >>> from sympy import expint

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

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

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

1843

1844 It can be rewritten in the form of sinc function (by definition):

1845

1846 >>> from sympy import sinc

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

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

1849

1850 See Also

1851 ========

1852

1853 Ci: Cosine integral.

1854 Shi: Hyperbolic sine integral.

1855 Chi: Hyperbolic cosine integral.

1856 Ei: Exponential integral.

1857 expint: Generalised exponential integral.

1858 sinc: unnormalized sinc function

1859 E1: Special case of the generalised exponential integral.

1860 li: Logarithmic integral.

1861 Li: Offset logarithmic integral.

1862

1863 References

1864 ==========

1865

1866 .. [1] https://en.wikipedia.org/wiki/Trigonometric_integral

1867

1868 """

1869

1870 _trigfunc = sin

1871 _atzero = S.Zero

1872

1873 @classmethod

1874 def _atinf(cls):

1875 return pi*S.Half

1876

1877 @classmethod

1878 def _atneginf(cls):

1879 return -pi*S.Half

1880

1881 @classmethod

1882 def _minusfactor(cls, z):

1883 return -Si(z)

1884

1885 @classmethod

1886 def _Ifactor(cls, z, sign):

1887 return I*Shi(z)*sign

1888

1889 def _eval_rewrite_as_expint(self, z, **kwargs):

1890 # XXX should we polarify z?

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

1892

1893 def _eval_rewrite_as_sinc(self, z, **kwargs):

1894 from sympy.integrals.integrals import Integral

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

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

1897

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

1899 from sympy.series.order import Order

1900 point = args0[0]

1901

1902 # Expansion at oo

1903 if point is S.Infinity:

1904 z = self.args[0]

1905 p = [S.NegativeOne**k * factorial(2*k) / z**(2*k)

1906 for k in range(int((n - 1)/2))] + [Order(1/z**n, x)]

1907 q = [S.NegativeOne**k * factorial(2*k + 1) / z**(2*k + 1)

1908 for k in range(int(n/2) - 1)] + [Order(1/z**n, x)]

1909 return pi/2 - (cos(z)/z)*Add(*p) - (sin(z)/z)*Add(*q)

1910

1911 # All other points are not handled

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

1913

1914 def _eval_is_zero(self):

1915 z = self.args[0]

1916 if z.is_zero:

1917 return True

1918

1919

1920class Ci(TrigonometricIntegral):

1921 r"""

1922 Cosine integral.

1923

1924 Explanation

1925 ===========

1926

1927 This function is defined for positive $x$ by

1928

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

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

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

1932

1933 where $\gamma$ is the Euler-Mascheroni constant.

1934

1935 We have

1936

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

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

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

1940

1941 which holds for all polar $z$ and thus provides an analytic

1942 continuation to the Riemann surface of the logarithm.

1943

1944 The formula also holds as stated

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

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

1947 cut complex plane.

1948

1949 Examples

1950 ========

1951

1952 >>> from sympy import Ci

1953 >>> from sympy.abc import z

1954

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

1956

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

1958 cos(z)/z

1959

1960 It has a logarithmic branch point at the origin:

1961

1962 >>> from sympy import exp_polar, I, pi

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

1964 Ci(z) + 2*I*pi

1965

1966 The cosine integral behaves somewhat like ordinary $\cos$ under

1967 multiplication by $i$:

1968

1969 >>> from sympy import polar_lift

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

1971 Chi(z) + I*pi/2

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

1973 Ci(z) + I*pi

1974

1975 It can also be expressed in terms of exponential integrals:

1976

1977 >>> from sympy import expint

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

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

1980

1981 See Also

1982 ========

1983

1984 Si: Sine integral.

1985 Shi: Hyperbolic sine integral.

1986 Chi: Hyperbolic cosine integral.

1987 Ei: Exponential integral.

1988 expint: Generalised exponential integral.

1989 E1: Special case of the generalised exponential integral.

1990 li: Logarithmic integral.

1991 Li: Offset logarithmic integral.

1992

1993 References

1994 ==========

1995

1996 .. [1] https://en.wikipedia.org/wiki/Trigonometric_integral

1997

1998 """

1999

2000 _trigfunc = cos

2001 _atzero = S.ComplexInfinity

2002

2003 @classmethod

2004 def _atinf(cls):

2005 return S.Zero

2006

2007 @classmethod

2008 def _atneginf(cls):

2009 return I*pi

2010

2011 @classmethod

2012 def _minusfactor(cls, z):

2013 return Ci(z) + I*pi

2014

2015 @classmethod

2016 def _Ifactor(cls, z, sign):

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

2018

2019 def _eval_rewrite_as_expint(self, z, **kwargs):

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

2021

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

2023 arg = self.args[0].as_leading_term(x, logx=logx, cdir=cdir)

2024 arg0 = arg.subs(x, 0)

2025

2026 if arg0 is S.NaN:

2027 arg0 = arg.limit(x, 0, dir='-' if re(cdir).is_negative else '+')

2028 if arg0.is_zero:

2029 c, e = arg.as_coeff_exponent(x)

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

2031 return log(c) + e*logx + EulerGamma

2032 elif arg0.is_finite:

2033 return self.func(arg0)

2034 else:

2035 return self

2036

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

2038 from sympy.series.order import Order

2039 point = args0[0]

2040

2041 # Expansion at oo

2042 if point is S.Infinity:

2043 z = self.args[0]

2044 p = [S.NegativeOne**k * factorial(2*k) / z**(2*k)

2045 for k in range(int((n - 1)/2))] + [Order(1/z**n, x)]

2046 q = [S.NegativeOne**k * factorial(2*k + 1) / z**(2*k + 1)

2047 for k in range(int(n/2) - 1)] + [Order(1/z**n, x)]

2048 return (sin(z)/z)*Add(*p) - (cos(z)/z)*Add(*q)

2049

2050 # All other points are not handled

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

2052

2053

2054class Shi(TrigonometricIntegral):

2055 r"""

2056 Sinh integral.

2057

2058 Explanation

2059 ===========

2060

2061 This function is defined by

2062

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

2064

2065 It is an entire function.

2066

2067 Examples

2068 ========

2069

2070 >>> from sympy import Shi

2071 >>> from sympy.abc import z

2072

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

2074

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

2076 sinh(z)/z

2077

2078 It is unbranched:

2079

2080 >>> from sympy import exp_polar, I, pi

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

2082 Shi(z)

2083

2084 The $\sinh$ integral behaves much like ordinary $\sinh$ under

2085 multiplication by $i$:

2086

2087 >>> Shi(I*z)

2088 I*Si(z)

2089 >>> Shi(-z)

2090 -Shi(z)

2091

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

2093 that the latter is branched:

2094

2095 >>> from sympy import expint

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

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

2098

2099 See Also

2100 ========

2101

2102 Si: Sine integral.

2103 Ci: Cosine integral.

2104 Chi: Hyperbolic cosine integral.

2105 Ei: Exponential integral.

2106 expint: Generalised exponential integral.

2107 E1: Special case of the generalised exponential integral.

2108 li: Logarithmic integral.

2109 Li: Offset logarithmic integral.

2110

2111 References

2112 ==========

2113

2114 .. [1] https://en.wikipedia.org/wiki/Trigonometric_integral

2115

2116 """

2117

2118 _trigfunc = sinh

2119 _atzero = S.Zero

2120

2121 @classmethod

2122 def _atinf(cls):

2123 return S.Infinity

2124

2125 @classmethod

2126 def _atneginf(cls):

2127 return S.NegativeInfinity

2128

2129 @classmethod

2130 def _minusfactor(cls, z):

2131 return -Shi(z)

2132

2133 @classmethod

2134 def _Ifactor(cls, z, sign):

2135 return I*Si(z)*sign

2136

2137 def _eval_rewrite_as_expint(self, z, **kwargs):

2138 # XXX should we polarify z?

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

2140

2141 def _eval_is_zero(self):

2142 z = self.args[0]

2143 if z.is_zero:

2144 return True

2145

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

2147 arg = self.args[0].as_leading_term(x)

2148 arg0 = arg.subs(x, 0)

2149

2150 if arg0 is S.NaN:

2151 arg0 = arg.limit(x, 0, dir='-' if re(cdir).is_negative else '+')

2152 if arg0.is_zero:

2153 return arg

2154 elif not arg0.is_infinite:

2155 return self.func(arg0)

2156 else:

2157 return self

2158

2159

2160class Chi(TrigonometricIntegral):

2161 r"""

2162 Cosh integral.

2163

2164 Explanation

2165 ===========

2166

2167 This function is defined for positive $x$ by

2168

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

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

2171

2172 where $\gamma$ is the Euler-Mascheroni constant.

2173

2174 We have

2175

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

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

2178

2179 which holds for all polar $z$ and thus provides an analytic

2180 continuation to the Riemann surface of the logarithm.

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

2182 cut complex plane.

2183

2184 Examples

2185 ========

2186

2187 >>> from sympy import Chi

2188 >>> from sympy.abc import z

2189

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

2191

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

2193 cosh(z)/z

2194

2195 It has a logarithmic branch point at the origin:

2196

2197 >>> from sympy import exp_polar, I, pi

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

2199 Chi(z) + 2*I*pi

2200

2201 The $\cosh$ integral behaves somewhat like ordinary $\cosh$ under

2202 multiplication by $i$:

2203

2204 >>> from sympy import polar_lift

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

2206 Ci(z) + I*pi/2

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

2208 Chi(z) + I*pi

2209

2210 It can also be expressed in terms of exponential integrals:

2211

2212 >>> from sympy import expint

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

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

2215

2216 See Also

2217 ========

2218

2219 Si: Sine integral.

2220 Ci: Cosine integral.

2221 Shi: Hyperbolic sine integral.

2222 Ei: Exponential integral.

2223 expint: Generalised exponential integral.

2224 E1: Special case of the generalised exponential integral.

2225 li: Logarithmic integral.

2226 Li: Offset logarithmic integral.

2227

2228 References

2229 ==========

2230

2231 .. [1] https://en.wikipedia.org/wiki/Trigonometric_integral

2232

2233 """

2234

2235 _trigfunc = cosh

2236 _atzero = S.ComplexInfinity

2237

2238 @classmethod

2239 def _atinf(cls):

2240 return S.Infinity

2241

2242 @classmethod

2243 def _atneginf(cls):

2244 return S.Infinity

2245

2246 @classmethod

2247 def _minusfactor(cls, z):

2248 return Chi(z) + I*pi

2249

2250 @classmethod

2251 def _Ifactor(cls, z, sign):

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

2253

2254 def _eval_rewrite_as_expint(self, z, **kwargs):

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

2256

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

2258 arg = self.args[0].as_leading_term(x, logx=logx, cdir=cdir)

2259 arg0 = arg.subs(x, 0)

2260

2261 if arg0 is S.NaN:

2262 arg0 = arg.limit(x, 0, dir='-' if re(cdir).is_negative else '+')

2263 if arg0.is_zero:

2264 c, e = arg.as_coeff_exponent(x)

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

2266 return log(c) + e*logx + EulerGamma

2267 elif arg0.is_finite:

2268 return self.func(arg0)

2269 else:

2270 return self

2271

2272

2273###############################################################################

2274#################### FRESNEL INTEGRALS ########################################

2275###############################################################################

2276

2277class FresnelIntegral(Function):

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

2279

2280 unbranched = True

2281

2282 @classmethod

2283 def eval(cls, z):

2284 # Values at positive infinities signs

2285 # if any were extracted automatically

2286 if z is S.Infinity:

2287 return S.Half

2288

2289 # Value at zero

2290 if z.is_zero:

2291 return S.Zero

2292

2293 # Try to pull out factors of -1 and I

2294 prefact = S.One

2295 newarg = z

2296 changed = False

2297

2298 nz = newarg.extract_multiplicatively(-1)

2299 if nz is not None:

2300 prefact = -prefact

2301 newarg = nz

2302 changed = True

2303

2304 nz = newarg.extract_multiplicatively(I)

2305 if nz is not None:

2306 prefact = cls._sign*I*prefact

2307 newarg = nz

2308 changed = True

2309

2310 if changed:

2311 return prefact*cls(newarg)

2312

2313 def fdiff(self, argindex=1):

2314 if argindex == 1:

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

2316 else:

2317 raise ArgumentIndexError(self, argindex)

2318

2319 def _eval_is_extended_real(self):

2320 return self.args[0].is_extended_real

2321

2322 _eval_is_finite = _eval_is_extended_real

2323

2324 def _eval_is_zero(self):

2325 return self.args[0].is_zero

2326

2327 def _eval_conjugate(self):

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

2329

2330 as_real_imag = real_to_real_as_real_imag

2331

2332

2333class fresnels(FresnelIntegral):

2334 r"""

2335 Fresnel integral S.

2336

2337 Explanation

2338 ===========

2339

2340 This function is defined by

2341

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

2343

2344 It is an entire function.

2345

2346 Examples

2347 ========

2348

2349 >>> from sympy import I, oo, fresnels

2350 >>> from sympy.abc import z

2351

2352 Several special values are known:

2353

2354 >>> fresnels(0)

2355 0

2356 >>> fresnels(oo)

2357 1/2

2358 >>> fresnels(-oo)

2359 -1/2

2360 >>> fresnels(I*oo)

2361 -I/2

2362 >>> fresnels(-I*oo)

2363 I/2

2364

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

2366

2367 >>> fresnels(-z)

2368 -fresnels(z)

2369 >>> fresnels(I*z)

2370 -I*fresnels(z)

2371

2372 The Fresnel S integral obeys the mirror symmetry

2373 $\overline{S(z)} = S(\bar{z})$:

2374

2375 >>> from sympy import conjugate

2376 >>> conjugate(fresnels(z))

2377 fresnels(conjugate(z))

2378

2379 Differentiation with respect to $z$ is supported:

2380

2381 >>> from sympy import diff

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

2383 sin(pi*z**2/2)

2384

2385 Defining the Fresnel functions via an integral:

2386

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

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

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

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

2391 fresnels(z)

2392

2393 We can numerically evaluate the Fresnel integral to arbitrary precision

2394 on the whole complex plane:

2395

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

2397 0.343415678363698242195300815958

2398

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

2400 0.343415678363698242195300815958*I

2401

2402 See Also

2403 ========

2404

2405 fresnelc: Fresnel cosine integral.

2406

2407 References

2408 ==========

2409

2410 .. [1] https://en.wikipedia.org/wiki/Fresnel_integral

2411 .. [2] https://dlmf.nist.gov/7

2412 .. [3] https://mathworld.wolfram.com/FresnelIntegrals.html

2413 .. [4] https://functions.wolfram.com/GammaBetaErf/FresnelS

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

2415 by John W. Wrench Jr. and Vicki Alley

2416

2417 """

2418 _trigfunc = sin

2419 _sign = -S.One

2420

2421 @staticmethod

2422 @cacheit

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

2424 if n < 0:

2425 return S.Zero

2426 else:

2427 x = sympify(x)

2428 if len(previous_terms) > 1:

2429 p = previous_terms[-1]

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

2431 else:

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

2433

2434 def _eval_rewrite_as_erf(self, z, **kwargs):

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

2436

2437 def _eval_rewrite_as_hyper(self, z, **kwargs):

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

2439

2440 def _eval_rewrite_as_meijerg(self, z, **kwargs):

2441 return (pi*z**Rational(9, 4) / (sqrt(2)*(z**2)**Rational(3, 4)*(-z)**Rational(3, 4))

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

2443

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

2445 from sympy.series.order import Order

2446 arg = self.args[0].as_leading_term(x, logx=logx, cdir=cdir)

2447 arg0 = arg.subs(x, 0)

2448

2449 if arg0 is S.ComplexInfinity:

2450 arg0 = arg.limit(x, 0, dir='-' if re(cdir).is_negative else '+')

2451 if arg0.is_zero:

2452 return pi*arg**3/6

2453 elif arg0 in [S.Infinity, S.NegativeInfinity]:

2454 s = 1 if arg0 is S.Infinity else -1

2455 return s*S.Half + Order(x, x)

2456 else:

2457 return self.func(arg0)

2458

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

2460 from sympy.series.order import Order

2461 point = args0[0]

2462

2463 # Expansion at oo and -oo

2464 if point in [S.Infinity, -S.Infinity]:

2465 z = self.args[0]

2466

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

2468 # as only real infinities are dealt with, sin and cos are O(1)

2469 p = [S.NegativeOne**k * factorial(4*k + 1) /

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

2471 for k in range(0, n) if 4*k + 3 < n]

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

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

2474 for k in range(1, n) if 4*k + 1 < n]

2475

2476 p = [-sqrt(2/pi)*t for t in p]

2477 q = [-sqrt(2/pi)*t for t in q]

2478 s = 1 if point is S.Infinity else -1

2479 # The expansion at oo is 1/2 + some odd powers of z

2480 # To get the expansion at -oo, replace z by -z and flip the sign

2481 # The result -1/2 + the same odd powers of z as before.

2482 return s*S.Half + (sin(z**2)*Add(*p) + cos(z**2)*Add(*q)

2483 ).subs(x, sqrt(2/pi)*x) + Order(1/z**n, x)

2484

2485 # All other points are not handled

2486 return super()._eval_aseries(n, args0, x, logx)

2487

2488

2489class fresnelc(FresnelIntegral):

2490 r"""

2491 Fresnel integral C.

2492

2493 Explanation

2494 ===========

2495

2496 This function is defined by

2497

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

2499

2500 It is an entire function.

2501

2502 Examples

2503 ========

2504

2505 >>> from sympy import I, oo, fresnelc

2506 >>> from sympy.abc import z

2507

2508 Several special values are known:

2509

2510 >>> fresnelc(0)

2511 0

2512 >>> fresnelc(oo)

2513 1/2

2514 >>> fresnelc(-oo)

2515 -1/2

2516 >>> fresnelc(I*oo)

2517 I/2

2518 >>> fresnelc(-I*oo)

2519 -I/2

2520

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

2522

2523 >>> fresnelc(-z)

2524 -fresnelc(z)

2525 >>> fresnelc(I*z)

2526 I*fresnelc(z)

2527

2528 The Fresnel C integral obeys the mirror symmetry

2529 $\overline{C(z)} = C(\bar{z})$:

2530

2531 >>> from sympy import conjugate

2532 >>> conjugate(fresnelc(z))

2533 fresnelc(conjugate(z))

2534

2535 Differentiation with respect to $z$ is supported:

2536

2537 >>> from sympy import diff

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

2539 cos(pi*z**2/2)

2540

2541 Defining the Fresnel functions via an integral:

2542

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

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

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

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

2547 fresnelc(z)

2548

2549 We can numerically evaluate the Fresnel integral to arbitrary precision

2550 on the whole complex plane:

2551

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

2553 0.488253406075340754500223503357

2554

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

2556 -0.488253406075340754500223503357*I

2557

2558 See Also

2559 ========

2560

2561 fresnels: Fresnel sine integral.

2562

2563 References

2564 ==========

2565

2566 .. [1] https://en.wikipedia.org/wiki/Fresnel_integral

2567 .. [2] https://dlmf.nist.gov/7

2568 .. [3] https://mathworld.wolfram.com/FresnelIntegrals.html

2569 .. [4] https://functions.wolfram.com/GammaBetaErf/FresnelC

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

2571 by John W. Wrench Jr. and Vicki Alley

2572

2573 """

2574 _trigfunc = cos

2575 _sign = S.One

2576

2577 @staticmethod

2578 @cacheit

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

2580 if n < 0:

2581 return S.Zero

2582 else:

2583 x = sympify(x)

2584 if len(previous_terms) > 1:

2585 p = previous_terms[-1]

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

2587 else:

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

2589

2590 def _eval_rewrite_as_erf(self, z, **kwargs):

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

2592

2593 def _eval_rewrite_as_hyper(self, z, **kwargs):

2594 return z * hyper([Rational(1, 4)], [S.Half, Rational(5, 4)], -pi**2*z**4/16)

2595

2596 def _eval_rewrite_as_meijerg(self, z, **kwargs):

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

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

2599

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

2601 from sympy.series.order import Order

2602 arg = self.args[0].as_leading_term(x, logx=logx, cdir=cdir)

2603 arg0 = arg.subs(x, 0)

2604

2605 if arg0 is S.ComplexInfinity:

2606 arg0 = arg.limit(x, 0, dir='-' if re(cdir).is_negative else '+')

2607 if arg0.is_zero:

2608 return arg

2609 elif arg0 in [S.Infinity, S.NegativeInfinity]:

2610 s = 1 if arg0 is S.Infinity else -1

2611 return s*S.Half + Order(x, x)

2612 else:

2613 return self.func(arg0)

2614

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

2616 from sympy.series.order import Order

2617 point = args0[0]

2618

2619 # Expansion at oo

2620 if point in [S.Infinity, -S.Infinity]:

2621 z = self.args[0]

2622

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

2624 # as only real infinities are dealt with, sin and cos are O(1)

2625 p = [S.NegativeOne**k * factorial(4*k + 1) /

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

2627 for k in range(n) if 4*k + 3 < n]

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

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

2630 for k in range(1, n) if 4*k + 1 < n]

2631

2632 p = [-sqrt(2/pi)*t for t in p]

2633 q = [ sqrt(2/pi)*t for t in q]

2634 s = 1 if point is S.Infinity else -1

2635 # The expansion at oo is 1/2 + some odd powers of z

2636 # To get the expansion at -oo, replace z by -z and flip the sign

2637 # The result -1/2 + the same odd powers of z as before.

2638 return s*S.Half + (cos(z**2)*Add(*p) + sin(z**2)*Add(*q)

2639 ).subs(x, sqrt(2/pi)*x) + Order(1/z**n, x)

2640

2641 # All other points are not handled

2642 return super()._eval_aseries(n, args0, x, logx)

2643

2644

2645###############################################################################

2646#################### HELPER FUNCTIONS #########################################

2647###############################################################################

2648

2649

2650class _erfs(Function):

2651 """

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

2653 tractable for the Gruntz algorithm.

2654

2655 """

2656 @classmethod

2657 def eval(cls, arg):

2658 if arg.is_zero:

2659 return S.One

2660

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

2662 from sympy.series.order import Order

2663 point = args0[0]

2664

2665 # Expansion at oo

2666 if point is S.Infinity:

2667 z = self.args[0]

2668 l = [1/sqrt(pi) * factorial(2*k)*(-S(

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

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

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

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

2673

2674 # Expansion at I*oo

2675 t = point.extract_multiplicatively(I)

2676 if t is S.Infinity:

2677 z = self.args[0]

2678 # TODO: is the series really correct?

2679 l = [1/sqrt(pi) * factorial(2*k)*(-S(

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

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

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

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

2684

2685 # All other points are not handled

2686 return super()._eval_aseries(n, args0, x, logx)

2687

2688 def fdiff(self, argindex=1):

2689 if argindex == 1:

2690 z = self.args[0]

2691 return -2/sqrt(pi) + 2*z*_erfs(z)

2692 else:

2693 raise ArgumentIndexError(self, argindex)

2694

2695 def _eval_rewrite_as_intractable(self, z, **kwargs):

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

2697

2698

2699class _eis(Function):

2700 """

2701 Helper function to make the $\\mathrm{Ei}(z)$ and $\\mathrm{li}(z)$

2702 functions tractable for the Gruntz algorithm.

2703

2704 """

2705

2706

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

2708 from sympy.series.order import Order

2709 if args0[0] != S.Infinity:

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

2711

2712 z = self.args[0]

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

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

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

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

2717

2718

2719 def fdiff(self, argindex=1):

2720 if argindex == 1:

2721 z = self.args[0]

2722 return S.One / z - _eis(z)

2723 else:

2724 raise ArgumentIndexError(self, argindex)

2725

2726 def _eval_rewrite_as_intractable(self, z, **kwargs):

2727 return exp(-z)*Ei(z)

2728

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

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

2731 if x0.is_zero:

2732 f = self._eval_rewrite_as_intractable(*self.args)

2733 return f._eval_as_leading_term(x, logx=logx, cdir=cdir)

2734 return super()._eval_as_leading_term(x, logx=logx, cdir=cdir)

2735

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

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

2738 if x0.is_zero:

2739 f = self._eval_rewrite_as_intractable(*self.args)

2740 return f._eval_nseries(x, n, logx)

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

Coverage for /usr/lib/python3/dist-packages/sympy/functions/special/error_functions.py: 29%

992 statements