Coverage for /usr/lib/python3/dist-packages/sympy/integrals/rationaltools.py: 9%

1"""This module implements tools for integrating rational functions. """

3from sympy.core.function import Lambda

4from sympy.core.numbers import I

5from sympy.core.singleton import S

6from sympy.core.symbol import (Dummy, Symbol, symbols)

7from sympy.functions.elementary.exponential import log

8from sympy.functions.elementary.trigonometric import atan

9from sympy.polys.polyroots import roots

10from sympy.polys.polytools import cancel

11from sympy.polys.rootoftools import RootSum

12from sympy.polys import Poly, resultant, ZZ

15def ratint(f, x, **flags):

16 """

17 Performs indefinite integration of rational functions.

19 Explanation

20 ===========

22 Given a field :math:`K` and a rational function :math:`f = p/q`,

23 where :math:`p` and :math:`q` are polynomials in :math:`K[x]`,

24 returns a function :math:`g` such that :math:`f = g'`.

26 Examples

27 ========

29 >>> from sympy.integrals.rationaltools import ratint

30 >>> from sympy.abc import x

32 >>> ratint(36/(x**5 - 2*x**4 - 2*x**3 + 4*x**2 + x - 2), x)

33 (12*x + 6)/(x**2 - 1) + 4*log(x - 2) - 4*log(x + 1)

35 References

36 ==========

38 .. [1] M. Bronstein, Symbolic Integration I: Transcendental

39 Functions, Second Edition, Springer-Verlag, 2005, pp. 35-70

41 See Also

42 ========

44 sympy.integrals.integrals.Integral.doit

45 sympy.integrals.rationaltools.ratint_logpart

46 sympy.integrals.rationaltools.ratint_ratpart

48 """

49 if isinstance(f, tuple):

50 p, q = f

51 else:

52 p, q = f.as_numer_denom()

54 p, q = Poly(p, x, composite=False, field=True), Poly(q, x, composite=False, field=True)

56 coeff, p, q = p.cancel(q)

57 poly, p = p.div(q)

59 result = poly.integrate(x).as_expr()

61 if p.is_zero:

62 return coeff*result

64 g, h = ratint_ratpart(p, q, x)

66 P, Q = h.as_numer_denom()

68 P = Poly(P, x)

69 Q = Poly(Q, x)

71 q, r = P.div(Q)

73 result += g + q.integrate(x).as_expr()

75 if not r.is_zero:

76 symbol = flags.get('symbol', 't')

78 if not isinstance(symbol, Symbol):

79 t = Dummy(symbol)

80 else:

81 t = symbol.as_dummy()

83 L = ratint_logpart(r, Q, x, t)

85 real = flags.get('real')

87 if real is None:

88 if isinstance(f, tuple):

89 p, q = f

90 atoms = p.atoms() | q.atoms()

91 else:

92 atoms = f.atoms()

94 for elt in atoms - {x}:

95 if not elt.is_extended_real:

96 real = False

97 break

98 else:

99 real = True

100

101 eps = S.Zero

102

103 if not real:

104 for h, q in L:

105 _, h = h.primitive()

106 eps += RootSum(

107 q, Lambda(t, t*log(h.as_expr())), quadratic=True)

108 else:

109 for h, q in L:

110 _, h = h.primitive()

111 R = log_to_real(h, q, x, t)

112

113 if R is not None:

114 eps += R

115 else:

116 eps += RootSum(

117 q, Lambda(t, t*log(h.as_expr())), quadratic=True)

118

119 result += eps

120

121 return coeff*result

122

123

124def ratint_ratpart(f, g, x):

125 """

126 Horowitz-Ostrogradsky algorithm.

127

128 Explanation

129 ===========

130

131 Given a field K and polynomials f and g in K[x], such that f and g

132 are coprime and deg(f) < deg(g), returns fractions A and B in K(x),

133 such that f/g = A' + B and B has square-free denominator.

134

135 Examples

136 ========

137

138 >>> from sympy.integrals.rationaltools import ratint_ratpart

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

140 >>> from sympy import Poly

141 >>> ratint_ratpart(Poly(1, x, domain='ZZ'),

142 ... Poly(x + 1, x, domain='ZZ'), x)

143 (0, 1/(x + 1))

144 >>> ratint_ratpart(Poly(1, x, domain='EX'),

145 ... Poly(x**2 + y**2, x, domain='EX'), x)

146 (0, 1/(x**2 + y**2))

147 >>> ratint_ratpart(Poly(36, x, domain='ZZ'),

148 ... Poly(x**5 - 2*x**4 - 2*x**3 + 4*x**2 + x - 2, x, domain='ZZ'), x)

149 ((12*x + 6)/(x**2 - 1), 12/(x**2 - x - 2))

150

151 See Also

152 ========

153

154 ratint, ratint_logpart

155 """

