Coverage for /usr/lib/python3/dist-packages/sympy/series/formal.py: 17%

1"""Formal Power Series"""

3from collections import defaultdict

5from sympy.core.numbers import (nan, oo, zoo)

6from sympy.core.add import Add

7from sympy.core.expr import Expr

8from sympy.core.function import Derivative, Function, expand

9from sympy.core.mul import Mul

10from sympy.core.numbers import Rational

11from sympy.core.relational import Eq

12from sympy.sets.sets import Interval

13from sympy.core.singleton import S

14from sympy.core.symbol import Wild, Dummy, symbols, Symbol

15from sympy.core.sympify import sympify

16from sympy.discrete.convolutions import convolution

17from sympy.functions.combinatorial.factorials import binomial, factorial, rf

18from sympy.functions.combinatorial.numbers import bell

19from sympy.functions.elementary.integers import floor, frac, ceiling

20from sympy.functions.elementary.miscellaneous import Min, Max

21from sympy.functions.elementary.piecewise import Piecewise

22from sympy.series.limits import Limit

23from sympy.series.order import Order

24from sympy.series.sequences import sequence

25from sympy.series.series_class import SeriesBase

26from sympy.utilities.iterables import iterable

30def rational_algorithm(f, x, k, order=4, full=False):

31 """

32 Rational algorithm for computing

33 formula of coefficients of Formal Power Series

34 of a function.

36 Explanation

37 ===========

39 Applicable when f(x) or some derivative of f(x)

40 is a rational function in x.

42 :func:`rational_algorithm` uses :func:`~.apart` function for partial fraction

43 decomposition. :func:`~.apart` by default uses 'undetermined coefficients

44 method'. By setting ``full=True``, 'Bronstein's algorithm' can be used

45 instead.

47 Looks for derivative of a function up to 4'th order (by default).

48 This can be overridden using order option.

50 Parameters

51 ==========

53 x : Symbol

54 order : int, optional

55 Order of the derivative of ``f``, Default is 4.

56 full : bool

58 Returns

59 =======

61 formula : Expr

62 ind : Expr

63 Independent terms.

64 order : int

65 full : bool

67 Examples

68 ========

70 >>> from sympy import log, atan

71 >>> from sympy.series.formal import rational_algorithm as ra

72 >>> from sympy.abc import x, k

74 >>> ra(1 / (1 - x), x, k)

75 (1, 0, 0)

76 >>> ra(log(1 + x), x, k)

77 (-1/((-1)**k*k), 0, 1)

79 >>> ra(atan(x), x, k, full=True)

80 ((-I/(2*(-I)**k) + I/(2*I**k))/k, 0, 1)

82 Notes

83 =====

85 By setting ``full=True``, range of admissible functions to be solved using

86 ``rational_algorithm`` can be increased. This option should be used

87 carefully as it can significantly slow down the computation as ``doit`` is

88 performed on the :class:`~.RootSum` object returned by the :func:`~.apart`

89 function. Use ``full=False`` whenever possible.

91 See Also

92 ========

94 sympy.polys.partfrac.apart

96 References

97 ==========

99 .. [1] Formal Power Series - Dominik Gruntz, Wolfram Koepf

100 .. [2] Power Series in Computer Algebra - Wolfram Koepf

101

102 """

103 from sympy.polys import RootSum, apart

104 from sympy.integrals import integrate

105

106 diff = f

107 ds = [] # list of diff

108

109 for i in range(order + 1):

110 if i:

111 diff = diff.diff(x)

112

113 if diff.is_rational_function(x):

114 coeff, sep = S.Zero, S.Zero

115

116 terms = apart(diff, x, full=full)

117 if terms.has(RootSum):

118 terms = terms.doit()

119

120 for t in Add.make_args(terms):

121 num, den = t.as_numer_denom()

122 if not den.has(x):

123 sep += t

124 else:

125 if isinstance(den, Mul):

126 # m*(n*x - a)**j -> (n*x - a)**j

127 ind = den.as_independent(x)

128 den = ind[1]

129 num /= ind[0]

130

131 # (n*x - a)**j -> (x - b)

132 den, j = den.as_base_exp()

133 a, xterm = den.as_coeff_add(x)

134

135 # term -> m/x**n

136 if not a:

137 sep += t

138 continue

139

140 xc = xterm[0].coeff(x)

141 a /= -xc

142 num /= xc**j

143

144 ak = ((-1)**j * num *

145 binomial(j + k - 1, k).rewrite(factorial) /

146 a**(j + k))

147 coeff += ak

148

149 # Hacky, better way?

150 if coeff.is_zero:

151 return None

152 if (coeff.has(x) or coeff.has(zoo) or coeff.has(oo) or

153 coeff.has(nan)):

154 return None

155

156 for j in range(i):

157 coeff = (coeff / (k + j + 1))

158 sep = integrate(sep, x)

159 sep += (ds.pop() - sep).limit(x, 0) # constant of integration

160 return (coeff.subs(k, k - i), sep, i)

161

162 else:

163 ds.append(diff)

164

165 return None

166

167

168def rational_independent(terms, x):

169 """

170 Returns a list of all the rationally independent terms.

171

172 Examples

173 ========

174

175 >>> from sympy import sin, cos

176 >>> from sympy.series.formal import rational_independent

177 >>> from sympy.abc import x

178

179 >>> rational_independent([cos(x), sin(x)], x)

180 [cos(x), sin(x)]

181 >>> rational_independent([x**2, sin(x), x*sin(x), x**3], x)

182 [x**3 + x**2, x*sin(x) + sin(x)]

183 """

184 if not terms:

185 return []

186

187 ind = terms[0:1]

188

189 for t in terms[1:]:

190 n = t.as_independent(x)[1]

191 for i, term in enumerate(ind):

192 d = term.as_independent(x)[1]

193 q = (n / d).cancel()

194 if q.is_rational_function(x):

195 ind[i] += t

196 break

197 else:

198 ind.append(t)

199 return ind

200

201

202def simpleDE(f, x, g, order=4):

203 r"""

204 Generates simple DE.

205

206 Explanation

207 ===========

208

209 DE is of the form

210

211 .. math::

212 f^k(x) + \sum\limits_{j=0}^{k-1} A_j f^j(x) = 0

213

214 where :math:`A_j` should be rational function in x.

215

216 Generates DE's upto order 4 (default). DE's can also have free parameters.

217

218 By increasing order, higher order DE's can be found.

219

220 Yields a tuple of (DE, order).

221 """

222 from sympy.solvers.solveset import linsolve

223

224 a = symbols('a:%d' % (order))

225

226 def _makeDE(k):

227 eq = f.diff(x, k) + Add(*[a[i]*f.diff(x, i) for i in range(0, k)])

228 DE = g(x).diff(x, k) + Add(*[a[i]*g(x).diff(x, i) for i in range(0, k)])

229 return eq, DE

230

231 found = False

232 for k in range(1, order + 1):

233 eq, DE = _makeDE(k)

234 eq = eq.expand()

235 terms = eq.as_ordered_terms()

236 ind = rational_independent(terms, x)

237 if found or len(ind) == k:

238 sol = dict(zip(a, (i for s in linsolve(ind, a[:k]) for i in s)))

239 if sol:

240 found = True

241 DE = DE.subs(sol)

242 DE = DE.as_numer_denom()[0]

243 DE = DE.factor().as_coeff_mul(Derivative)[1][0]

244 yield DE.collect(Derivative(g(x))), k

245

246

247def exp_re(DE, r, k):

