Coverage for /usr/lib/python3/dist-packages/sympy/polys/partfrac.py: 10%

1"""Algorithms for partial fraction decomposition of rational functions. """

4from sympy.core import S, Add, sympify, Function, Lambda, Dummy

5from sympy.core.traversal import preorder_traversal

6from sympy.polys import Poly, RootSum, cancel, factor

7from sympy.polys.polyerrors import PolynomialError

8from sympy.polys.polyoptions import allowed_flags, set_defaults

9from sympy.polys.polytools import parallel_poly_from_expr

10from sympy.utilities import numbered_symbols, take, xthreaded, public

13@xthreaded

14@public

15def apart(f, x=None, full=False, **options):

16 """

17 Compute partial fraction decomposition of a rational function.

19 Given a rational function ``f``, computes the partial fraction

20 decomposition of ``f``. Two algorithms are available: One is based on the

21 undertermined coefficients method, the other is Bronstein's full partial

22 fraction decomposition algorithm.

24 The undetermined coefficients method (selected by ``full=False``) uses

25 polynomial factorization (and therefore accepts the same options as

26 factor) for the denominator. Per default it works over the rational

27 numbers, therefore decomposition of denominators with non-rational roots

28 (e.g. irrational, complex roots) is not supported by default (see options

29 of factor).

31 Bronstein's algorithm can be selected by using ``full=True`` and allows a

32 decomposition of denominators with non-rational roots. A human-readable

33 result can be obtained via ``doit()`` (see examples below).

35 Examples

36 ========

38 >>> from sympy.polys.partfrac import apart

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

41 By default, using the undetermined coefficients method:

43 >>> apart(y/(x + 2)/(x + 1), x)

44 -y/(x + 2) + y/(x + 1)

46 The undetermined coefficients method does not provide a result when the

47 denominators roots are not rational:

49 >>> apart(y/(x**2 + x + 1), x)

50 y/(x**2 + x + 1)

52 You can choose Bronstein's algorithm by setting ``full=True``:

54 >>> apart(y/(x**2 + x + 1), x, full=True)

55 RootSum(_w**2 + _w + 1, Lambda(_a, (-2*_a*y/3 - y/3)/(-_a + x)))

57 Calling ``doit()`` yields a human-readable result:

59 >>> apart(y/(x**2 + x + 1), x, full=True).doit()

60 (-y/3 - 2*y*(-1/2 - sqrt(3)*I/2)/3)/(x + 1/2 + sqrt(3)*I/2) + (-y/3 -

61 2*y*(-1/2 + sqrt(3)*I/2)/3)/(x + 1/2 - sqrt(3)*I/2)

64 See Also

65 ========

67 apart_list, assemble_partfrac_list

68 """

69 allowed_flags(options, [])

71 f = sympify(f)

73 if f.is_Atom:

74 return f

75 else:

76 P, Q = f.as_numer_denom()

78 _options = options.copy()

79 options = set_defaults(options, extension=True)

80 try:

81 (P, Q), opt = parallel_poly_from_expr((P, Q), x, **options)

82 except PolynomialError as msg:

83 if f.is_commutative:

84 raise PolynomialError(msg)

85 # non-commutative

86 if f.is_Mul:

87 c, nc = f.args_cnc(split_1=False)

88 nc = f.func(*nc)

89 if c:

90 c = apart(f.func._from_args(c), x=x, full=full, **_options)

91 return c*nc

92 else:

93 return nc

94 elif f.is_Add:

95 c = []

96 nc = []

97 for i in f.args:

98 if i.is_commutative:

99 c.append(i)

100 else:

101 try:

102 nc.append(apart(i, x=x, full=full, **_options))

103 except NotImplementedError:

104 nc.append(i)

105 return apart(f.func(*c), x=x, full=full, **_options) + f.func(*nc)

106 else:

107 reps = []

108 pot = preorder_traversal(f)

109 next(pot)

110 for e in pot:

111 try:

112 reps.append((e, apart(e, x=x, full=full, **_options)))

113 pot.skip() # this was handled successfully

114 except NotImplementedError:

115 pass

116 return f.xreplace(dict(reps))

117

118 if P.is_multivariate:

119 fc = f.cancel()

120 if fc != f:

121 return apart(fc, x=x, full=full, **_options)

122

123 raise NotImplementedError(

124 "multivariate partial fraction decomposition")

125

