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

1"""Minimal polynomials for algebraic numbers.""" 

2 

3from functools import reduce 

4 

5from sympy.core.add import Add 

6from sympy.core.exprtools import Factors 

7from sympy.core.function import expand_mul, expand_multinomial, _mexpand 

8from sympy.core.mul import Mul 

9from sympy.core.numbers import (I, Rational, pi, _illegal) 

10from sympy.core.singleton import S 

11from sympy.core.symbol import Dummy 

12from sympy.core.sympify import sympify 

13from sympy.core.traversal import preorder_traversal 

14from sympy.functions.elementary.exponential import exp 

15from sympy.functions.elementary.miscellaneous import sqrt, cbrt 

16from sympy.functions.elementary.trigonometric import cos, sin, tan 

17from sympy.ntheory.factor_ import divisors 

18from sympy.utilities.iterables import subsets 

19 

20from sympy.polys.domains import ZZ, QQ, FractionField 

21from sympy.polys.orthopolys import dup_chebyshevt 

22from sympy.polys.polyerrors import ( 

23 NotAlgebraic, 

24 GeneratorsError, 

25) 

26from sympy.polys.polytools import ( 

27 Poly, PurePoly, invert, factor_list, groebner, resultant, 

28 degree, poly_from_expr, parallel_poly_from_expr, lcm 

29) 

30from sympy.polys.polyutils import dict_from_expr, expr_from_dict 

31from sympy.polys.ring_series import rs_compose_add 

32from sympy.polys.rings import ring 

33from sympy.polys.rootoftools import CRootOf 

34from sympy.polys.specialpolys import cyclotomic_poly 

35from sympy.utilities import ( 

36 numbered_symbols, public, sift 

37) 

38 

39 

40def _choose_factor(factors, x, v, dom=QQ, prec=200, bound=5): 

41 """ 

42 Return a factor having root ``v`` 

43 It is assumed that one of the factors has root ``v``. 

44 """ 

45 

46 if isinstance(factors[0], tuple): 

47 factors = [f[0] for f in factors] 

48 if len(factors) == 1: 

49 return factors[0] 

50 

51 prec1 = 10 

52 points = {} 

53 symbols = dom.symbols if hasattr(dom, 'symbols') else [] 

54 while prec1 <= prec: 

55 # when dealing with non-Rational numbers we usually evaluate 

56 # with `subs` argument but we only need a ballpark evaluation 

57 fe = [f.as_expr().xreplace({x:v}) for f in factors] 

58 if v.is_number: 

59 fe = [f.n(prec) for f in fe] 

60 

61 # assign integers [0, n) to symbols (if any) 

62 for n in subsets(range(bound), k=len(symbols), repetition=True): 

63 for s, i in zip(symbols, n): 

64 points[s] = i 

65 

66 # evaluate the expression at these points 

67 candidates = [(abs(f.subs(points).n(prec1)), i) 

68 for i,f in enumerate(fe)] 

69 

70 # if we get invalid numbers (e.g. from division by zero) 

71 # we try again 

72 if any(i in _illegal for i, _ in candidates): 

73 continue 

74 

75 # find the smallest two -- if they differ significantly 

76 # then we assume we have found the factor that becomes 

77 # 0 when v is substituted into it 

78 can = sorted(candidates) 

79 (a, ix), (b, _) = can[:2] 

80 if b > a * 10**6: # XXX what to use? 

81 return factors[ix] 

82 

83 prec1 *= 2 

84 

85 raise NotImplementedError("multiple candidates for the minimal polynomial of %s" % v) 

86 

87 

88def _is_sum_surds(p): 

89 args = p.args if p.is_Add else [p] 

90 for y in args: 

91 if not ((y**2).is_Rational and y.is_extended_real): 

92 return False 

93 return True 

94 

95 

96def _separate_sq(p): 