156 from sympy.solvers.solvers import solve

157

158 f = Poly(f, x)

159 g = Poly(g, x)

160

161 u, v, _ = g.cofactors(g.diff())

162

163 n = u.degree()

164 m = v.degree()

165

166 A_coeffs = [ Dummy('a' + str(n - i)) for i in range(0, n) ]

167 B_coeffs = [ Dummy('b' + str(m - i)) for i in range(0, m) ]

168

169 C_coeffs = A_coeffs + B_coeffs

170

171 A = Poly(A_coeffs, x, domain=ZZ[C_coeffs])

172 B = Poly(B_coeffs, x, domain=ZZ[C_coeffs])

173

174 H = f - A.diff()*v + A*(u.diff()*v).quo(u) - B*u

175

176 result = solve(H.coeffs(), C_coeffs)

177

178 A = A.as_expr().subs(result)

179 B = B.as_expr().subs(result)

180

181 rat_part = cancel(A/u.as_expr(), x)

182 log_part = cancel(B/v.as_expr(), x)

183

184 return rat_part, log_part

185

186

187def ratint_logpart(f, g, x, t=None):

188 r"""

189 Lazard-Rioboo-Trager algorithm.

190

191 Explanation

192 ===========

193

194 Given a field K and polynomials f and g in K[x], such that f and g

195 are coprime, deg(f) < deg(g) and g is square-free, returns a list

196 of tuples (s_i, q_i) of polynomials, for i = 1..n, such that s_i

197 in K[t, x] and q_i in K[t], and::

198

199 ___ ___

200 d f d \ ` \ `

201 -- - = -- ) ) a log(s_i(a, x))

202 dx g dx /__, /__,

203 i=1..n a | q_i(a) = 0

204

205 Examples

206 ========

207

208 >>> from sympy.integrals.rationaltools import ratint_logpart

209 >>> from sympy.abc import x

210 >>> from sympy import Poly

211 >>> ratint_logpart(Poly(1, x, domain='ZZ'),

212 ... Poly(x**2 + x + 1, x, domain='ZZ'), x)

213 [(Poly(x + 3*_t/2 + 1/2, x, domain='QQ[_t]'),

214 ...Poly(3*_t**2 + 1, _t, domain='ZZ'))]

215 >>> ratint_logpart(Poly(12, x, domain='ZZ'),

216 ... Poly(x**2 - x - 2, x, domain='ZZ'), x)

217 [(Poly(x - 3*_t/8 - 1/2, x, domain='QQ[_t]'),

218 ...Poly(-_t**2 + 16, _t, domain='ZZ'))]

219

220 See Also

221 ========

222

223 ratint, ratint_ratpart

224 """

225 f, g = Poly(f, x), Poly(g, x)

226

227 t = t or Dummy('t')

228 a, b = g, f - g.diff()*Poly(t, x)

229

230 res, R = resultant(a, b, includePRS=True)

231 res = Poly(res, t, composite=False)

232

233 assert res, "BUG: resultant(%s, %s) cannot be zero" % (a, b)

234

235 R_map, H = {}, []

236

237 for r in R:

238 R_map[r.degree()] = r

239

240 def _include_sign(c, sqf):

241 if c.is_extended_real and (c < 0) == True:

242 h, k = sqf[0]

243 c_poly = c.as_poly(h.gens)

244 sqf[0] = h*c_poly, k

245

246 C, res_sqf = res.sqf_list()

247 _include_sign(C, res_sqf)

248

249 for q, i in res_sqf:

250 _, q = q.primitive()

251

252 if g.degree() == i:

253 H.append((g, q))

254 else:

255 h = R_map[i]

256 h_lc = Poly(h.LC(), t, field=True)

257

258 c, h_lc_sqf = h_lc.sqf_list(all=True)

259 _include_sign(c, h_lc_sqf)

260

261 for a, j in h_lc_sqf:

262 h = h.quo(Poly(a.gcd(q)**j, x))

263

264 inv, coeffs = h_lc.invert(q), [S.One]

265

266 for coeff in h.coeffs()[1:]:

267 coeff = coeff.as_poly(inv.gens)

268 T = (inv*coeff).rem(q)

269 coeffs.append(T.as_expr())

270

271 h = Poly(dict(list(zip(h.monoms(), coeffs))), x)

273 H.append((h, q))

275 return H