126 common, P, Q = P.cancel(Q)

127

128 poly, P = P.div(Q, auto=True)

129 P, Q = P.rat_clear_denoms(Q)

130

131 if Q.degree() <= 1:

132 partial = P/Q

133 else:

134 if not full:

135 partial = apart_undetermined_coeffs(P, Q)

136 else:

137 partial = apart_full_decomposition(P, Q)

138

139 terms = S.Zero

140

141 for term in Add.make_args(partial):

142 if term.has(RootSum):

143 terms += term

144 else:

145 terms += factor(term)

146

147 return common*(poly.as_expr() + terms)

148

149

150def apart_undetermined_coeffs(P, Q):

151 """Partial fractions via method of undetermined coefficients. """

152 X = numbered_symbols(cls=Dummy)

153 partial, symbols = [], []

154

155 _, factors = Q.factor_list()

156

157 for f, k in factors:

158 n, q = f.degree(), Q

159

160 for i in range(1, k + 1):

161 coeffs, q = take(X, n), q.quo(f)

162 partial.append((coeffs, q, f, i))

163 symbols.extend(coeffs)

164

165 dom = Q.get_domain().inject(*symbols)

166 F = Poly(0, Q.gen, domain=dom)

167

168 for i, (coeffs, q, f, k) in enumerate(partial):

169 h = Poly(coeffs, Q.gen, domain=dom)

170 partial[i] = (h, f, k)

171 q = q.set_domain(dom)

172 F += h*q

173

174 system, result = [], S.Zero

175

176 for (k,), coeff in F.terms():

177 system.append(coeff - P.nth(k))

178

179 from sympy.solvers import solve

180 solution = solve(system, symbols)

181

182 for h, f, k in partial:

183 h = h.as_expr().subs(solution)

184 result += h/f.as_expr()**k

185

186 return result

187

188

189def apart_full_decomposition(P, Q):

190 """

191 Bronstein's full partial fraction decomposition algorithm.

192

193 Given a univariate rational function ``f``, performing only GCD

194 operations over the algebraic closure of the initial ground domain

195 of definition, compute full partial fraction decomposition with

196 fractions having linear denominators.

197

198 Note that no factorization of the initial denominator of ``f`` is

199 performed. The final decomposition is formed in terms of a sum of

200 :class:`RootSum` instances.

201

202 References

203 ==========

204

205 .. [1] [Bronstein93]_

206

207 """

208 return assemble_partfrac_list(apart_list(P/Q, P.gens[0]))

209

210

211@public

212def apart_list(f, x=None, dummies=None, **options):