97 """ 

98 helper function for ``_minimal_polynomial_sq`` 

99 

100 It selects a rational ``g`` such that the polynomial ``p`` 

101 consists of a sum of terms whose surds squared have gcd equal to ``g`` 

102 and a sum of terms with surds squared prime with ``g``; 

103 then it takes the field norm to eliminate ``sqrt(g)`` 

104 

105 See simplify.simplify.split_surds and polytools.sqf_norm. 

106 

107 Examples 

108 ======== 

109 

110 >>> from sympy import sqrt 

111 >>> from sympy.abc import x 

112 >>> from sympy.polys.numberfields.minpoly import _separate_sq 

113 >>> p= -x + sqrt(2) + sqrt(3) + sqrt(7) 

114 >>> p = _separate_sq(p); p 

115 -x**2 + 2*sqrt(3)*x + 2*sqrt(7)*x - 2*sqrt(21) - 8 

116 >>> p = _separate_sq(p); p 

117 -x**4 + 4*sqrt(7)*x**3 - 32*x**2 + 8*sqrt(7)*x + 20 

118 >>> p = _separate_sq(p); p 

119 -x**8 + 48*x**6 - 536*x**4 + 1728*x**2 - 400 

120 

121 """ 

122 def is_sqrt(expr): 

123 return expr.is_Pow and expr.exp is S.Half 

124 # p = c1*sqrt(q1) + ... + cn*sqrt(qn) -> a = [(c1, q1), .., (cn, qn)] 

125 a = [] 

126 for y in p.args: 

127 if not y.is_Mul: 

128 if is_sqrt(y): 

129 a.append((S.One, y**2)) 

130 elif y.is_Atom: 

131 a.append((y, S.One)) 

132 elif y.is_Pow and y.exp.is_integer: 

133 a.append((y, S.One)) 

134 else: 

135 raise NotImplementedError 

136 else: 

137 T, F = sift(y.args, is_sqrt, binary=True) 

138 a.append((Mul(*F), Mul(*T)**2)) 

139 a.sort(key=lambda z: z[1]) 

140 if a[-1][1] is S.One: 

141 # there are no surds 

142 return p 

143 surds = [z for y, z in a] 

144 for i in range(len(surds)): 

145 if surds[i] != 1: 

146 break 

147 from sympy.simplify.radsimp import _split_gcd 

148 g, b1, b2 = _split_gcd(*surds[i:]) 

149 a1 = [] 

150 a2 = [] 

151 for y, z in a: 

152 if z in b1: 

153 a1.append(y*z**S.Half) 

154 else: 

155 a2.append(y*z**S.Half) 

156 p1 = Add(*a1) 

157 p2 = Add(*a2) 

158 p = _mexpand(p1**2) - _mexpand(p2**2) 

159 return p 

160 

161def _minimal_polynomial_sq(p, n, x): 

162 """ 

163 Returns the minimal polynomial for the ``nth-root`` of a sum of surds 

164 or ``None`` if it fails. 

165 

166 Parameters 

167 ========== 

168 

169 p : sum of surds 

170 n : positive integer 

171 x : variable of the returned polynomial 

172 

173 Examples 

174 ======== 

175 

176 >>> from sympy.polys.numberfields.minpoly import _minimal_polynomial_sq 

177 >>> from sympy import sqrt 

178 >>> from sympy.abc import x 

179 >>> q = 1 + sqrt(2) + sqrt(3) 

180 >>> _minimal_polynomial_sq(q, 3, x) 

181 x**12 - 4*x**9 - 4*x**6 + 16*x**3 - 8 

182 

183 """ 

184 p = sympify(p) 

185 n = sympify(n) 

186 if not n.is_Integer or not n > 0 or not _is_sum_surds(p): 

187 return None 

188 pn = p**Rational(1, n) 

189 # eliminate the square roots 

190 p -= x 

191 while 1: 

192 p1 = _separate_sq(p) 

193 if p1 is p: 

194 p = p1.subs({x:x**n}) 

195 break 

196 else: 

197 p = p1 

198 

199 # _separate_sq eliminates field extensions in a minimal way, so that 