248 """Converts a DE with constant coefficients (explike) into a RE.

249

250 Explanation

251 ===========

252

253 Performs the substitution:

254

255 .. math::

256 f^j(x) \\to r(k + j)

257

258 Normalises the terms so that lowest order of a term is always r(k).

259

260 Examples

261 ========

262

263 >>> from sympy import Function, Derivative

264 >>> from sympy.series.formal import exp_re

265 >>> from sympy.abc import x, k

266 >>> f, r = Function('f'), Function('r')

267

268 >>> exp_re(-f(x) + Derivative(f(x)), r, k)

269 -r(k) + r(k + 1)

270 >>> exp_re(Derivative(f(x), x) + Derivative(f(x), (x, 2)), r, k)

271 r(k) + r(k + 1)

272

273 See Also

274 ========

275

276 sympy.series.formal.hyper_re

277 """

278 RE = S.Zero

279

280 g = DE.atoms(Function).pop()

281

282 mini = None

283 for t in Add.make_args(DE):

284 coeff, d = t.as_independent(g)

285 if isinstance(d, Derivative):

286 j = d.derivative_count

287 else:

288 j = 0

289 if mini is None or j < mini:

290 mini = j

291 RE += coeff * r(k + j)

292 if mini:

293 RE = RE.subs(k, k - mini)

294 return RE

295

296

297def hyper_re(DE, r, k):

298 """

299 Converts a DE into a RE.

300

301 Explanation

302 ===========

303

304 Performs the substitution:

305

306 .. math::

307 x^l f^j(x) \\to (k + 1 - l)_j . a_{k + j - l}

308

309 Normalises the terms so that lowest order of a term is always r(k).

310

311 Examples

312 ========

313

314 >>> from sympy import Function, Derivative

315 >>> from sympy.series.formal import hyper_re

316 >>> from sympy.abc import x, k

317 >>> f, r = Function('f'), Function('r')

318

319 >>> hyper_re(-f(x) + Derivative(f(x)), r, k)

320 (k + 1)*r(k + 1) - r(k)

321 >>> hyper_re(-x*f(x) + Derivative(f(x), (x, 2)), r, k)

322 (k + 2)*(k + 3)*r(k + 3) - r(k)

323

324 See Also

325 ========

326

327 sympy.series.formal.exp_re

328 """

329 RE = S.Zero

330

331 g = DE.atoms(Function).pop()

332 x = g.atoms(Symbol).pop()

333

334 mini = None

335 for t in Add.make_args(DE.expand()):

336 coeff, d = t.as_independent(g)

337 c, v = coeff.as_independent(x)

338 l = v.as_coeff_exponent(x)[1]

339 if isinstance(d, Derivative):

340 j = d.derivative_count

341 else:

342 j = 0

343 RE += c * rf(k + 1 - l, j) * r(k + j - l)

344 if mini is None or j - l < mini:

345 mini = j - l

346

347 RE = RE.subs(k, k - mini)

348

349 m = Wild('m')

350 return RE.collect(r(k + m))

351

352

353def _transformation_a(f, x, P, Q, k, m, shift):

354 f *= x**(-shift)

355 P = P.subs(k, k + shift)

356 Q = Q.subs(k, k + shift)

357 return f, P, Q, m

358

359

360def _transformation_c(f, x, P, Q, k, m, scale):

361 f = f.subs(x, x**scale)

362 P = P.subs(k, k / scale)

363 Q = Q.subs(k, k / scale)

364 m *= scale

365 return f, P, Q, m

366

367

368def _transformation_e(f, x, P, Q, k, m):

369 f = f.diff(x)

370 P = P.subs(k, k + 1) * (k + m + 1)

371 Q = Q.subs(k, k + 1) * (k + 1)

372 return f, P, Q, m

373

374

375def _apply_shift(sol, shift):

376 return [(res, cond + shift) for res, cond in sol]

377

378

379def _apply_scale(sol, scale):

380 return [(res, cond / scale) for res, cond in sol]

381

382

383def _apply_integrate(sol, x, k):

384 return [(res / ((cond + 1)*(cond.as_coeff_Add()[1].coeff(k))), cond + 1)

385 for res, cond in sol]

386

387

388def _compute_formula(f, x, P, Q, k, m, k_max):

389 """Computes the formula for f."""

390 from sympy.polys import roots

391

392 sol = []

393 for i in range(k_max + 1, k_max + m + 1):

394 if (i < 0) == True:

395 continue

396 r = f.diff(x, i).limit(x, 0) / factorial(i)

397 if r.is_zero:

398 continue

399

400 kterm = m*k + i

401 res = r

402

403 p = P.subs(k, kterm)

404 q = Q.subs(k, kterm)

405 c1 = p.subs(k, 1/k).leadterm(k)[0]

406 c2 = q.subs(k, 1/k).leadterm(k)[0]

407 res *= (-c1 / c2)**k

408

409 res *= Mul(*[rf(-r, k)**mul for r, mul in roots(p, k).items()])

410 res /= Mul(*[rf(-r, k)**mul for r, mul in roots(q, k).items()])

411

412 sol.append((res, kterm))

413

414 return sol

415

416

417def _rsolve_hypergeometric(f, x, P, Q, k, m):

418 """

419 Recursive wrapper to rsolve_hypergeometric.

420

421 Explanation

422 ===========

423

424 Returns a Tuple of (formula, series independent terms,

425 maximum power of x in independent terms) if successful

426 otherwise ``None``.

427

428 See :func:`rsolve_hypergeometric` for details.

429 """

430 from sympy.polys import lcm, roots

431 from sympy.integrals import integrate

432

433 # transformation - c

434 proots, qroots = roots(P, k), roots(Q, k)

435 all_roots = dict(proots)

436 all_roots.update(qroots)

437 scale = lcm([r.as_numer_denom()[1] for r, t in all_roots.items()

438 if r.is_rational])

439 f, P, Q, m = _transformation_c(f, x, P, Q, k, m, scale)

440

441 # transformation - a

442 qroots = roots(Q, k)

443 if qroots:

444 k_min = Min(*qroots.keys())

445 else:

446 k_min = S.Zero

447 shift = k_min + m

448 f, P, Q, m = _transformation_a(f, x, P, Q, k, m, shift)

449

450 l = (x*f).limit(x, 0)

451 if not isinstance(l, Limit) and l != 0: # Ideally should only be l != 0

452 return None

453

454 qroots = roots(Q, k)

455 if qroots:

456 k_max = Max(*qroots.keys())

457 else:

458 k_max = S.Zero

459

460 ind, mp = S.Zero, -oo

461 for i in range(k_max + m + 1):

462 r = f.diff(x, i).limit(x, 0) / factorial(i)

463 if r.is_finite is False:

464 old_f = f

465 f, P, Q, m = _transformation_a(f, x, P, Q, k, m, i)

466 f, P, Q, m = _transformation_e(f, x, P, Q, k, m)

467 sol, ind, mp = _rsolve_hypergeometric(f, x, P, Q, k, m)

468 sol = _apply_integrate(sol, x, k)

469 sol = _apply_shift(sol, i)

470 ind = integrate(ind, x)

471 ind += (old_f - ind).limit(x, 0) # constant of integration

472 mp += 1

473 return sol, ind, mp

474 elif r:

475 ind += r*x**(i + shift)

476 pow_x = Rational((i + shift), scale)

477 if pow_x > mp:

478 mp = pow_x # maximum power of x

479 ind = ind.subs(x, x**(1/scale))

480

481 sol = _compute_formula(f, x, P, Q, k, m, k_max)

482 sol = _apply_shift(sol, shift)