213 """

214 Compute partial fraction decomposition of a rational function

215 and return the result in structured form.

216

217 Given a rational function ``f`` compute the partial fraction decomposition

218 of ``f``. Only Bronstein's full partial fraction decomposition algorithm

219 is supported by this method. The return value is highly structured and

220 perfectly suited for further algorithmic treatment rather than being

221 human-readable. The function returns a tuple holding three elements:

222

223 * The first item is the common coefficient, free of the variable `x` used

224 for decomposition. (It is an element of the base field `K`.)

225

226 * The second item is the polynomial part of the decomposition. This can be

227 the zero polynomial. (It is an element of `K[x]`.)

228

229 * The third part itself is a list of quadruples. Each quadruple

230 has the following elements in this order:

231

232 - The (not necessarily irreducible) polynomial `D` whose roots `w_i` appear

233 in the linear denominator of a bunch of related fraction terms. (This item

234 can also be a list of explicit roots. However, at the moment ``apart_list``

235 never returns a result this way, but the related ``assemble_partfrac_list``

236 function accepts this format as input.)

237

238 - The numerator of the fraction, written as a function of the root `w`

239

240 - The linear denominator of the fraction *excluding its power exponent*,

241 written as a function of the root `w`.

242

243 - The power to which the denominator has to be raised.

244

245 On can always rebuild a plain expression by using the function ``assemble_partfrac_list``.

246

247 Examples

248 ========

249

250 A first example:

251

252 >>> from sympy.polys.partfrac import apart_list, assemble_partfrac_list

253 >>> from sympy.abc import x, t

254

255 >>> f = (2*x**3 - 2*x) / (x**2 - 2*x + 1)

256 >>> pfd = apart_list(f)

257 >>> pfd

258 (1,

259 Poly(2*x + 4, x, domain='ZZ'),

260 [(Poly(_w - 1, _w, domain='ZZ'), Lambda(_a, 4), Lambda(_a, -_a + x), 1)])

261

262 >>> assemble_partfrac_list(pfd)

263 2*x + 4 + 4/(x - 1)

264

265 Second example:

266

267 >>> f = (-2*x - 2*x**2) / (3*x**2 - 6*x)

268 >>> pfd = apart_list(f)

269 >>> pfd

270 (-1,

271 Poly(2/3, x, domain='QQ'),

272 [(Poly(_w - 2, _w, domain='ZZ'), Lambda(_a, 2), Lambda(_a, -_a + x), 1)])

273

274 >>> assemble_partfrac_list(pfd)

275 -2/3 - 2/(x - 2)

276

277 Another example, showing symbolic parameters:

278

279 >>> pfd = apart_list(t/(x**2 + x + t), x)

280 >>> pfd

281 (1,

282 Poly(0, x, domain='ZZ[t]'),

283 [(Poly(_w**2 + _w + t, _w, domain='ZZ[t]'),

284 Lambda(_a, -2*_a*t/(4*t - 1) - t/(4*t - 1)),

285 Lambda(_a, -_a + x),

286 1)])

287

288 >>> assemble_partfrac_list(pfd)

289 RootSum(_w**2 + _w + t, Lambda(_a, (-2*_a*t/(4*t - 1) - t/(4*t - 1))/(-_a + x)))

290

291 This example is taken from Bronstein's original paper:

292

293 >>> f = 36 / (x**5 - 2*x**4 - 2*x**3 + 4*x**2 + x - 2)

294 >>> pfd = apart_list(f)

295 >>> pfd

296 (1,

297 Poly(0, x, domain='ZZ'),

298 [(Poly(_w - 2, _w, domain='ZZ'), Lambda(_a, 4), Lambda(_a, -_a + x), 1),

299 (Poly(_w**2 - 1, _w, domain='ZZ'), Lambda(_a, -3*_a - 6), Lambda(_a, -_a + x), 2),

300 (Poly(_w + 1, _w, domain='ZZ'), Lambda(_a, -4), Lambda(_a, -_a + x), 1)])

301

302 >>> assemble_partfrac_list(pfd)

303 -4/(x + 1) - 3/(x + 1)**2 - 9/(x - 1)**2 + 4/(x - 2)

304

305 See also

306 ========

307

308 apart, assemble_partfrac_list

309

310 References

311 ==========

312

313 .. [1] [Bronstein93]_

314

315 """

316 allowed_flags(options, [])

317

318 f = sympify(f)

319

320 if f.is_Atom:

321 return f

322 else:

323 P, Q = f.as_numer_denom()

324

325 options = set_defaults(options, extension=True)

326 (P, Q), opt = parallel_poly_from_expr((P, Q), x, **options)

327

328 if P.is_multivariate:

329 raise NotImplementedError(

330 "multivariate partial fraction decomposition")

331

332 common, P, Q = P.cancel(Q)

333

334 poly, P = P.div(Q, auto=True)

335 P, Q = P.rat_clear_denoms(Q)

336

337 polypart = poly

338

339 if dummies is None:

340 def dummies(name):

341 d = Dummy(name)

342 while True:

343 yield d

344

345 dummies = dummies("w")

346

347 rationalpart = apart_list_full_decomposition(P, Q, dummies)

348

349 return (common, polypart, rationalpart)

350

351

352def apart_list_full_decomposition(P, Q, dummygen):

353 """

354 Bronstein's full partial fraction decomposition algorithm.

355

356 Given a univariate rational function ``f``, performing only GCD

357 operations over the algebraic closure of the initial ground domain

358 of definition, compute full partial fraction decomposition with

359 fractions having linear denominators.

360

361 Note that no factorization of the initial denominator of ``f`` is

362 performed. The final decomposition is formed in terms of a sum of

363 :class:`RootSum` instances.

364

365 References

366 ==========

367

368 .. [1] [Bronstein93]_

369

370 """

371 f, x, U = P/Q, P.gen, []

372

373 u = Function('u')(x)

374 a = Dummy('a')

375

376 partial = []

377

378 for d, n in Q.sqf_list_include(all=True):

379 b = d.as_expr()

380 U += [ u.diff(x, n - 1) ]

381

382 h = cancel(f*b**n) / u**n

383

384 H, subs = [h], []

385