200 # if n = 1 then `p = constant*(minimal_polynomial(p))` 

201 # if n > 1 it contains the minimal polynomial as a factor. 

202 if n == 1: 

203 p1 = Poly(p) 

204 if p.coeff(x**p1.degree(x)) < 0: 

205 p = -p 

206 p = p.primitive()[1] 

207 return p 

208 # by construction `p` has root `pn` 

209 # the minimal polynomial is the factor vanishing in x = pn 

210 factors = factor_list(p)[1] 

211 

212 result = _choose_factor(factors, x, pn) 

213 return result 

214 

215def _minpoly_op_algebraic_element(op, ex1, ex2, x, dom, mp1=None, mp2=None): 

216 """ 

217 return the minimal polynomial for ``op(ex1, ex2)`` 

218 

219 Parameters 

220 ========== 

221 

222 op : operation ``Add`` or ``Mul`` 

223 ex1, ex2 : expressions for the algebraic elements 

224 x : indeterminate of the polynomials 

225 dom: ground domain 

226 mp1, mp2 : minimal polynomials for ``ex1`` and ``ex2`` or None 

227 

228 Examples 

229 ======== 

230 

231 >>> from sympy import sqrt, Add, Mul, QQ 

232 >>> from sympy.polys.numberfields.minpoly import _minpoly_op_algebraic_element 

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

234 >>> p1 = sqrt(sqrt(2) + 1) 

235 >>> p2 = sqrt(sqrt(2) - 1) 

236 >>> _minpoly_op_algebraic_element(Mul, p1, p2, x, QQ) 

237 x - 1 

238 >>> q1 = sqrt(y) 

239 >>> q2 = 1 / y 

240 >>> _minpoly_op_algebraic_element(Add, q1, q2, x, QQ.frac_field(y)) 

241 x**2*y**2 - 2*x*y - y**3 + 1 

242 

243 References 

244 ========== 

245 

246 .. [1] https://en.wikipedia.org/wiki/Resultant 

247 .. [2] I.M. Isaacs, Proc. Amer. Math. Soc. 25 (1970), 638 

248 "Degrees of sums in a separable field extension". 

249 

250 """ 

251 y = Dummy(str(x)) 

252 if mp1 is None: 

253 mp1 = _minpoly_compose(ex1, x, dom) 

254 if mp2 is None: 

255 mp2 = _minpoly_compose(ex2, y, dom) 

256 else: 

257 mp2 = mp2.subs({x: y}) 

258 

259 if op is Add: 

260 # mp1a = mp1.subs({x: x - y}) 

261 if dom == QQ: 

262 R, X = ring('X', QQ) 

263 p1 = R(dict_from_expr(mp1)[0]) 

264 p2 = R(dict_from_expr(mp2)[0]) 

265 else: 

266 (p1, p2), _ = parallel_poly_from_expr((mp1, x - y), x, y) 

267 r = p1.compose(p2) 

268 mp1a = r.as_expr() 

269 

270 elif op is Mul: 

271 mp1a = _muly(mp1, x, y) 

272 else: 

273 raise NotImplementedError('option not available') 

274 

275 if op is Mul or dom != QQ: 

276 r = resultant(mp1a, mp2, gens=[y, x]) 

277 else: 

278 r = rs_compose_add(p1, p2) 

279 r = expr_from_dict(r.as_expr_dict(), x) 

280 

281 deg1 = degree(mp1, x) 

282 deg2 = degree(mp2, y) 

283 if op is Mul and deg1 == 1 or deg2 == 1: 

284 # if deg1 = 1, then mp1 = x - a; mp1a = x - y - a; 

285 # r = mp2(x - a), so that `r` is irreducible 

286 return r 

287 

288 r = Poly(r, x, domain=dom) 

289 _, factors = r.factor_list() 

290 res = _choose_factor(factors, x, op(ex1, ex2), dom) 

291 return res.as_expr() 

292 

293 

294def _invertx(p, x): 

295 """ 

296 Returns ``expand_mul(x**degree(p, x)*p.subs(x, 1/x))`` 

297 """ 