483 sol = _apply_scale(sol, scale)

484

485 return sol, ind, mp

486

487

488def rsolve_hypergeometric(f, x, P, Q, k, m):

489 """

490 Solves RE of hypergeometric type.

491

492 Explanation

493 ===========

494

495 Attempts to solve RE of the form

496

497 Q(k)*a(k + m) - P(k)*a(k)

498

499 Transformations that preserve Hypergeometric type:

500

501 a. x**n*f(x): b(k + m) = R(k - n)*b(k)

502 b. f(A*x): b(k + m) = A**m*R(k)*b(k)

503 c. f(x**n): b(k + n*m) = R(k/n)*b(k)

504 d. f(x**(1/m)): b(k + 1) = R(k*m)*b(k)

505 e. f'(x): b(k + m) = ((k + m + 1)/(k + 1))*R(k + 1)*b(k)

506

507 Some of these transformations have been used to solve the RE.

508

509 Returns

510 =======

511

512 formula : Expr

513 ind : Expr

514 Independent terms.

515 order : int

516

517 Examples

518 ========

519

520 >>> from sympy import exp, ln, S

521 >>> from sympy.series.formal import rsolve_hypergeometric as rh

522 >>> from sympy.abc import x, k

523

524 >>> rh(exp(x), x, -S.One, (k + 1), k, 1)

525 (Piecewise((1/factorial(k), Eq(Mod(k, 1), 0)), (0, True)), 1, 1)

526

527 >>> rh(ln(1 + x), x, k**2, k*(k + 1), k, 1)

528 (Piecewise(((-1)**(k - 1)*factorial(k - 1)/RisingFactorial(2, k - 1),

529 Eq(Mod(k, 1), 0)), (0, True)), x, 2)

530

531 References

532 ==========

533

534 .. [1] Formal Power Series - Dominik Gruntz, Wolfram Koepf

535 .. [2] Power Series in Computer Algebra - Wolfram Koepf

536 """

537 result = _rsolve_hypergeometric(f, x, P, Q, k, m)

538

539 if result is None:

540 return None

541

542 sol_list, ind, mp = result

543

544 sol_dict = defaultdict(lambda: S.Zero)

545 for res, cond in sol_list:

546 j, mk = cond.as_coeff_Add()

547 c = mk.coeff(k)

548

549 if j.is_integer is False:

550 res *= x**frac(j)

551 j = floor(j)

552

553 res = res.subs(k, (k - j) / c)

554 cond = Eq(k % c, j % c)

555 sol_dict[cond] += res # Group together formula for same conditions

556

557 sol = []

558 for cond, res in sol_dict.items():

559 sol.append((res, cond))

560 sol.append((S.Zero, True))

561 sol = Piecewise(*sol)

562

563 if mp is -oo:

564 s = S.Zero

565 elif mp.is_integer is False:

566 s = ceiling(mp)

567 else:

568 s = mp + 1

569

570 # save all the terms of

571 # form 1/x**k in ind

572 if s < 0:

573 ind += sum(sequence(sol * x**k, (k, s, -1)))

574 s = S.Zero

575

576 return (sol, ind, s)

577

578

579def _solve_hyper_RE(f, x, RE, g, k):

580 """See docstring of :func:`rsolve_hypergeometric` for details."""

581 terms = Add.make_args(RE)

582

583 if len(terms) == 2:

584 gs = list(RE.atoms(Function))

585 P, Q = map(RE.coeff, gs)

586 m = gs[1].args[0] - gs[0].args[0]

587 if m < 0:

588 P, Q = Q, P

589 m = abs(m)

590 return rsolve_hypergeometric(f, x, P, Q, k, m)

591

592

593def _solve_explike_DE(f, x, DE, g, k):

594 """Solves DE with constant coefficients."""

595 from sympy.solvers import rsolve

596

597 for t in Add.make_args(DE):

598 coeff, d = t.as_independent(g)

599 if coeff.free_symbols:

600 return

601

602 RE = exp_re(DE, g, k)

603

604 init = {}

605 for i in range(len(Add.make_args(RE))):

606 if i:

607 f = f.diff(x)

608 init[g(k).subs(k, i)] = f.limit(x, 0)

609

610 sol = rsolve(RE, g(k), init)

611

612 if sol:

613 return (sol / factorial(k), S.Zero, S.Zero)

614

615

616def _solve_simple(f, x, DE, g, k):

617 """Converts DE into RE and solves using :func:`rsolve`."""

618 from sympy.solvers import rsolve

619

620 RE = hyper_re(DE, g, k)

621

622 init = {}

623 for i in range(len(Add.make_args(RE))):

624 if i:

625 f = f.diff(x)

626 init[g(k).subs(k, i)] = f.limit(x, 0) / factorial(i)

627

628 sol = rsolve(RE, g(k), init)

629

630 if sol:

631 return (sol, S.Zero, S.Zero)

632

633

634def _transform_explike_DE(DE, g, x, order, syms):

635 """Converts DE with free parameters into DE with constant coefficients."""

636 from sympy.solvers.solveset import linsolve

637

638 eq = []

639 highest_coeff = DE.coeff(Derivative(g(x), x, order))

640 for i in range(order):

641 coeff = DE.coeff(Derivative(g(x), x, i))

642 coeff = (coeff / highest_coeff).expand().collect(x)

643 for t in Add.make_args(coeff):

644 eq.append(t)

645 temp = []

646 for e in eq:

647 if e.has(x):

648 break

649 elif e.has(Symbol):

650 temp.append(e)

651 else:

652 eq = temp

653 if eq:

654 sol = dict(zip(syms, (i for s in linsolve(eq, list(syms)) for i in s)))

655 if sol:

656 DE = DE.subs(sol)

657 DE = DE.factor().as_coeff_mul(Derivative)[1][0]

658 DE = DE.collect(Derivative(g(x)))

659 return DE

660

661

662def _transform_DE_RE(DE, g, k, order, syms):

663 """Converts DE with free parameters into RE of hypergeometric type."""

664 from sympy.solvers.solveset import linsolve

665

666 RE = hyper_re(DE, g, k)

667

668 eq = []

669 for i in range(1, order):

670 coeff = RE.coeff(g(k + i))

671 eq.append(coeff)

672 sol = dict(zip(syms, (i for s in linsolve(eq, list(syms)) for i in s)))

673 if sol:

674 m = Wild('m')

675 RE = RE.subs(sol)

676 RE = RE.factor().as_numer_denom()[0].collect(g(k + m))

677 RE = RE.as_coeff_mul(g)[1][0]

678 for i in range(order): # smallest order should be g(k)

679 if RE.coeff(g(k + i)) and i:

680 RE = RE.subs(k, k - i)

681 break

682 return RE

683

684

685def solve_de(f, x, DE, order, g, k):

686 """

687 Solves the DE.

688

689 Explanation

690 ===========

691

692 Tries to solve DE by either converting into a RE containing two terms or

693 converting into a DE having constant coefficients.

694

695 Returns

696 =======

697

698 formula : Expr

699 ind : Expr

700 Independent terms.

701 order : int

702

703 Examples

704 ========

705

706 >>> from sympy import Derivative as D, Function

707 >>> from sympy import exp, ln

708 >>> from sympy.series.formal import solve_de

709 >>> from sympy.abc import x, k

710 >>> f = Function('f')

711

712 >>> solve_de(exp(x), x, D(f(x), x) - f(x), 1, f, k)

713 (Piecewise((1/factorial(k), Eq(Mod(k, 1), 0)), (0, True)), 1, 1)

714

715 >>> solve_de(ln(1 + x), x, (x + 1)*D(f(x), x, 2) + D(f(x)), 2, f, k)

716 (Piecewise(((-1)**(k - 1)*factorial(k - 1)/RisingFactorial(2, k - 1),

717 Eq(Mod(k, 1), 0)), (0, True)), x, 2)

718 """