278def log_to_atan(f, g):

279 """

280 Convert complex logarithms to real arctangents.

281

282 Explanation

283 ===========

284

285 Given a real field K and polynomials f and g in K[x], with g != 0,

286 returns a sum h of arctangents of polynomials in K[x], such that:

287

288 dh d f + I g

289 -- = -- I log( ------- )

290 dx dx f - I g

291

292 Examples

293 ========

294

295 >>> from sympy.integrals.rationaltools import log_to_atan

296 >>> from sympy.abc import x

297 >>> from sympy import Poly, sqrt, S

298 >>> log_to_atan(Poly(x, x, domain='ZZ'), Poly(1, x, domain='ZZ'))

299 2*atan(x)

300 >>> log_to_atan(Poly(x + S(1)/2, x, domain='QQ'),

301 ... Poly(sqrt(3)/2, x, domain='EX'))

302 2*atan(2*sqrt(3)*x/3 + sqrt(3)/3)

303

304 See Also

305 ========

306

307 log_to_real

308 """

309 if f.degree() < g.degree():

310 f, g = -g, f

311

312 f = f.to_field()

313 g = g.to_field()

314

315 p, q = f.div(g)

316

317 if q.is_zero:

318 return 2*atan(p.as_expr())

319 else:

320 s, t, h = g.gcdex(-f)

321 u = (f*s + g*t).quo(h)

322 A = 2*atan(u.as_expr())

323

324 return A + log_to_atan(s, t)

325

326

327def log_to_real(h, q, x, t):

328 r"""

329 Convert complex logarithms to real functions.

330

331 Explanation

332 ===========

333

334 Given real field K and polynomials h in K[t,x] and q in K[t],

335 returns real function f such that:

336 ___

337 df d \ `

338 -- = -- ) a log(h(a, x))

339 dx dx /__,

340 a | q(a) = 0

341

342 Examples

343 ========

344

345 >>> from sympy.integrals.rationaltools import log_to_real

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

347 >>> from sympy import Poly, S

348 >>> log_to_real(Poly(x + 3*y/2 + S(1)/2, x, domain='QQ[y]'),

349 ... Poly(3*y**2 + 1, y, domain='ZZ'), x, y)

350 2*sqrt(3)*atan(2*sqrt(3)*x/3 + sqrt(3)/3)/3

351 >>> log_to_real(Poly(x**2 - 1, x, domain='ZZ'),

352 ... Poly(-2*y + 1, y, domain='ZZ'), x, y)

353 log(x**2 - 1)/2

354

355 See Also

356 ========

357

358 log_to_atan

359 """

360 from sympy.simplify.radsimp import collect

361 u, v = symbols('u,v', cls=Dummy)

362

363 H = h.as_expr().subs({t: u + I*v}).expand()

364 Q = q.as_expr().subs({t: u + I*v}).expand()

365

366 H_map = collect(H, I, evaluate=False)

367 Q_map = collect(Q, I, evaluate=False)

368

369 a, b = H_map.get(S.One, S.Zero), H_map.get(I, S.Zero)

370 c, d = Q_map.get(S.One, S.Zero), Q_map.get(I, S.Zero)

371

372 R = Poly(resultant(c, d, v), u)

373

374 R_u = roots(R, filter='R')

375

376 if len(R_u) != R.count_roots():

377 return None

378

379 result = S.Zero

380

381 for r_u in R_u.keys():

382 C = Poly(c.subs({u: r_u}), v)

383 R_v = roots(C, filter='R')

384

385 if len(R_v) != C.count_roots():

386 return None

387

388 R_v_paired = [] # take one from each pair of conjugate roots

389 for r_v in R_v:

390 if r_v not in R_v_paired and -r_v not in R_v_paired:

391 if r_v.is_negative or r_v.could_extract_minus_sign():

392 R_v_paired.append(-r_v)

393 elif not r_v.is_zero:

394 R_v_paired.append(r_v)

395

396 for r_v in R_v_paired:

397

398 D = d.subs({u: r_u, v: r_v})

399

400 if D.evalf(chop=True) != 0:

401 continue

402

403 A = Poly(a.subs({u: r_u, v: r_v}), x)

404 B = Poly(b.subs({u: r_u, v: r_v}), x)

405

406 AB = (A**2 + B**2).as_expr()

407

408 result += r_u*log(AB) + r_v*log_to_atan(A, B)

409

410 R_q = roots(q, filter='R')

411

412 if len(R_q) != q.count_roots():

413 return None

414

415 for r in R_q.keys():

416 result += r*log(h.as_expr().subs(t, r))

417

418 return result