298 p1 = poly_from_expr(p, x)[0] 

299 

300 n = degree(p1) 

301 a = [c * x**(n - i) for (i,), c in p1.terms()] 

302 return Add(*a) 

303 

304 

305def _muly(p, x, y): 

306 """ 

307 Returns ``_mexpand(y**deg*p.subs({x:x / y}))`` 

308 """ 

309 p1 = poly_from_expr(p, x)[0] 

310 

311 n = degree(p1) 

312 a = [c * x**i * y**(n - i) for (i,), c in p1.terms()] 

313 return Add(*a) 

314 

315 

316def _minpoly_pow(ex, pw, x, dom, mp=None): 

317 """ 

318 Returns ``minpoly(ex**pw, x)`` 

319 

320 Parameters 

321 ========== 

322 

323 ex : algebraic element 

324 pw : rational number 

325 x : indeterminate of the polynomial 

326 dom: ground domain 

327 mp : minimal polynomial of ``p`` 

328 

329 Examples 

330 ======== 

331 

332 >>> from sympy import sqrt, QQ, Rational 

333 >>> from sympy.polys.numberfields.minpoly import _minpoly_pow, minpoly 

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

335 >>> p = sqrt(1 + sqrt(2)) 

336 >>> _minpoly_pow(p, 2, x, QQ) 

337 x**2 - 2*x - 1 

338 >>> minpoly(p**2, x) 

339 x**2 - 2*x - 1 

340 >>> _minpoly_pow(y, Rational(1, 3), x, QQ.frac_field(y)) 

341 x**3 - y 

342 >>> minpoly(y**Rational(1, 3), x) 

343 x**3 - y 

344 

345 """ 

346 pw = sympify(pw) 

347 if not mp: 

348 mp = _minpoly_compose(ex, x, dom) 

349 if not pw.is_rational: 

350 raise NotAlgebraic("%s does not seem to be an algebraic element" % ex) 

351 if pw < 0: 

352 if mp == x: 

353 raise ZeroDivisionError('%s is zero' % ex) 

354 mp = _invertx(mp, x) 

355 if pw == -1: 

356 return mp 

357 pw = -pw 

358 ex = 1/ex 

359 

360 y = Dummy(str(x)) 

361 mp = mp.subs({x: y}) 

362 n, d = pw.as_numer_denom() 

363 res = Poly(resultant(mp, x**d - y**n, gens=[y]), x, domain=dom) 

364 _, factors = res.factor_list() 

365 res = _choose_factor(factors, x, ex**pw, dom) 

366 return res.as_expr() 

367 

368 

369def _minpoly_add(x, dom, *a): 

370 """ 

371 returns ``minpoly(Add(*a), dom, x)`` 

372 """ 

373 mp = _minpoly_op_algebraic_element(Add, a[0], a[1], x, dom) 

374 p = a[0] + a[1] 

375 for px in a[2:]: 

376 mp = _minpoly_op_algebraic_element(Add, p, px, x, dom, mp1=mp) 

377 p = p + px 

378 return mp 

379 

380 

381def _minpoly_mul(x, dom, *a): 

382 """ 

383 returns ``minpoly(Mul(*a), dom, x)`` 

384 """ 

385 mp = _minpoly_op_algebraic_element(Mul, a[0], a[1], x, dom) 

386 p = a[0] * a[1] 

387 for px in a[2:]: 

388 mp = _minpoly_op_algebraic_element(Mul, p, px, x, dom, mp1=mp) 

389 p = p * px 

390 return mp 

391 

392 

393def _minpoly_sin(ex, x): 

394 """ 

395 Returns the minimal polynomial of ``sin(ex)`` 

396 see https://mathworld.wolfram.com/TrigonometryAngles.html 

397 """ 

398 c, a = ex.args[0].as_coeff_Mul() 

399 if a is pi: 

400 if c.is_rational: 

401 n = c.q 

402 q = sympify(n) 

403 if q.is_prime: 