719 sol = None

720 syms = DE.free_symbols.difference({g, x})

721

722 if syms:

723 RE = _transform_DE_RE(DE, g, k, order, syms)

724 else:

725 RE = hyper_re(DE, g, k)

726 if not RE.free_symbols.difference({k}):

727 sol = _solve_hyper_RE(f, x, RE, g, k)

728

729 if sol:

730 return sol

731

732 if syms:

733 DE = _transform_explike_DE(DE, g, x, order, syms)

734 if not DE.free_symbols.difference({x}):

735 sol = _solve_explike_DE(f, x, DE, g, k)

736

737 if sol:

738 return sol

739

740

741def hyper_algorithm(f, x, k, order=4):

742 """

743 Hypergeometric algorithm for computing Formal Power Series.

744

745 Explanation

746 ===========

747

748 Steps:

749 * Generates DE

750 * Convert the DE into RE

751 * Solves the RE

752

753 Examples

754 ========

755

756 >>> from sympy import exp, ln

757 >>> from sympy.series.formal import hyper_algorithm

758

759 >>> from sympy.abc import x, k

760

761 >>> hyper_algorithm(exp(x), x, k)

762 (Piecewise((1/factorial(k), Eq(Mod(k, 1), 0)), (0, True)), 1, 1)

763

764 >>> hyper_algorithm(ln(1 + x), x, k)

765 (Piecewise(((-1)**(k - 1)*factorial(k - 1)/RisingFactorial(2, k - 1),

766 Eq(Mod(k, 1), 0)), (0, True)), x, 2)

767

768 See Also

769 ========

770

771 sympy.series.formal.simpleDE

772 sympy.series.formal.solve_de

773 """

774 g = Function('g')

775

776 des = [] # list of DE's

777 sol = None

778 for DE, i in simpleDE(f, x, g, order):

779 if DE is not None:

780 sol = solve_de(f, x, DE, i, g, k)

781 if sol:

782 return sol

783 if not DE.free_symbols.difference({x}):

784 des.append(DE)

785

786 # If nothing works

787 # Try plain rsolve

788 for DE in des:

789 sol = _solve_simple(f, x, DE, g, k)

790 if sol:

791 return sol

792

793

794def _compute_fps(f, x, x0, dir, hyper, order, rational, full):

795 """Recursive wrapper to compute fps.

796

797 See :func:`compute_fps` for details.

798 """

799 if x0 in [S.Infinity, S.NegativeInfinity]:

800 dir = S.One if x0 is S.Infinity else -S.One

801 temp = f.subs(x, 1/x)

802 result = _compute_fps(temp, x, 0, dir, hyper, order, rational, full)

803 if result is None:

804 return None

805 return (result[0], result[1].subs(x, 1/x), result[2].subs(x, 1/x))

806 elif x0 or dir == -S.One:

807 if dir == -S.One:

808 rep = -x + x0

809 rep2 = -x

810 rep2b = x0

811 else:

812 rep = x + x0

813 rep2 = x

814 rep2b = -x0

815 temp = f.subs(x, rep)

816 result = _compute_fps(temp, x, 0, S.One, hyper, order, rational, full)

817 if result is None:

818 return None

819 return (result[0], result[1].subs(x, rep2 + rep2b),

820 result[2].subs(x, rep2 + rep2b))

821

822 if f.is_polynomial(x):

823 k = Dummy('k')

824 ak = sequence(Coeff(f, x, k), (k, 1, oo))

825 xk = sequence(x**k, (k, 0, oo))

826 ind = f.coeff(x, 0)

827 return ak, xk, ind

828

829 # Break instances of Add

830 # this allows application of different

831 # algorithms on different terms increasing the

832 # range of admissible functions.

833 if isinstance(f, Add):

834 result = False

835 ak = sequence(S.Zero, (0, oo))

836 ind, xk = S.Zero, None

837 for t in Add.make_args(f):

838 res = _compute_fps(t, x, 0, S.One, hyper, order, rational, full)

839 if res:

840 if not result:

841 result = True

842 xk = res[1]

843 if res[0].start > ak.start:

844 seq = ak

845 s, f = ak.start, res[0].start

846 else:

847 seq = res[0]

848 s, f = res[0].start, ak.start

849 save = Add(*[z[0]*z[1] for z in zip(seq[0:(f - s)], xk[s:f])])

850 ak += res[0]

851 ind += res[2] + save

852 else:

853 ind += t

854 if result:

855 return ak, xk, ind

856 return None

857

858 # The symbolic term - symb, if present, is being separated from the function

859 # Otherwise symb is being set to S.One

860 syms = f.free_symbols.difference({x})

861 (f, symb) = expand(f).as_independent(*syms)

862

863 result = None

864

865 # from here on it's x0=0 and dir=1 handling

866 k = Dummy('k')

867 if rational:

868 result = rational_algorithm(f, x, k, order, full)

869

870 if result is None and hyper:

871 result = hyper_algorithm(f, x, k, order)

872

873 if result is None:

874 return None

875

876 from sympy.simplify.powsimp import powsimp

877 if symb.is_zero:

878 symb = S.One

879 else:

880 symb = powsimp(symb)

881 ak = sequence(result[0], (k, result[2], oo))

882 xk_formula = powsimp(x**k * symb)

883 xk = sequence(xk_formula, (k, 0, oo))

884 ind = powsimp(result[1] * symb)

885

886 return ak, xk, ind

887

888

889def compute_fps(f, x, x0=0, dir=1, hyper=True, order=4, rational=True,

890 full=False):

891 """

892 Computes the formula for Formal Power Series of a function.

893

894 Explanation

895 ===========

896

897 Tries to compute the formula by applying the following techniques

898 (in order):

899

900 * rational_algorithm

901 * Hypergeometric algorithm

902

903 Parameters

904 ==========

905

906 x : Symbol

907 x0 : number, optional

908 Point to perform series expansion about. Default is 0.

909 dir : {1, -1, '+', '-'}, optional

910 If dir is 1 or '+' the series is calculated from the right and

911 for -1 or '-' the series is calculated from the left. For smooth

912 functions this flag will not alter the results. Default is 1.

913 hyper : {True, False}, optional

914 Set hyper to False to skip the hypergeometric algorithm.

915 By default it is set to False.

916 order : int, optional

917 Order of the derivative of ``f``, Default is 4.

918 rational : {True, False}, optional

919 Set rational to False to skip rational algorithm. By default it is set

920 to True.

921 full : {True, False}, optional

922 Set full to True to increase the range of rational algorithm.

923 See :func:`rational_algorithm` for details. By default it is set to

924 False.

925

926 Returns

927 =======

928

929 ak : sequence

930 Sequence of coefficients.

931 xk : sequence

932 Sequence of powers of x.

933 ind : Expr

934 Independent terms.

935 mul : Pow

936 Common terms.

937

938 See Also

939 ========

940

941 sympy.series.formal.rational_algorithm

942 sympy.series.formal.hyper_algorithm

943 """

944 f = sympify(f)

945 x = sympify(x)

946

947 if not f.has(x):

948 return None

949

950 x0 = sympify(x0)

951

952 if dir == '+':

953 dir = S.One

954 elif dir == '-':

955 dir = -S.One