386 for j in range(1, n):

387 H += [ H[-1].diff(x) / j ]

388

389 for j in range(1, n + 1):

390 subs += [ (U[j - 1], b.diff(x, j) / j) ]

391

392 for j in range(0, n):

393 P, Q = cancel(H[j]).as_numer_denom()

394

395 for i in range(0, j + 1):

396 P = P.subs(*subs[j - i])

397

398 Q = Q.subs(*subs[0])

399

400 P = Poly(P, x)

401 Q = Poly(Q, x)

402

403 G = P.gcd(d)

404 D = d.quo(G)

405

406 B, g = Q.half_gcdex(D)

407 b = (P * B.quo(g)).rem(D)

408

409 Dw = D.subs(x, next(dummygen))

410 numer = Lambda(a, b.as_expr().subs(x, a))

411 denom = Lambda(a, (x - a))

412 exponent = n-j

413

414 partial.append((Dw, numer, denom, exponent))

415

416 return partial

417

418

419@public

420def assemble_partfrac_list(partial_list):

421 r"""Reassemble a full partial fraction decomposition

422 from a structured result obtained by the function ``apart_list``.

423

424 Examples

425 ========

426

427 This example is taken from Bronstein's original paper:

428

429 >>> from sympy.polys.partfrac import apart_list, assemble_partfrac_list

430 >>> from sympy.abc import x

431

432 >>> f = 36 / (x**5 - 2*x**4 - 2*x**3 + 4*x**2 + x - 2)

433 >>> pfd = apart_list(f)

434 >>> pfd

435 (1,

436 Poly(0, x, domain='ZZ'),

437 [(Poly(_w - 2, _w, domain='ZZ'), Lambda(_a, 4), Lambda(_a, -_a + x), 1),

438 (Poly(_w**2 - 1, _w, domain='ZZ'), Lambda(_a, -3*_a - 6), Lambda(_a, -_a + x), 2),

439 (Poly(_w + 1, _w, domain='ZZ'), Lambda(_a, -4), Lambda(_a, -_a + x), 1)])

440

441 >>> assemble_partfrac_list(pfd)

442 -4/(x + 1) - 3/(x + 1)**2 - 9/(x - 1)**2 + 4/(x - 2)

443

444 If we happen to know some roots we can provide them easily inside the structure:

445

446 >>> pfd = apart_list(2/(x**2-2))

447 >>> pfd

448 (1,

449 Poly(0, x, domain='ZZ'),

450 [(Poly(_w**2 - 2, _w, domain='ZZ'),

451 Lambda(_a, _a/2),

452 Lambda(_a, -_a + x),

453 1)])

454

455 >>> pfda = assemble_partfrac_list(pfd)

456 >>> pfda

457 RootSum(_w**2 - 2, Lambda(_a, _a/(-_a + x)))/2

458

459 >>> pfda.doit()

460 -sqrt(2)/(2*(x + sqrt(2))) + sqrt(2)/(2*(x - sqrt(2)))

461

462 >>> from sympy import Dummy, Poly, Lambda, sqrt

463 >>> a = Dummy("a")

464 >>> pfd = (1, Poly(0, x, domain='ZZ'), [([sqrt(2),-sqrt(2)], Lambda(a, a/2), Lambda(a, -a + x), 1)])

465

466 >>> assemble_partfrac_list(pfd)

467 -sqrt(2)/(2*(x + sqrt(2))) + sqrt(2)/(2*(x - sqrt(2)))

468

469 See Also

470 ========

471

472 apart, apart_list

473 """

474 # Common factor

475 common = partial_list[0]

476

477 # Polynomial part

478 polypart = partial_list[1]

479 pfd = polypart.as_expr()

480

481 # Rational parts

482 for r, nf, df, ex in partial_list[2]:

483 if isinstance(r, Poly):

484 # Assemble in case the roots are given implicitly by a polynomials

485 an, nu = nf.variables, nf.expr

486 ad, de = df.variables, df.expr

487 # Hack to make dummies equal because Lambda created new Dummies

488 de = de.subs(ad[0], an[0])

489 func = Lambda(tuple(an), nu/de**ex)

490 pfd += RootSum(r, func, auto=False, quadratic=False)

491 else:

492 # Assemble in case the roots are given explicitly by a list of algebraic numbers

493 for root in r:

494 pfd += nf(root)/df(root)**ex

495

496 return common*pfd