404 # for a = pi*p/q with q odd prime, using chebyshevt 

405 # write sin(q*a) = mp(sin(a))*sin(a); 

406 # the roots of mp(x) are sin(pi*p/q) for p = 1,..., q - 1 

407 a = dup_chebyshevt(n, ZZ) 

408 return Add(*[x**(n - i - 1)*a[i] for i in range(n)]) 

409 if c.p == 1: 

410 if q == 9: 

411 return 64*x**6 - 96*x**4 + 36*x**2 - 3 

412 

413 if n % 2 == 1: 

414 # for a = pi*p/q with q odd, use 

415 # sin(q*a) = 0 to see that the minimal polynomial must be 

416 # a factor of dup_chebyshevt(n, ZZ) 

417 a = dup_chebyshevt(n, ZZ) 

418 a = [x**(n - i)*a[i] for i in range(n + 1)] 

419 r = Add(*a) 

420 _, factors = factor_list(r) 

421 res = _choose_factor(factors, x, ex) 

422 return res 

423 

424 expr = ((1 - cos(2*c*pi))/2)**S.Half 

425 res = _minpoly_compose(expr, x, QQ) 

426 return res 

427 

428 raise NotAlgebraic("%s does not seem to be an algebraic element" % ex) 

429 

430 

431def _minpoly_cos(ex, x): 

432 """ 

433 Returns the minimal polynomial of ``cos(ex)`` 

434 see https://mathworld.wolfram.com/TrigonometryAngles.html 

435 """ 

436 c, a = ex.args[0].as_coeff_Mul() 

437 if a is pi: 

438 if c.is_rational: 

439 if c.p == 1: 

440 if c.q == 7: 

441 return 8*x**3 - 4*x**2 - 4*x + 1 

442 if c.q == 9: 

443 return 8*x**3 - 6*x - 1 

444 elif c.p == 2: 

445 q = sympify(c.q) 

446 if q.is_prime: 

447 s = _minpoly_sin(ex, x) 

448 return _mexpand(s.subs({x:sqrt((1 - x)/2)})) 

449 

450 # for a = pi*p/q, cos(q*a) =T_q(cos(a)) = (-1)**p 

451 n = int(c.q) 

452 a = dup_chebyshevt(n, ZZ) 

453 a = [x**(n - i)*a[i] for i in range(n + 1)] 

454 r = Add(*a) - (-1)**c.p 

455 _, factors = factor_list(r) 

456 res = _choose_factor(factors, x, ex) 

457 return res 

458 

459 raise NotAlgebraic("%s does not seem to be an algebraic element" % ex) 

460 

461 

462def _minpoly_tan(ex, x): 

463 """ 

464 Returns the minimal polynomial of ``tan(ex)`` 

465 see https://github.com/sympy/sympy/issues/21430 

466 """ 

467 c, a = ex.args[0].as_coeff_Mul() 

468 if a is pi: 

469 if c.is_rational: 

470 c = c * 2 

471 n = int(c.q) 

472 a = n if c.p % 2 == 0 else 1 

473 terms = [] 

474 for k in range((c.p+1)%2, n+1, 2): 

475 terms.append(a*x**k) 

476 a = -(a*(n-k-1)*(n-k)) // ((k+1)*(k+2)) 

477 

478 r = Add(*terms) 

479 _, factors = factor_list(r) 

480 res = _choose_factor(factors, x, ex) 

481 return res 

482 

483 raise NotAlgebraic("%s does not seem to be an algebraic element" % ex) 

484 

485 

486def _minpoly_exp(ex, x): 

487 """ 

488 Returns the minimal polynomial of ``exp(ex)`` 

489 """ 

490 c, a = ex.args[0].as_coeff_Mul() 

491 if a == I*pi: 

492 if c.is_rational: 

493 q = sympify(c.q) 

494 if c.p == 1 or c.p == -1: 

495 if q == 3: 

496 return x**2 - x + 1 

497 if q == 4: 

498 return x**4 + 1 