956 elif dir not in [S.One, -S.One]:

957 raise ValueError("Dir must be '+' or '-'")

958 else:

959 dir = sympify(dir)

960

961 return _compute_fps(f, x, x0, dir, hyper, order, rational, full)

962

963

964class Coeff(Function):

965 """

966 Coeff(p, x, n) represents the nth coefficient of the polynomial p in x

967 """

968 @classmethod

969 def eval(cls, p, x, n):

970 if p.is_polynomial(x) and n.is_integer:

971 return p.coeff(x, n)

972

973

974class FormalPowerSeries(SeriesBase):

975 """

976 Represents Formal Power Series of a function.

977

978 Explanation

979 ===========

980

981 No computation is performed. This class should only to be used to represent

982 a series. No checks are performed.

983

984 For computing a series use :func:`fps`.

985

986 See Also

987 ========

988

989 sympy.series.formal.fps

990 """

991 def __new__(cls, *args):

992 args = map(sympify, args)

993 return Expr.__new__(cls, *args)

994

995 def __init__(self, *args):

996 ak = args[4][0]

997 k = ak.variables[0]

998 self.ak_seq = sequence(ak.formula, (k, 1, oo))

999 self.fact_seq = sequence(factorial(k), (k, 1, oo))

1000 self.bell_coeff_seq = self.ak_seq * self.fact_seq

1001 self.sign_seq = sequence((-1, 1), (k, 1, oo))

1002

1003 @property

1004 def function(self):

1005 return self.args[0]

1006

1007 @property

1008 def x(self):

1009 return self.args[1]

1010

1011 @property

1012 def x0(self):

1013 return self.args[2]

1014

1015 @property

1016 def dir(self):

1017 return self.args[3]

1018

1019 @property

1020 def ak(self):

1021 return self.args[4][0]

1022

1023 @property

1024 def xk(self):

1025 return self.args[4][1]

1026

1027 @property

1028 def ind(self):

1029 return self.args[4][2]

1030

1031 @property

1032 def interval(self):

1033 return Interval(0, oo)

1034

1035 @property

1036 def start(self):

1037 return self.interval.inf

1038

1039 @property

1040 def stop(self):

1041 return self.interval.sup

1042

1043 @property

1044 def length(self):

1045 return oo

1046

1047 @property

1048 def infinite(self):

1049 """Returns an infinite representation of the series"""

1050 from sympy.concrete import Sum

1051 ak, xk = self.ak, self.xk

1052 k = ak.variables[0]

1053 inf_sum = Sum(ak.formula * xk.formula, (k, ak.start, ak.stop))

1054

1055 return self.ind + inf_sum

1056

1057 def _get_pow_x(self, term):

1058 """Returns the power of x in a term."""

1059 xterm, pow_x = term.as_independent(self.x)[1].as_base_exp()

1060 if not xterm.has(self.x):

1061 return S.Zero

1062 return pow_x

1063

1064 def polynomial(self, n=6):

1065 """

1066 Truncated series as polynomial.

1067

1068 Explanation

1069 ===========

1070

1071 Returns series expansion of ``f`` upto order ``O(x**n)``

1072 as a polynomial(without ``O`` term).

1073 """

1074 terms = []

1075 sym = self.free_symbols

1076 for i, t in enumerate(self):

1077 xp = self._get_pow_x(t)

1078 if xp.has(*sym):

1079 xp = xp.as_coeff_add(*sym)[0]

1080 if xp >= n:

1081 break

1082 elif xp.is_integer is True and i == n + 1:

1083 break

1084 elif t is not S.Zero:

1085 terms.append(t)

1086

1087 return Add(*terms)

1088

1089 def truncate(self, n=6):

1090 """

1091 Truncated series.

1092

1093 Explanation

1094 ===========

1095

1096 Returns truncated series expansion of f upto

1097 order ``O(x**n)``.

1098

1099 If n is ``None``, returns an infinite iterator.

1100 """

1101 if n is None:

1102 return iter(self)

1103

1104 x, x0 = self.x, self.x0

1105 pt_xk = self.xk.coeff(n)

1106 if x0 is S.NegativeInfinity:

1107 x0 = S.Infinity

1108

1109 return self.polynomial(n) + Order(pt_xk, (x, x0))

1110

1111 def zero_coeff(self):

1112 return self._eval_term(0)

1113

1114 def _eval_term(self, pt):

1115 try:

1116 pt_xk = self.xk.coeff(pt)

1117 pt_ak = self.ak.coeff(pt).simplify() # Simplify the coefficients

1118 except IndexError:

1119 term = S.Zero

1120 else:

1121 term = (pt_ak * pt_xk)

1122

1123 if self.ind:

1124 ind = S.Zero

1125 sym = self.free_symbols

1126 for t in Add.make_args(self.ind):

1127 pow_x = self._get_pow_x(t)

1128 if pow_x.has(*sym):

1129 pow_x = pow_x.as_coeff_add(*sym)[0]

1130 if pt == 0 and pow_x < 1:

1131 ind += t

1132 elif pow_x >= pt and pow_x < pt + 1:

1133 ind += t

1134 term += ind

1135

1136 return term.collect(self.x)

1137

1138 def _eval_subs(self, old, new):

1139 x = self.x

1140 if old.has(x):

1141 return self

1142

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

1144 for t in self:

1145 if t is not S.Zero:

1146 return t

1147

1148 def _eval_derivative(self, x):

1149 f = self.function.diff(x)

1150 ind = self.ind.diff(x)

1151

1152 pow_xk = self._get_pow_x(self.xk.formula)

1153 ak = self.ak

1154 k = ak.variables[0]

1155 if ak.formula.has(x):

1156 form = []

1157 for e, c in ak.formula.args:

1158 temp = S.Zero

1159 for t in Add.make_args(e):

1160 pow_x = self._get_pow_x(t)

1161 temp += t * (pow_xk + pow_x)

1162 form.append((temp, c))

1163 form = Piecewise(*form)

1164 ak = sequence(form.subs(k, k + 1), (k, ak.start - 1, ak.stop))

1165 else:

1166 ak = sequence((ak.formula * pow_xk).subs(k, k + 1),

1167 (k, ak.start - 1, ak.stop))

1168

1169 return self.func(f, self.x, self.x0, self.dir, (ak, self.xk, ind))

1170

1171 def integrate(self, x=None, **kwargs):

1172 """

1173 Integrate Formal Power Series.

1174

1175 Examples

1176 ========

1177

1178 >>> from sympy import fps, sin, integrate

1179 >>> from sympy.abc import x

1180 >>> f = fps(sin(x))

1181 >>> f.integrate(x).truncate()

1182 -1 + x**2/2 - x**4/24 + O(x**6)

1183 >>> integrate(f, (x, 0, 1))

1184 1 - cos(1)

1185 """

1186 from sympy.integrals import integrate

1187

1188 if x is None:

1189 x = self.x

1190 elif iterable(x):

1191 return integrate(self.function, x)

1192

1193 f = integrate(self.function, x)

1194 ind = integrate(self.ind, x)

1195 ind += (f - ind).limit(x, 0) # constant of integration

1196

1197 pow_xk = self._get_pow_x(self.xk.formula)

1198 ak = self.ak

1199 k = ak.variables[0]

1200 if ak.formula.has(x):

1201 form = []

1202 for e, c in ak.formula.args:

1203 temp = S.Zero

1204 for t in Add.make_args(e):

1205 pow_x = self._get_pow_x(t)

1206 temp += t / (pow_xk + pow_x + 1)

1207 form.append((temp, c))