499 if q == 6: 

500 return x**4 - x**2 + 1 

501 if q == 8: 

502 return x**8 + 1 

503 if q == 9: 

504 return x**6 - x**3 + 1 

505 if q == 10: 

506 return x**8 - x**6 + x**4 - x**2 + 1 

507 if q.is_prime: 

508 s = 0 

509 for i in range(q): 

510 s += (-x)**i 

511 return s 

512 

513 # x**(2*q) = product(factors) 

514 factors = [cyclotomic_poly(i, x) for i in divisors(2*q)] 

515 mp = _choose_factor(factors, x, ex) 

516 return mp 

517 else: 

518 raise NotAlgebraic("%s does not seem to be an algebraic element" % ex) 

519 raise NotAlgebraic("%s does not seem to be an algebraic element" % ex) 

520 

521 

522def _minpoly_rootof(ex, x): 

523 """ 

524 Returns the minimal polynomial of a ``CRootOf`` object. 

525 """ 

526 p = ex.expr 

527 p = p.subs({ex.poly.gens[0]:x}) 

528 _, factors = factor_list(p, x) 

529 result = _choose_factor(factors, x, ex) 

530 return result 

531 

532 

533def _minpoly_compose(ex, x, dom): 

534 """ 

535 Computes the minimal polynomial of an algebraic element 

536 using operations on minimal polynomials 

537 

538 Examples 

539 ======== 

540 

541 >>> from sympy import minimal_polynomial, sqrt, Rational 

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

543 >>> minimal_polynomial(sqrt(2) + 3*Rational(1, 3), x, compose=True) 

544 x**2 - 2*x - 1 

545 >>> minimal_polynomial(sqrt(y) + 1/y, x, compose=True) 

546 x**2*y**2 - 2*x*y - y**3 + 1 

547 

548 """ 

549 if ex.is_Rational: 

550 return ex.q*x - ex.p 

551 if ex is I: 

552 _, factors = factor_list(x**2 + 1, x, domain=dom) 

553 return x**2 + 1 if len(factors) == 1 else x - I 

554 

555 if ex is S.GoldenRatio: 

556 _, factors = factor_list(x**2 - x - 1, x, domain=dom) 

557 if len(factors) == 1: 

558 return x**2 - x - 1 

559 else: 

560 return _choose_factor(factors, x, (1 + sqrt(5))/2, dom=dom) 

561 

562 if ex is S.TribonacciConstant: 

563 _, factors = factor_list(x**3 - x**2 - x - 1, x, domain=dom) 

564 if len(factors) == 1: 

565 return x**3 - x**2 - x - 1 

566 else: 

567 fac = (1 + cbrt(19 - 3*sqrt(33)) + cbrt(19 + 3*sqrt(33))) / 3 

568 return _choose_factor(factors, x, fac, dom=dom) 

569 

570 if hasattr(dom, 'symbols') and ex in dom.symbols: 

571 return x - ex 

572 

573 if dom.is_QQ and _is_sum_surds(ex): 

574 # eliminate the square roots 

575 ex -= x 

576 while 1: 

577 ex1 = _separate_sq(ex) 

578 if ex1 is ex: 

579 return ex 

580 else: 

581 ex = ex1 

582 

583 if ex.is_Add: 

584 res = _minpoly_add(x, dom, *ex.args) 

585 elif ex.is_Mul: 

586 f = Factors(ex).factors 

587 r = sift(f.items(), lambda itx: itx[0].is_Rational and itx[1].is_Rational) 

588 if r[True] and dom == QQ: 

589 ex1 = Mul(*[bx**ex for bx, ex in r[False] + r[None]]) 

590 r1 = dict(r[True]) 

591 dens = [y.q for y in r1.values()] 

592 lcmdens = reduce(lcm, dens, 1) 

593 neg1 = S.NegativeOne 

594 expn1 = r1.pop(neg1, S.Zero) 

595 nums = [base**(y.p*lcmdens // y.q) for base, y in r1.items()] 

596 ex2 = Mul(*nums) 

597 mp1 = minimal_polynomial(ex1, x) 

598 # use the fact that in SymPy canonicalization products of integers 

599 # raised to rational powers are organized in relatively prime 

600 # bases, and that in ``base**(n/d)`` a perfect power is 

601 # simplified with the root 

602 # Powers of -1 have to be treated separately to preserve sign. 

603 mp2 = ex2.q*x**lcmdens - ex2.p*neg1**(expn1*lcmdens) 

604 ex2 = neg1**expn1 * ex2**Rational(1, lcmdens) 

605 res = _minpoly_op_algebraic_element(Mul, ex1, ex2, x, dom, mp1=mp1, mp2=mp2) 

606 else: 

607 res = _minpoly_mul(x, dom, *ex.args) 

608 elif ex.is_Pow: 

609 res = _minpoly_pow(ex.base, ex.exp, x, dom) 

610 elif ex.__class__ is sin: 

611 res = _minpoly_sin(ex, x) 

612 elif ex.__class__ is cos: 

613 res = _minpoly_cos(ex, x) 

614 elif ex.__class__ is tan: 

615 res = _minpoly_tan(ex, x) 

616 elif ex.__class__ is exp: 

617 res = _minpoly_exp(ex, x) 

618 elif ex.__class__ is CRootOf: 

619 res = _minpoly_rootof(ex, x) 

620 else: 

621 raise NotAlgebraic("%s does not seem to be an algebraic element" % ex) 

622 return res 

623 

624 

625@public 

626def minimal_polynomial(ex, x=None, compose=True, polys=False, domain=None): 

627 """ 

628 Computes the minimal polynomial of an algebraic element. 

629 

630 Parameters 

631 ========== 

632 

633 ex : Expr 

634 Element or expression whose minimal polynomial is to be calculated. 

635 

636 x : Symbol, optional 

637 Independent variable of the minimal polynomial 

638 

639 compose : boolean, optional (default=True) 

640 Method to use for computing minimal polynomial. If ``compose=True`` 

641 (default) then ``_minpoly_compose`` is used, if ``compose=False`` then 

642 groebner bases are used. 

643 

644 polys : boolean, optional (default=False) 

645 If ``True`` returns a ``Poly`` object else an ``Expr`` object. 

646 

647 domain : Domain, optional 

648 Ground domain 

649 

650 Notes 

651 ===== 

652 

653 By default ``compose=True``, the minimal polynomial of the subexpressions of ``ex`` 

654 are computed, then the arithmetic operations on them are performed using the resultant 

655 and factorization. 

656 If ``compose=False``, a bottom-up algorithm is used with ``groebner``. 

657 The default algorithm stalls less frequently. 

658 

659 If no ground domain is given, it will be generated automatically from the expression. 

660 

661 Examples 

662 ======== 

663 

664 >>> from sympy import minimal_polynomial, sqrt, solve, QQ 

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

666 

667 >>> minimal_polynomial(sqrt(2), x) 

668 x**2 - 2 

669 >>> minimal_polynomial(sqrt(2), x, domain=QQ.algebraic_field(sqrt(2))) 

670 x - sqrt(2) 

671 >>> minimal_polynomial(sqrt(2) + sqrt(3), x) 

672 x**4 - 10*x**2 + 1 

673 >>> minimal_polynomial(solve(x**3 + x + 3)[0], x) 

674 x**3 + x + 3 

675 >>> minimal_polynomial(sqrt(y), x) 

676 x**2 - y 

677 

678 """ 

679 

680 ex = sympify(ex) 

681 if ex.is_number: 

682 # not sure if it's always needed but try it for numbers (issue 8354) 

683 ex = _mexpand(ex, recursive=True) 

684 for expr in preorder_traversal(ex): 

685 if expr.is_AlgebraicNumber: 

686 compose = False 

687 break 

688 

689 if x is not None: 

690 x, cls = sympify(x), Poly 

691 else: 

692 x, cls = Dummy('x'), PurePoly 

693 

694 if not domain: 

695 if ex.free_symbols: 

696 domain = FractionField(QQ, list(ex.free_symbols)) 

697 else: 

698 domain = QQ 

699 if hasattr(domain, 'symbols') and x in domain.symbols: 

700 raise GeneratorsError("the variable %s is an element of the ground " 

701 "domain %s" % (x, domain)) 

702 

703 if compose: 

704 result = _minpoly_compose(ex, x, domain) 

705 result = result.primitive()[1] 

706 c = result.coeff(x**degree(result, x)) 

707 if c.is_negative: 

708 result = expand_mul(-result) 

709 return cls(result, x, field=True) if polys else result.collect(x) 

710 

711 if not domain.is_QQ: 

712 raise NotImplementedError("groebner method only works for QQ") 

713 

714 result = _minpoly_groebner(ex, x, cls) 

715 return cls(result, x, field=True) if polys else result.collect(x) 

716 

717 

718def _minpoly_groebner(ex, x, cls): 

719 """ 

720 Computes the minimal polynomial of an algebraic number 

721 using Groebner bases 

722 

723 Examples 

724 ======== 

725 

726 >>> from sympy import minimal_polynomial, sqrt, Rational 

727 >>> from sympy.abc import x 

728 >>> minimal_polynomial(sqrt(2) + 3*Rational(1, 3), x, compose=False) 

729 x**2 - 2*x - 1 

730 

731 """ 

732 

733 generator = numbered_symbols('a', cls=Dummy) 

734 mapping, symbols = {}, {} 

735 

736 def update_mapping(ex, exp, base=None): 

737 a = next(generator) 

738 symbols[ex] = a 

739 

740 if base is not None: 

741 mapping[ex] = a**exp + base 

742 else: 

743 mapping[ex] = exp.as_expr(a) 

744 

745 return a 

746 

747 def bottom_up_scan(ex): 

748 """ 

749 Transform a given algebraic expression *ex* into a multivariate 

750 polynomial, by introducing fresh variables with defining equations. 

751 

752 Explanation 

753 =========== 

754 

755 The critical elements of the algebraic expression *ex* are root 

756 extractions, instances of :py:class:`~.AlgebraicNumber`, and negative 

757 powers. 

758 

759 When we encounter a root extraction or an :py:class:`~.AlgebraicNumber` 

760 we replace this expression with a fresh variable ``a_i``, and record 

761 the defining polynomial for ``a_i``. For example, if ``a_0**(1/3)`` 

762 occurs, we will replace it with ``a_1``, and record the new defining 

763 polynomial ``a_1**3 - a_0``. 

764 

765 When we encounter a negative power we transform it into a positive 

766 power by algebraically inverting the base. This means computing the 

767 minimal polynomial in ``x`` for the base, inverting ``x`` modulo this 

768 poly (which generates a new polynomial) and then substituting the 

769 original base expression for ``x`` in this last polynomial. 

770 

771 We return the transformed expression, and we record the defining 

772 equations for new symbols using the ``update_mapping()`` function. 

773 

774 """ 

775 if ex.is_Atom: 

776 if ex is S.ImaginaryUnit: 

777 if ex not in mapping: 

778 return update_mapping(ex, 2, 1) 

779 else: 

780 return symbols[ex] 

781 elif ex.is_Rational: 

782 return ex 

783 elif ex.is_Add: 

784 return Add(*[ bottom_up_scan(g) for g in ex.args ]) 

785 elif ex.is_Mul: 

786 return Mul(*[ bottom_up_scan(g) for g in ex.args ]) 

787 elif ex.is_Pow: 

788 if ex.exp.is_Rational: 

789 if ex.exp < 0: 

790 minpoly_base = _minpoly_groebner(ex.base, x, cls) 

791 inverse = invert(x, minpoly_base).as_expr() 

792 base_inv = inverse.subs(x, ex.base).expand() 

793 

794 if ex.exp == -1: 

795 return bottom_up_scan(base_inv)