1208 form = Piecewise(*form)

1209 ak = sequence(form.subs(k, k - 1), (k, ak.start + 1, ak.stop))

1210 else:

1211 ak = sequence((ak.formula / (pow_xk + 1)).subs(k, k - 1),

1212 (k, ak.start + 1, ak.stop))

1213

1214 return self.func(f, self.x, self.x0, self.dir, (ak, self.xk, ind))

1215

1216 def product(self, other, x=None, n=6):

1217 """

1218 Multiplies two Formal Power Series, using discrete convolution and

1219 return the truncated terms upto specified order.

1220

1221 Parameters

1222 ==========

1223

1224 n : Number, optional

1225 Specifies the order of the term up to which the polynomial should

1226 be truncated.

1227

1228 Examples

1229 ========

1230

1231 >>> from sympy import fps, sin, exp

1232 >>> from sympy.abc import x

1233 >>> f1 = fps(sin(x))

1234 >>> f2 = fps(exp(x))

1235

1236 >>> f1.product(f2, x).truncate(4)

1237 x + x**2 + x**3/3 + O(x**4)

1238

1239 See Also

1240 ========

1241

1242 sympy.discrete.convolutions

1243 sympy.series.formal.FormalPowerSeriesProduct

1244

1245 """

1246

1247 if n is None:

1248 return iter(self)

1249

1250 other = sympify(other)

1251

1252 if not isinstance(other, FormalPowerSeries):

1253 raise ValueError("Both series should be an instance of FormalPowerSeries"

1254 " class.")

1255

1256 if self.dir != other.dir:

1257 raise ValueError("Both series should be calculated from the"

1258 " same direction.")

1259 elif self.x0 != other.x0:

1260 raise ValueError("Both series should be calculated about the"

1261 " same point.")

1262

1263 elif self.x != other.x:

1264 raise ValueError("Both series should have the same symbol.")

1265

1266 return FormalPowerSeriesProduct(self, other)

1267

1268 def coeff_bell(self, n):

1269 r"""

1270 self.coeff_bell(n) returns a sequence of Bell polynomials of the second kind.

1271 Note that ``n`` should be a integer.

1272

1273 The second kind of Bell polynomials (are sometimes called "partial" Bell

1274 polynomials or incomplete Bell polynomials) are defined as

1275

1276 .. math::

1277 B_{n,k}(x_1, x_2,\dotsc x_{n-k+1}) =

1278 \sum_{j_1+j_2+j_2+\dotsb=k \atop j_1+2j_2+3j_2+\dotsb=n}

1279 \frac{n!}{j_1!j_2!\dotsb j_{n-k+1}!}

1280 \left(\frac{x_1}{1!} \right)^{j_1}

1281 \left(\frac{x_2}{2!} \right)^{j_2} \dotsb

1282 \left(\frac{x_{n-k+1}}{(n-k+1)!} \right) ^{j_{n-k+1}}.

1283

1284 * ``bell(n, k, (x1, x2, ...))`` gives Bell polynomials of the second kind,

1285 `B_{n,k}(x_1, x_2, \dotsc, x_{n-k+1})`.

1286

1287 See Also

1288 ========

1289

1290 sympy.functions.combinatorial.numbers.bell

1291

1292 """

1293

1294 inner_coeffs = [bell(n, j, tuple(self.bell_coeff_seq[:n-j+1])) for j in range(1, n+1)]

1295

1296 k = Dummy('k')

1297 return sequence(tuple(inner_coeffs), (k, 1, oo))

1298

1299 def compose(self, other, x=None, n=6):

1300 r"""

1301 Returns the truncated terms of the formal power series of the composed function,

1302 up to specified ``n``.

1303

1304 Explanation

1305 ===========

1306

1307 If ``f`` and ``g`` are two formal power series of two different functions,

1308 then the coefficient sequence ``ak`` of the composed formal power series `fp`

1309 will be as follows.

1310

1311 .. math::

1312 \sum\limits_{k=0}^{n} b_k B_{n,k}(x_1, x_2, \dotsc, x_{n-k+1})

1313

1314 Parameters

1315 ==========

1316

1317 n : Number, optional

1318 Specifies the order of the term up to which the polynomial should

1319 be truncated.

1320

1321 Examples

1322 ========

1323

1324 >>> from sympy import fps, sin, exp

1325 >>> from sympy.abc import x

1326 >>> f1 = fps(exp(x))

1327 >>> f2 = fps(sin(x))

1328

1329 >>> f1.compose(f2, x).truncate()

1330 1 + x + x**2/2 - x**4/8 - x**5/15 + O(x**6)

1331

1332 >>> f1.compose(f2, x).truncate(8)

1333 1 + x + x**2/2 - x**4/8 - x**5/15 - x**6/240 + x**7/90 + O(x**8)

1334

1335 See Also

1336 ========

1337

1338 sympy.functions.combinatorial.numbers.bell

1339 sympy.series.formal.FormalPowerSeriesCompose

1340

1341 References

1342 ==========

1343

1344 .. [1] Comtet, Louis: Advanced combinatorics; the art of finite and infinite expansions. Reidel, 1974.

1345

1346 """

1347

1348 if n is None:

1349 return iter(self)

1350

1351 other = sympify(other)

1352

1353 if not isinstance(other, FormalPowerSeries):

1354 raise ValueError("Both series should be an instance of FormalPowerSeries"

1355 " class.")

1356

1357 if self.dir != other.dir:

1358 raise ValueError("Both series should be calculated from the"

1359 " same direction.")

1360 elif self.x0 != other.x0:

1361 raise ValueError("Both series should be calculated about the"

1362 " same point.")

1363

1364 elif self.x != other.x:

1365 raise ValueError("Both series should have the same symbol.")

1366

1367 if other._eval_term(0).as_coeff_mul(other.x)[0] is not S.Zero:

1368 raise ValueError("The formal power series of the inner function should not have any "

1369 "constant coefficient term.")

1370

1371 return FormalPowerSeriesCompose(self, other)

1372

1373 def inverse(self, x=None, n=6):

1374 r"""

1375 Returns the truncated terms of the inverse of the formal power series,

1376 up to specified ``n``.

1377

1378 Explanation

1379 ===========

1380

1381 If ``f`` and ``g`` are two formal power series of two different functions,

1382 then the coefficient sequence ``ak`` of the composed formal power series ``fp``

1383 will be as follows.

1384

1385 .. math::

1386 \sum\limits_{k=0}^{n} (-1)^{k} x_0^{-k-1} B_{n,k}(x_1, x_2, \dotsc, x_{n-k+1})

1387

1388 Parameters

1389 ==========

1390

1391 n : Number, optional

1392 Specifies the order of the term up to which the polynomial should

1393 be truncated.

1394

1395 Examples

1396 ========

1397

1398 >>> from sympy import fps, exp, cos

1399 >>> from sympy.abc import x

1400 >>> f1 = fps(exp(x))

1401 >>> f2 = fps(cos(x))

1402

1403 >>> f1.inverse(x).truncate()

1404 1 - x + x**2/2 - x**3/6 + x**4/24 - x**5/120 + O(x**6)

1405

1406 >>> f2.inverse(x).truncate(8)

1407 1 + x**2/2 + 5*x**4/24 + 61*x**6/720 + O(x**8)

1408

1409 See Also

1410 ========

1411

1412 sympy.functions.combinatorial.numbers.bell

1413 sympy.series.formal.FormalPowerSeriesInverse

1414

1415 References

1416 ==========

1417

1418 .. [1] Comtet, Louis: Advanced combinatorics; the art of finite and infinite expansions. Reidel, 1974.

1419

1420 """

1421

1422 if n is None:

1423 return iter(self)

1424

1425 if self._eval_term(0).is_zero:

1426 raise ValueError("Constant coefficient should exist for an inverse of a formal"

1427 " power series to exist.")

1428

1429 return FormalPowerSeriesInverse(self)

1430

1431 def __add__(self, other):

1432 other = sympify(other)

1433

1434 if isinstance(other, FormalPowerSeries):

1435 if self.dir != other.dir:

1436 raise ValueError("Both series should be calculated from the"

1437 " same direction.")

1438 elif self.x0 != other.x0:

1439 raise ValueError("Both series should be calculated about the"

1440 " same point.")

1441

1442 x, y = self.x, other.x

1443 f = self.function + other.function.subs(y, x)

1444

1445 if self.x not in f.free_symbols:

1446 return f

1447

1448 ak = self.ak + other.ak

1449 if self.ak.start > other.ak.start:

1450 seq = other.ak

1451 s, e = other.ak.start, self.ak.start

1452 else:

1453 seq = self.ak

1454 s, e = self.ak.start, other.ak.start

1455 save = Add(*[z[0]*z[1] for z in zip(seq[0:(e - s)], self.xk[s:e])])

1456 ind = self.ind + other.ind + save

1457

1458 return self.func(f, x, self.x0, self.dir, (ak, self.xk, ind))

1459

1460 elif not other.has(self.x):

1461 f = self.function + other

1462 ind = self.ind + other

1463

1464 return self.func(f, self.x, self.x0, self.dir,

1465 (self.ak, self.xk, ind))

1466

1467 return Add(self, other)

1468

1469 def __radd__(self, other):

1470 return self.__add__(other)

1471

1472 def __neg__(self):

1473 return self.func(-self.function, self.x, self.x0, self.dir,

1474 (-self.ak, self.xk, -self.ind))

1475

1476 def __sub__(self, other):

1477 return self.__add__(-other)

1478

1479 def __rsub__(self, other):

1480 return (-self).__add__(other)

1481

1482 def __mul__(self, other):

1483 other = sympify(other)

1484

1485 if other.has(self.x):

1486 return Mul(self, other)

1487

1488 f = self.function * other

1489 ak = self.ak.coeff_mul(other)

1490 ind = self.ind * other

1491

1492 return self.func(f, self.x, self.x0, self.dir, (ak, self.xk, ind))

1493

1494 def __rmul__(self, other):

1495 return self.__mul__(other)

1496

1497

1498class FiniteFormalPowerSeries(FormalPowerSeries):

1499 """Base Class for Product, Compose and Inverse classes"""

1500

1501 def __init__(self, *args):

1502 pass

1503

1504 @property

1505 def ffps(self):

1506 return self.args[0]

1507

1508 @property

1509 def gfps(self):

1510 return self.args[1]

1511

1512 @property

1513 def f(self):

1514 return self.ffps.function

1515

1516 @property

1517 def g(self):

1518 return self.gfps.function

1519

1520 @property

1521 def infinite(self):

1522 raise NotImplementedError("No infinite version for an object of"

1523 " FiniteFormalPowerSeries class.")

1524

1525 def _eval_terms(self, n):

1526 raise NotImplementedError("(%s)._eval_terms()" % self)

1527

1528 def _eval_term(self, pt):

1529 raise NotImplementedError("By the current logic, one can get terms"

1530 "upto a certain order, instead of getting term by term.")

1531

1532 def polynomial(self, n):

1533 return self._eval_terms(n)

1534

1535 def truncate(self, n=6):

1536 ffps = self.ffps

1537 pt_xk = ffps.xk.coeff(n)

1538 x, x0 = ffps.x, ffps.x0

1539

1540 return self.polynomial(n) + Order(pt_xk, (x, x0))

1541

1542 def _eval_derivative(self, x):

1543 raise NotImplementedError

1544

1545 def integrate(self, x):

1546 raise NotImplementedError

1547

1548

1549class FormalPowerSeriesProduct(FiniteFormalPowerSeries):

1550 """Represents the product of two formal power series of two functions.

1551

1552 Explanation

1553 ===========

1554

1555 No computation is performed. Terms are calculated using a term by term logic,

1556 instead of a point by point logic.

1557

1558 There are two differences between a :obj:`FormalPowerSeries` object and a

1559 :obj:`FormalPowerSeriesProduct` object. The first argument contains the two

1560 functions involved in the product. Also, the coefficient sequence contains

1561 both the coefficient sequence of the formal power series of the involved functions.

1562

1563 See Also

1564 ========

1565

1566 sympy.series.formal.FormalPowerSeries

1567 sympy.series.formal.FiniteFormalPowerSeries

1568

1569 """

1570

1571 def __init__(self, *args):

1572 ffps, gfps = self.ffps, self.gfps

1573

1574 k = ffps.ak.variables[0]

1575 self.coeff1 = sequence(ffps.ak.formula, (k, 0, oo))

1576

1577 k = gfps.ak.variables[0]

1578 self.coeff2 = sequence(gfps.ak.formula, (k, 0, oo))

1579

1580 @property

1581 def function(self):

1582 """Function of the product of two formal power series."""

1583 return self.f * self.g

1584

1585 def _eval_terms(self, n):

1586 """

1587 Returns the first ``n`` terms of the product formal power series.

1588 Term by term logic is implemented here.

1589

1590 Examples

1591 ========

1592

1593 >>> from sympy import fps, sin, exp

1594 >>> from sympy.abc import x

1595 >>> f1 = fps(sin(x))

1596 >>> f2 = fps(exp(x))

1597 >>> fprod = f1.product(f2, x)

1598

1599 >>> fprod._eval_terms(4)

1600 x**3/3 + x**2 + x

1601

1602 See Also

1603 ========

1604

1605 sympy.series.formal.FormalPowerSeries.product

1606

1607 """

1608 coeff1, coeff2 = self.coeff1, self.coeff2

1609

1610 aks = convolution(coeff1[:n], coeff2[:n])

1611

1612 terms = []

1613 for i in range(0, n):

1614 terms.append(aks[i] * self.ffps.xk.coeff(i))

1615

1616 return Add(*terms)

1617

1618

1619class FormalPowerSeriesCompose(FiniteFormalPowerSeries):

1620 """

1621 Represents the composed formal power series of two functions.

1622

1623 Explanation

1624 ===========

1625

1626 No computation is performed. Terms are calculated using a term by term logic,

1627 instead of a point by point logic.

1628

1629 There are two differences between a :obj:`FormalPowerSeries` object and a

1630 :obj:`FormalPowerSeriesCompose` object. The first argument contains the outer

1631 function and the inner function involved in the omposition. Also, the

1632 coefficient sequence contains the generic sequence which is to be multiplied

1633 by a custom ``bell_seq`` finite sequence. The finite terms will then be added up to

1634 get the final terms.

1635

1636 See Also

1637 ========

1638

1639 sympy.series.formal.FormalPowerSeries

1640 sympy.series.formal.FiniteFormalPowerSeries

1641

1642 """

1643

1644 @property

1645 def function(self):

1646 """Function for the composed formal power series."""

1647 f, g, x = self.f, self.g, self.ffps.x

1648 return f.subs(x, g)

1649

1650 def _eval_terms(self, n):

1651 """

1652 Returns the first `n` terms of the composed formal power series.

1653 Term by term logic is implemented here.

1654

1655 Explanation

1656 ===========

1657

1658 The coefficient sequence of the :obj:`FormalPowerSeriesCompose` object is the generic sequence.

1659 It is multiplied by ``bell_seq`` to get a sequence, whose terms are added up to get

1660 the final terms for the polynomial.

1661

1662 Examples

1663 ========

1664

1665 >>> from sympy import fps, sin, exp

1666 >>> from sympy.abc import x

1667 >>> f1 = fps(exp(x))

1668 >>> f2 = fps(sin(x))

1669 >>> fcomp = f1.compose(f2, x)

1670

1671 >>> fcomp._eval_terms(6)

1672 -x**5/15 - x**4/8 + x**2/2 + x + 1

1673

1674 >>> fcomp._eval_terms(8)

1675 x**7/90 - x**6/240 - x**5/15 - x**4/8 + x**2/2 + x + 1

1676

1677 See Also

1678 ========

1679

1680 sympy.series.formal.FormalPowerSeries.compose

1681 sympy.series.formal.FormalPowerSeries.coeff_bell

1682

1683 """

1684

1685 ffps, gfps = self.ffps, self.gfps

1686 terms = [ffps.zero_coeff()]

1687

1688 for i in range(1, n):

1689 bell_seq = gfps.coeff_bell(i)

1690 seq = (ffps.bell_coeff_seq * bell_seq)

1691 terms.append(Add(*(seq[:i])) / ffps.fact_seq[i-1] * ffps.xk.coeff(i))

1692

1693 return Add(*terms)

1694

1695

1696class FormalPowerSeriesInverse(FiniteFormalPowerSeries):

1697 """

1698 Represents the Inverse of a formal power series.

1699

1700 Explanation

1701 ===========

1702

1703 No computation is performed. Terms are calculated using a term by term logic,

1704 instead of a point by point logic.

1705

1706 There is a single difference between a :obj:`FormalPowerSeries` object and a

1707 :obj:`FormalPowerSeriesInverse` object. The coefficient sequence contains the

1708 generic sequence which is to be multiplied by a custom ``bell_seq`` finite sequence.

1709 The finite terms will then be added up to get the final terms.

1710

1711 See Also

1712 ========

1713

1714 sympy.series.formal.FormalPowerSeries

1715 sympy.series.formal.FiniteFormalPowerSeries

1716

1717 """

1718 def __init__(self, *args):

1719 ffps = self.ffps

1720 k = ffps.xk.variables[0]

1721

1722 inv = ffps.zero_coeff()

1723 inv_seq = sequence(inv ** (-(k + 1)), (k, 1, oo))

1724 self.aux_seq = ffps.sign_seq * ffps.fact_seq * inv_seq

1725

1726 @property

1727 def function(self):

1728 """Function for the inverse of a formal power series."""

1729 f = self.f

1730 return 1 / f

1731

1732 @property

1733 def g(self):

1734 raise ValueError("Only one function is considered while performing"

1735 "inverse of a formal power series.")

1736

1737 @property

1738 def gfps(self):

1739 raise ValueError("Only one function is considered while performing"

1740 "inverse of a formal power series.")

1741

1742 def _eval_terms(self, n):

1743 """

1744 Returns the first ``n`` terms of the composed formal power series.

1745 Term by term logic is implemented here.

1746

1747 Explanation

1748 ===========

1749

1750 The coefficient sequence of the `FormalPowerSeriesInverse` object is the generic sequence.

1751 It is multiplied by ``bell_seq`` to get a sequence, whose terms are added up to get

1752 the final terms for the polynomial.

1753

1754 Examples

1755 ========

1756

1757 >>> from sympy import fps, exp, cos

1758 >>> from sympy.abc import x

1759 >>> f1 = fps(exp(x))

1760 >>> f2 = fps(cos(x))

1761 >>> finv1, finv2 = f1.inverse(), f2.inverse()

1762

1763 >>> finv1._eval_terms(6)

1764 -x**5/120 + x**4/24 - x**3/6 + x**2/2 - x + 1

1765

1766 >>> finv2._eval_terms(8)

1767 61*x**6/720 + 5*x**4/24 + x**2/2 + 1

1768

1769 See Also

1770 ========

1771

1772 sympy.series.formal.FormalPowerSeries.inverse

1773 sympy.series.formal.FormalPowerSeries.coeff_bell

1774

1775 """

1776 ffps = self.ffps

1777 terms = [ffps.zero_coeff()]

1778

1779 for i in range(1, n):

1780 bell_seq = ffps.coeff_bell(i)

1781 seq = (self.aux_seq * bell_seq)

1782 terms.append(Add(*(seq[:i])) / ffps.fact_seq[i-1] * ffps.xk.coeff(i))

1783

1784 return Add(*terms)

1785

1786

1787def fps(f, x=None, x0=0, dir=1, hyper=True, order=4, rational=True, full=False):

1788 """

1789 Generates Formal Power Series of ``f``.

1790

1791 Explanation

1792 ===========

1793

1794 Returns the formal series expansion of ``f`` around ``x = x0``

1795 with respect to ``x`` in the form of a ``FormalPowerSeries`` object.

1796

1797 Formal Power Series is represented using an explicit formula

1798 computed using different algorithms.

1799

1800 See :func:`compute_fps` for the more details regarding the computation

1801 of formula.

1802

1803 Parameters

1804 ==========

1805

1806 x : Symbol, optional

1807 If x is None and ``f`` is univariate, the univariate symbols will be

1808 supplied, otherwise an error will be raised.

1809 x0 : number, optional

1810 Point to perform series expansion about. Default is 0.

1811 dir : {1, -1, '+', '-'}, optional

1812 If dir is 1 or '+' the series is calculated from the right and

1813 for -1 or '-' the series is calculated from the left. For smooth

1814 functions this flag will not alter the results. Default is 1.

1815 hyper : {True, False}, optional

1816 Set hyper to False to skip the hypergeometric algorithm.

1817 By default it is set to False.

1818 order : int, optional

1819 Order of the derivative of ``f``, Default is 4.

1820 rational : {True, False}, optional

1821 Set rational to False to skip rational algorithm. By default it is set

1822 to True.

1823 full : {True, False}, optional

1824 Set full to True to increase the range of rational algorithm.

1825 See :func:`rational_algorithm` for details. By default it is set to

1826 False.

1827

1828 Examples

1829 ========

1830

1831 >>> from sympy import fps, ln, atan, sin

1832 >>> from sympy.abc import x, n

1833

1834 Rational Functions

1835

1836 >>> fps(ln(1 + x)).truncate()

1837 x - x**2/2 + x**3/3 - x**4/4 + x**5/5 + O(x**6)

1838

1839 >>> fps(atan(x), full=True).truncate()

1840 x - x**3/3 + x**5/5 + O(x**6)

1841

1842 Symbolic Functions

1843

1844 >>> fps(x**n*sin(x**2), x).truncate(8)

1845 -x**(n + 6)/6 + x**(n + 2) + O(x**(n + 8))

1846

1847 See Also

1848 ========

1849

1850 sympy.series.formal.FormalPowerSeries

1851 sympy.series.formal.compute_fps

1852 """

1853 f = sympify(f)

1854

1855 if x is None:

1856 free = f.free_symbols

1857 if len(free) == 1:

1858 x = free.pop()

1859 elif not free:

1860 return f

1861 else:

1862 raise NotImplementedError("multivariate formal power series")

1863

1864 result = compute_fps(f, x, x0, dir, hyper, order, rational, full)

1865

1866 if result is None:

1867 return f

1868

1869 return FormalPowerSeries(f, x, x0, dir, result)