Coverage for /usr/lib/python3/dist-packages/sympy/polys/polyclasses.py: 27%

1"""OO layer for several polynomial representations. """

4from sympy.core.numbers import oo

5from sympy.core.sympify import CantSympify

6from sympy.polys.polyerrors import CoercionFailed, NotReversible, NotInvertible

7from sympy.polys.polyutils import PicklableWithSlots

10class GenericPoly(PicklableWithSlots):

11 """Base class for low-level polynomial representations. """

13 def ground_to_ring(f):

14 """Make the ground domain a ring. """

15 return f.set_domain(f.dom.get_ring())

17 def ground_to_field(f):

18 """Make the ground domain a field. """

19 return f.set_domain(f.dom.get_field())

21 def ground_to_exact(f):

22 """Make the ground domain exact. """

23 return f.set_domain(f.dom.get_exact())

25 @classmethod

26 def _perify_factors(per, result, include):

27 if include:

28 coeff, factors = result

30 factors = [ (per(g), k) for g, k in factors ]

32 if include:

33 return coeff, factors

34 else:

35 return factors

37from sympy.polys.densebasic import (

38 dmp_validate,

39 dup_normal, dmp_normal,

40 dup_convert, dmp_convert,

41 dmp_from_sympy,

42 dup_strip,

43 dup_degree, dmp_degree_in,

44 dmp_degree_list,

45 dmp_negative_p,

46 dup_LC, dmp_ground_LC,

47 dup_TC, dmp_ground_TC,

48 dmp_ground_nth,

49 dmp_one, dmp_ground,

50 dmp_zero_p, dmp_one_p, dmp_ground_p,

51 dup_from_dict, dmp_from_dict,

52 dmp_to_dict,

53 dmp_deflate,

54 dmp_inject, dmp_eject,

55 dmp_terms_gcd,

56 dmp_list_terms, dmp_exclude,

57 dmp_slice_in, dmp_permute,

58 dmp_to_tuple,)

60from sympy.polys.densearith import (

61 dmp_add_ground,

62 dmp_sub_ground,

63 dmp_mul_ground,

64 dmp_quo_ground,

65 dmp_exquo_ground,

66 dmp_abs,

67 dup_neg, dmp_neg,

68 dup_add, dmp_add,

69 dup_sub, dmp_sub,

70 dup_mul, dmp_mul,

71 dmp_sqr,

72 dup_pow, dmp_pow,

73 dmp_pdiv,

74 dmp_prem,

75 dmp_pquo,

76 dmp_pexquo,

77 dmp_div,

78 dup_rem, dmp_rem,

79 dmp_quo,

80 dmp_exquo,

81 dmp_add_mul, dmp_sub_mul,

82 dmp_max_norm,

83 dmp_l1_norm,

84 dmp_l2_norm_squared)

86from sympy.polys.densetools import (

87 dmp_clear_denoms,

88 dmp_integrate_in,

89 dmp_diff_in,

90 dmp_eval_in,

91 dup_revert,

92 dmp_ground_trunc,

93 dmp_ground_content,

94 dmp_ground_primitive,

95 dmp_ground_monic,

96 dmp_compose,

97 dup_decompose,

98 dup_shift,

99 dup_transform,

100 dmp_lift)

101

102from sympy.polys.euclidtools import (

103 dup_half_gcdex, dup_gcdex, dup_invert,

104 dmp_subresultants,

105 dmp_resultant,

106 dmp_discriminant,

107 dmp_inner_gcd,

108 dmp_gcd,

109 dmp_lcm,

110 dmp_cancel)

111

112from sympy.polys.sqfreetools import (

113 dup_gff_list,

114 dmp_norm,

115 dmp_sqf_p,

116 dmp_sqf_norm,

117 dmp_sqf_part,

118 dmp_sqf_list, dmp_sqf_list_include)

119

120from sympy.polys.factortools import (

121 dup_cyclotomic_p, dmp_irreducible_p,

122 dmp_factor_list, dmp_factor_list_include)

123

124from sympy.polys.rootisolation import (

125 dup_isolate_real_roots_sqf,

126 dup_isolate_real_roots,

127 dup_isolate_all_roots_sqf,

128 dup_isolate_all_roots,

129 dup_refine_real_root,

130 dup_count_real_roots,

131 dup_count_complex_roots,

132 dup_sturm,

133 dup_cauchy_upper_bound,

134 dup_cauchy_lower_bound,

135 dup_mignotte_sep_bound_squared)

136

137from sympy.polys.polyerrors import (

138 UnificationFailed,

139 PolynomialError)

140

141

142def init_normal_DMP(rep, lev, dom):

143 return DMP(dmp_normal(rep, lev, dom), dom, lev)

144

145

146class DMP(PicklableWithSlots, CantSympify):

147 """Dense Multivariate Polynomials over `K`. """

148

149 __slots__ = ('rep', 'lev', 'dom', 'ring')

150

151 def __init__(self, rep, dom, lev=None, ring=None):

152 if lev is not None:

153 # Not possible to check with isinstance

154 if type(rep) is dict:

155 rep = dmp_from_dict(rep, lev, dom)

156 elif not isinstance(rep, list):

157 rep = dmp_ground(dom.convert(rep), lev)

158 else:

159 rep, lev = dmp_validate(rep)

160

161 self.rep = rep

162 self.lev = lev

163 self.dom = dom

164 self.ring = ring

165

166 def __repr__(f):

167 return "%s(%s, %s, %s)" % (f.__class__.__name__, f.rep, f.dom, f.ring)

168

169 def __hash__(f):

170 return hash((f.__class__.__name__, f.to_tuple(), f.lev, f.dom, f.ring))

171

172 def unify(f, g):

173 """Unify representations of two multivariate polynomials. """

174 if not isinstance(g, DMP) or f.lev != g.lev:

175 raise UnificationFailed("Cannot unify %s with %s" % (f, g))

176

177 if f.dom == g.dom and f.ring == g.ring:

178 return f.lev, f.dom, f.per, f.rep, g.rep

179 else:

180 lev, dom = f.lev, f.dom.unify(g.dom)

181 ring = f.ring

182 if g.ring is not None:

183 if ring is not None:

184 ring = ring.unify(g.ring)

185 else:

186 ring = g.ring

187

188 F = dmp_convert(f.rep, lev, f.dom, dom)

189 G = dmp_convert(g.rep, lev, g.dom, dom)

190

191 def per(rep, dom=dom, lev=lev, kill=False):

192 if kill:

193 if not lev:

194 return rep

195 else:

196 lev -= 1

197

198 return DMP(rep, dom, lev, ring)

199

200 return lev, dom, per, F, G

201

202 def per(f, rep, dom=None, kill=False, ring=None):

203 """Create a DMP out of the given representation. """

204 lev = f.lev

205

206 if kill:

207 if not lev:

208 return rep

209 else:

210 lev -= 1

211

212 if dom is None:

213 dom = f.dom

214

215 if ring is None:

216 ring = f.ring

217

218 return DMP(rep, dom, lev, ring)

219

220 @classmethod

221 def zero(cls, lev, dom, ring=None):

222 return DMP(0, dom, lev, ring)

223

224 @classmethod

225 def one(cls, lev, dom, ring=None):

226 return DMP(1, dom, lev, ring)

227

228 @classmethod

229 def from_list(cls, rep, lev, dom):

230 """Create an instance of ``cls`` given a list of native coefficients. """

231 return cls(dmp_convert(rep, lev, None, dom), dom, lev)

232

233 @classmethod

234 def from_sympy_list(cls, rep, lev, dom):

235 """Create an instance of ``cls`` given a list of SymPy coefficients. """

236 return cls(dmp_from_sympy(rep, lev, dom), dom, lev)

237

238 def to_dict(f, zero=False):

239 """Convert ``f`` to a dict representation with native coefficients. """

240 return dmp_to_dict(f.rep, f.lev, f.dom, zero=zero)

241

242 def to_sympy_dict(f, zero=False):

243 """Convert ``f`` to a dict representation with SymPy coefficients. """

244 rep = dmp_to_dict(f.rep, f.lev, f.dom, zero=zero)

245

246 for k, v in rep.items():

247 rep[k] = f.dom.to_sympy(v)

248

249 return rep

250

251 def to_list(f):

252 """Convert ``f`` to a list representation with native coefficients. """

253 return f.rep

254

255 def to_sympy_list(f):

256 """Convert ``f`` to a list representation with SymPy coefficients. """

257 def sympify_nested_list(rep):

258 out = []

259 for val in rep:

260 if isinstance(val, list):

261 out.append(sympify_nested_list(val))

262 else:

263 out.append(f.dom.to_sympy(val))

264 return out

265

266 return sympify_nested_list(f.rep)

267

268 def to_tuple(f):

269 """

270 Convert ``f`` to a tuple representation with native coefficients.

271

272 This is needed for hashing.

273 """

274 return dmp_to_tuple(f.rep, f.lev)

275

276 @classmethod

277 def from_dict(cls, rep, lev, dom):

278 """Construct and instance of ``cls`` from a ``dict`` representation. """

279 return cls(dmp_from_dict(rep, lev, dom), dom, lev)

280

281 @classmethod

282 def from_monoms_coeffs(cls, monoms, coeffs, lev, dom, ring=None):

283 return DMP(dict(list(zip(monoms, coeffs))), dom, lev, ring)

284

285 def to_ring(f):

286 """Make the ground domain a ring. """

287 return f.convert(f.dom.get_ring())

288

289 def to_field(f):

290 """Make the ground domain a field. """

291 return f.convert(f.dom.get_field())

292

293 def to_exact(f):

294 """Make the ground domain exact. """

295 return f.convert(f.dom.get_exact())

296

297 def convert(f, dom):

298 """Convert the ground domain of ``f``. """

299 if f.dom == dom:

300 return f

301 else:

302 return DMP(dmp_convert(f.rep, f.lev, f.dom, dom), dom, f.lev)

303

304 def slice(f, m, n, j=0):

305 """Take a continuous subsequence of terms of ``f``. """

306 return f.per(dmp_slice_in(f.rep, m, n, j, f.lev, f.dom))

307

308 def coeffs(f, order=None):

309 """Returns all non-zero coefficients from ``f`` in lex order. """

310 return [ c for _, c in dmp_list_terms(f.rep, f.lev, f.dom, order=order) ]

311

312 def monoms(f, order=None):

313 """Returns all non-zero monomials from ``f`` in lex order. """

314 return [ m for m, _ in dmp_list_terms(f.rep, f.lev, f.dom, order=order) ]

315

316 def terms(f, order=None):

317 """Returns all non-zero terms from ``f`` in lex order. """

318 return dmp_list_terms(f.rep, f.lev, f.dom, order=order)

319

320 def all_coeffs(f):

321 """Returns all coefficients from ``f``. """

322 if not f.lev:

323 if not f:

324 return [f.dom.zero]

325 else:

326 return list(f.rep)

327 else:

328 raise PolynomialError('multivariate polynomials not supported')

329

330 def all_monoms(f):

331 """Returns all monomials from ``f``. """

332 if not f.lev:

333 n = dup_degree(f.rep)

334

335 if n < 0:

336 return [(0,)]

337 else:

338 return [ (n - i,) for i, c in enumerate(f.rep) ]

339 else:

340 raise PolynomialError('multivariate polynomials not supported')

341

342 def all_terms(f):

343 """Returns all terms from a ``f``. """

344 if not f.lev:

345 n = dup_degree(f.rep)

346

347 if n < 0:

348 return [((0,), f.dom.zero)]

349 else:

350 return [ ((n - i,), c) for i, c in enumerate(f.rep) ]

351 else:

352 raise PolynomialError('multivariate polynomials not supported')

353

354 def lift(f):

355 """Convert algebraic coefficients to rationals. """

356 return f.per(dmp_lift(f.rep, f.lev, f.dom), dom=f.dom.dom)

357

358 def deflate(f):

359 """Reduce degree of `f` by mapping `x_i^m` to `y_i`. """

360 J, F = dmp_deflate(f.rep, f.lev, f.dom)

361 return J, f.per(F)

362

363 def inject(f, front=False):

364 """Inject ground domain generators into ``f``. """

365 F, lev = dmp_inject(f.rep, f.lev, f.dom, front=front)

366 return f.__class__(F, f.dom.dom, lev)

367

368 def eject(f, dom, front=False):

369 """Eject selected generators into the ground domain. """

370 F = dmp_eject(f.rep, f.lev, dom, front=front)

371 return f.__class__(F, dom, f.lev - len(dom.symbols))

372

373 def exclude(f):

374 r"""

375 Remove useless generators from ``f``.

376

377 Returns the removed generators and the new excluded ``f``.

378

379 Examples

380 ========

381

382 >>> from sympy.polys.polyclasses import DMP

383 >>> from sympy.polys.domains import ZZ

384

385 >>> DMP([[[ZZ(1)]], [[ZZ(1)], [ZZ(2)]]], ZZ).exclude()

386 ([2], DMP([[1], [1, 2]], ZZ, None))

387

388 """

389 J, F, u = dmp_exclude(f.rep, f.lev, f.dom)

390 return J, f.__class__(F, f.dom, u)

391

392 def permute(f, P):

393 r"""

394 Returns a polynomial in `K[x_{P(1)}, ..., x_{P(n)}]`.

395

396 Examples

397 ========

398

399 >>> from sympy.polys.polyclasses import DMP

400 >>> from sympy.polys.domains import ZZ

401

402 >>> DMP([[[ZZ(2)], [ZZ(1), ZZ(0)]], [[]]], ZZ).permute([1, 0, 2])

403 DMP([[[2], []], [[1, 0], []]], ZZ, None)

404

405 >>> DMP([[[ZZ(2)], [ZZ(1), ZZ(0)]], [[]]], ZZ).permute([1, 2, 0])

406 DMP([[[1], []], [[2, 0], []]], ZZ, None)

407

408 """

409 return f.per(dmp_permute(f.rep, P, f.lev, f.dom))

410

411 def terms_gcd(f):

412 """Remove GCD of terms from the polynomial ``f``. """

413 J, F = dmp_terms_gcd(f.rep, f.lev, f.dom)

414 return J, f.per(F)

415

416 def add_ground(f, c):

417 """Add an element of the ground domain to ``f``. """

418 return f.per(dmp_add_ground(f.rep, f.dom.convert(c), f.lev, f.dom))

419

420 def sub_ground(f, c):

421 """Subtract an element of the ground domain from ``f``. """

422 return f.per(dmp_sub_ground(f.rep, f.dom.convert(c), f.lev, f.dom))

423

424 def mul_ground(f, c):

425 """Multiply ``f`` by a an element of the ground domain. """

426 return f.per(dmp_mul_ground(f.rep, f.dom.convert(c), f.lev, f.dom))

427

428 def quo_ground(f, c):

429 """Quotient of ``f`` by a an element of the ground domain. """

430 return f.per(dmp_quo_ground(f.rep, f.dom.convert(c), f.lev, f.dom))

431

432 def exquo_ground(f, c):

433 """Exact quotient of ``f`` by a an element of the ground domain. """

434 return f.per(dmp_exquo_ground(f.rep, f.dom.convert(c), f.lev, f.dom))

435

436 def abs(f):

437 """Make all coefficients in ``f`` positive. """

438 return f.per(dmp_abs(f.rep, f.lev, f.dom))

439

440 def neg(f):

441 """Negate all coefficients in ``f``. """

442 return f.per(dmp_neg(f.rep, f.lev, f.dom))

443

444 def add(f, g):

445 """Add two multivariate polynomials ``f`` and ``g``. """

446 lev, dom, per, F, G = f.unify(g)

447 return per(dmp_add(F, G, lev, dom))

448

449 def sub(f, g):

450 """Subtract two multivariate polynomials ``f`` and ``g``. """

451 lev, dom, per, F, G = f.unify(g)

452 return per(dmp_sub(F, G, lev, dom))

453

454 def mul(f, g):

455 """Multiply two multivariate polynomials ``f`` and ``g``. """

456 lev, dom, per, F, G = f.unify(g)

457 return per(dmp_mul(F, G, lev, dom))

458

459 def sqr(f):

460 """Square a multivariate polynomial ``f``. """

461 return f.per(dmp_sqr(f.rep, f.lev, f.dom))

462

463 def pow(f, n):

464 """Raise ``f`` to a non-negative power ``n``. """

465 if isinstance(n, int):

466 return f.per(dmp_pow(f.rep, n, f.lev, f.dom))

467 else:

468 raise TypeError("``int`` expected, got %s" % type(n))

469

470 def pdiv(f, g):

471 """Polynomial pseudo-division of ``f`` and ``g``. """

472 lev, dom, per, F, G = f.unify(g)

473 q, r = dmp_pdiv(F, G, lev, dom)

474 return per(q), per(r)

475

476 def prem(f, g):

477 """Polynomial pseudo-remainder of ``f`` and ``g``. """

478 lev, dom, per, F, G = f.unify(g)

479 return per(dmp_prem(F, G, lev, dom))

480

481 def pquo(f, g):

482 """Polynomial pseudo-quotient of ``f`` and ``g``. """

483 lev, dom, per, F, G = f.unify(g)

484 return per(dmp_pquo(F, G, lev, dom))

485

486 def pexquo(f, g):

487 """Polynomial exact pseudo-quotient of ``f`` and ``g``. """

488 lev, dom, per, F, G = f.unify(g)

489 return per(dmp_pexquo(F, G, lev, dom))

490

491 def div(f, g):

492 """Polynomial division with remainder of ``f`` and ``g``. """

493 lev, dom, per, F, G = f.unify(g)

494 q, r = dmp_div(F, G, lev, dom)

495 return per(q), per(r)

496

497 def rem(f, g):

498 """Computes polynomial remainder of ``f`` and ``g``. """

499 lev, dom, per, F, G = f.unify(g)

500 return per(dmp_rem(F, G, lev, dom))

501

502 def quo(f, g):

503 """Computes polynomial quotient of ``f`` and ``g``. """

504 lev, dom, per, F, G = f.unify(g)

505 return per(dmp_quo(F, G, lev, dom))

506

507 def exquo(f, g):

508 """Computes polynomial exact quotient of ``f`` and ``g``. """

509 lev, dom, per, F, G = f.unify(g)

510 res = per(dmp_exquo(F, G, lev, dom))

511 if f.ring and res not in f.ring:

512 from sympy.polys.polyerrors import ExactQuotientFailed

513 raise ExactQuotientFailed(f, g, f.ring)

514 return res

515

516 def degree(f, j=0):

517 """Returns the leading degree of ``f`` in ``x_j``. """

518 if isinstance(j, int):

519 return dmp_degree_in(f.rep, j, f.lev)

520 else:

521 raise TypeError("``int`` expected, got %s" % type(j))

522

523 def degree_list(f):

524 """Returns a list of degrees of ``f``. """

525 return dmp_degree_list(f.rep, f.lev)

526

527 def total_degree(f):

528 """Returns the total degree of ``f``. """

529 return max(sum(m) for m in f.monoms())

530

531 def homogenize(f, s):

532 """Return homogeneous polynomial of ``f``"""

533 td = f.total_degree()

534 result = {}

535 new_symbol = (s == len(f.terms()[0][0]))

536 for term in f.terms():

537 d = sum(term[0])

538 if d < td:

539 i = td - d

540 else:

541 i = 0

542 if new_symbol:

543 result[term[0] + (i,)] = term[1]

544 else:

545 l = list(term[0])

546 l[s] += i

547 result[tuple(l)] = term[1]

548 return DMP(result, f.dom, f.lev + int(new_symbol), f.ring)

549

550 def homogeneous_order(f):

551 """Returns the homogeneous order of ``f``. """

552 if f.is_zero:

553 return -oo

554

555 monoms = f.monoms()

556 tdeg = sum(monoms[0])

557

558 for monom in monoms:

559 _tdeg = sum(monom)

560

561 if _tdeg != tdeg:

562 return None

563

564 return tdeg

565

566 def LC(f):

567 """Returns the leading coefficient of ``f``. """

568 return dmp_ground_LC(f.rep, f.lev, f.dom)

569

570 def TC(f):

571 """Returns the trailing coefficient of ``f``. """

572 return dmp_ground_TC(f.rep, f.lev, f.dom)

573

574 def nth(f, *N):

575 """Returns the ``n``-th coefficient of ``f``. """

576 if all(isinstance(n, int) for n in N):

577 return dmp_ground_nth(f.rep, N, f.lev, f.dom)

578 else:

579 raise TypeError("a sequence of integers expected")

580

581 def max_norm(f):

582 """Returns maximum norm of ``f``. """

583 return dmp_max_norm(f.rep, f.lev, f.dom)

584

585 def l1_norm(f):

586 """Returns l1 norm of ``f``. """

587 return dmp_l1_norm(f.rep, f.lev, f.dom)

588

589 def l2_norm_squared(f):

590 """Return squared l2 norm of ``f``. """

591 return dmp_l2_norm_squared(f.rep, f.lev, f.dom)

592

593 def clear_denoms(f):

594 """Clear denominators, but keep the ground domain. """

595 coeff, F = dmp_clear_denoms(f.rep, f.lev, f.dom)

596 return coeff, f.per(F)

597

598 def integrate(f, m=1, j=0):

599 """Computes the ``m``-th order indefinite integral of ``f`` in ``x_j``. """

600 if not isinstance(m, int):

601 raise TypeError("``int`` expected, got %s" % type(m))

602

603 if not isinstance(j, int):

604 raise TypeError("``int`` expected, got %s" % type(j))

605

606 return f.per(dmp_integrate_in(f.rep, m, j, f.lev, f.dom))

607

608 def diff(f, m=1, j=0):

609 """Computes the ``m``-th order derivative of ``f`` in ``x_j``. """

610 if not isinstance(m, int):

611 raise TypeError("``int`` expected, got %s" % type(m))

612

613 if not isinstance(j, int):

614 raise TypeError("``int`` expected, got %s" % type(j))

615

616 return f.per(dmp_diff_in(f.rep, m, j, f.lev, f.dom))

617

618 def eval(f, a, j=0):

619 """Evaluates ``f`` at the given point ``a`` in ``x_j``. """

620 if not isinstance(j, int):

621 raise TypeError("``int`` expected, got %s" % type(j))

622

623 return f.per(dmp_eval_in(f.rep,

624 f.dom.convert(a), j, f.lev, f.dom), kill=True)

625

626 def half_gcdex(f, g):

627 """Half extended Euclidean algorithm, if univariate. """

628 lev, dom, per, F, G = f.unify(g)

629

630 if not lev:

631 s, h = dup_half_gcdex(F, G, dom)

632 return per(s), per(h)

633 else:

634 raise ValueError('univariate polynomial expected')

635

636 def gcdex(f, g):

637 """Extended Euclidean algorithm, if univariate. """

638 lev, dom, per, F, G = f.unify(g)

639

640 if not lev:

641 s, t, h = dup_gcdex(F, G, dom)

642 return per(s), per(t), per(h)

643 else:

644 raise ValueError('univariate polynomial expected')

645

646 def invert(f, g):

647 """Invert ``f`` modulo ``g``, if possible. """

648 lev, dom, per, F, G = f.unify(g)

649

650 if not lev:

651 return per(dup_invert(F, G, dom))

652 else:

653 raise ValueError('univariate polynomial expected')

654

655 def revert(f, n):

656 """Compute ``f**(-1)`` mod ``x**n``. """

657 if not f.lev:

658 return f.per(dup_revert(f.rep, n, f.dom))

659 else:

660 raise ValueError('univariate polynomial expected')

661

662 def subresultants(f, g):

663 """Computes subresultant PRS sequence of ``f`` and ``g``. """

664 lev, dom, per, F, G = f.unify(g)

665 R = dmp_subresultants(F, G, lev, dom)

666 return list(map(per, R))

667

668 def resultant(f, g, includePRS=False):

669 """Computes resultant of ``f`` and ``g`` via PRS. """

670 lev, dom, per, F, G = f.unify(g)

671 if includePRS:

672 res, R = dmp_resultant(F, G, lev, dom, includePRS=includePRS)

673 return per(res, kill=True), list(map(per, R))

674 return per(dmp_resultant(F, G, lev, dom), kill=True)

675

676 def discriminant(f):

677 """Computes discriminant of ``f``. """

678 return f.per(dmp_discriminant(f.rep, f.lev, f.dom), kill=True)

679

680 def cofactors(f, g):

681 """Returns GCD of ``f`` and ``g`` and their cofactors. """

682 lev, dom, per, F, G = f.unify(g)

683 h, cff, cfg = dmp_inner_gcd(F, G, lev, dom)

684 return per(h), per(cff), per(cfg)

685

686 def gcd(f, g):

687 """Returns polynomial GCD of ``f`` and ``g``. """

688 lev, dom, per, F, G = f.unify(g)

689 return per(dmp_gcd(F, G, lev, dom))

690

691 def lcm(f, g):

692 """Returns polynomial LCM of ``f`` and ``g``. """

693 lev, dom, per, F, G = f.unify(g)

694 return per(dmp_lcm(F, G, lev, dom))

695

696 def cancel(f, g, include=True):

697 """Cancel common factors in a rational function ``f/g``. """

698 lev, dom, per, F, G = f.unify(g)

699

700 if include:

701 F, G = dmp_cancel(F, G, lev, dom, include=True)

702 else:

703 cF, cG, F, G = dmp_cancel(F, G, lev, dom, include=False)

704

705 F, G = per(F), per(G)

706

707 if include:

708 return F, G

709 else:

710 return cF, cG, F, G

711

712 def trunc(f, p):

713 """Reduce ``f`` modulo a constant ``p``. """

714 return f.per(dmp_ground_trunc(f.rep, f.dom.convert(p), f.lev, f.dom))

715

716 def monic(f):

717 """Divides all coefficients by ``LC(f)``. """

718 return f.per(dmp_ground_monic(f.rep, f.lev, f.dom))

719

720 def content(f):

721 """Returns GCD of polynomial coefficients. """

722 return dmp_ground_content(f.rep, f.lev, f.dom)

723

724 def primitive(f):

725 """Returns content and a primitive form of ``f``. """

726 cont, F = dmp_ground_primitive(f.rep, f.lev, f.dom)

727 return cont, f.per(F)

728

729 def compose(f, g):

730 """Computes functional composition of ``f`` and ``g``. """

731 lev, dom, per, F, G = f.unify(g)

732 return per(dmp_compose(F, G, lev, dom))

733

734 def decompose(f):

735 """Computes functional decomposition of ``f``. """

736 if not f.lev:

737 return list(map(f.per, dup_decompose(f.rep, f.dom)))

738 else:

739 raise ValueError('univariate polynomial expected')

740

741 def shift(f, a):

742 """Efficiently compute Taylor shift ``f(x + a)``. """

743 if not f.lev:

744 return f.per(dup_shift(f.rep, f.dom.convert(a), f.dom))

745 else:

746 raise ValueError('univariate polynomial expected')

747

748 def transform(f, p, q):

749 """Evaluate functional transformation ``q**n * f(p/q)``."""

750 if f.lev:

751 raise ValueError('univariate polynomial expected')

752

753 lev, dom, per, P, Q = p.unify(q)

754 lev, dom, per, F, P = f.unify(per(P, dom, lev))

755 lev, dom, per, F, Q = per(F, dom, lev).unify(per(Q, dom, lev))

756

757 if not lev:

758 return per(dup_transform(F, P, Q, dom))

759 else:

760 raise ValueError('univariate polynomial expected')

761

762 def sturm(f):

763 """Computes the Sturm sequence of ``f``. """

764 if not f.lev:

765 return list(map(f.per, dup_sturm(f.rep, f.dom)))

766 else:

767 raise ValueError('univariate polynomial expected')

768

769 def cauchy_upper_bound(f):

770 """Computes the Cauchy upper bound on the roots of ``f``. """

771 if not f.lev:

772 return dup_cauchy_upper_bound(f.rep, f.dom)

773 else:

774 raise ValueError('univariate polynomial expected')

775

776 def cauchy_lower_bound(f):

777 """Computes the Cauchy lower bound on the nonzero roots of ``f``. """

778 if not f.lev:

779 return dup_cauchy_lower_bound(f.rep, f.dom)

780 else:

781 raise ValueError('univariate polynomial expected')

782

783 def mignotte_sep_bound_squared(f):

784 """Computes the squared Mignotte bound on root separations of ``f``. """

785 if not f.lev:

786 return dup_mignotte_sep_bound_squared(f.rep, f.dom)

787 else:

788 raise ValueError('univariate polynomial expected')

789

790 def gff_list(f):

791 """Computes greatest factorial factorization of ``f``. """

792 if not f.lev:

793 return [ (f.per(g), k) for g, k in dup_gff_list(f.rep, f.dom) ]

794 else:

795 raise ValueError('univariate polynomial expected')

796

797 def norm(f):

798 """Computes ``Norm(f)``."""

799 r = dmp_norm(f.rep, f.lev, f.dom)

800 return f.per(r, dom=f.dom.dom)

801

802 def sqf_norm(f):

803 """Computes square-free norm of ``f``. """

804 s, g, r = dmp_sqf_norm(f.rep, f.lev, f.dom)

805 return s, f.per(g), f.per(r, dom=f.dom.dom)

806

807 def sqf_part(f):

808 """Computes square-free part of ``f``. """

809 return f.per(dmp_sqf_part(f.rep, f.lev, f.dom))

810

811 def sqf_list(f, all=False):

812 """Returns a list of square-free factors of ``f``. """

813 coeff, factors = dmp_sqf_list(f.rep, f.lev, f.dom, all)

814 return coeff, [ (f.per(g), k) for g, k in factors ]

815

816 def sqf_list_include(f, all=False):

817 """Returns a list of square-free factors of ``f``. """

818 factors = dmp_sqf_list_include(f.rep, f.lev, f.dom, all)

819 return [ (f.per(g), k) for g, k in factors ]

820

821 def factor_list(f):

822 """Returns a list of irreducible factors of ``f``. """

823 coeff, factors = dmp_factor_list(f.rep, f.lev, f.dom)

824 return coeff, [ (f.per(g), k) for g, k in factors ]

825

826 def factor_list_include(f):

827 """Returns a list of irreducible factors of ``f``. """

828 factors = dmp_factor_list_include(f.rep, f.lev, f.dom)

829 return [ (f.per(g), k) for g, k in factors ]

830

831 def intervals(f, all=False, eps=None, inf=None, sup=None, fast=False, sqf=False):

832 """Compute isolating intervals for roots of ``f``. """

833 if not f.lev:

834 if not all:

835 if not sqf:

836 return dup_isolate_real_roots(f.rep, f.dom, eps=eps, inf=inf, sup=sup, fast=fast)

837 else:

838 return dup_isolate_real_roots_sqf(f.rep, f.dom, eps=eps, inf=inf, sup=sup, fast=fast)

839 else:

840 if not sqf:

841 return dup_isolate_all_roots(f.rep, f.dom, eps=eps, inf=inf, sup=sup, fast=fast)

842 else:

843 return dup_isolate_all_roots_sqf(f.rep, f.dom, eps=eps, inf=inf, sup=sup, fast=fast)

844 else:

845 raise PolynomialError(

846 "Cannot isolate roots of a multivariate polynomial")

847

848 def refine_root(f, s, t, eps=None, steps=None, fast=False):

849 """

850 Refine an isolating interval to the given precision.

851

852 ``eps`` should be a rational number.

853

854 """

855 if not f.lev:

856 return dup_refine_real_root(f.rep, s, t, f.dom, eps=eps, steps=steps, fast=fast)

857 else:

858 raise PolynomialError(

859 "Cannot refine a root of a multivariate polynomial")

860

861 def count_real_roots(f, inf=None, sup=None):

862 """Return the number of real roots of ``f`` in ``[inf, sup]``. """

863 return dup_count_real_roots(f.rep, f.dom, inf=inf, sup=sup)

864

865 def count_complex_roots(f, inf=None, sup=None):

866 """Return the number of complex roots of ``f`` in ``[inf, sup]``. """

867 return dup_count_complex_roots(f.rep, f.dom, inf=inf, sup=sup)

868

869 @property

870 def is_zero(f):

871 """Returns ``True`` if ``f`` is a zero polynomial. """

872 return dmp_zero_p(f.rep, f.lev)

873

874 @property

875 def is_one(f):

876 """Returns ``True`` if ``f`` is a unit polynomial. """

877 return dmp_one_p(f.rep, f.lev, f.dom)

878

879 @property

880 def is_ground(f):

881 """Returns ``True`` if ``f`` is an element of the ground domain. """

882 return dmp_ground_p(f.rep, None, f.lev)

883

884 @property

885 def is_sqf(f):

886 """Returns ``True`` if ``f`` is a square-free polynomial. """

887 return dmp_sqf_p(f.rep, f.lev, f.dom)

888

889 @property

890 def is_monic(f):

891 """Returns ``True`` if the leading coefficient of ``f`` is one. """

892 return f.dom.is_one(dmp_ground_LC(f.rep, f.lev, f.dom))

893

894 @property

895 def is_primitive(f):

896 """Returns ``True`` if the GCD of the coefficients of ``f`` is one. """

897 return f.dom.is_one(dmp_ground_content(f.rep, f.lev, f.dom))

898

899 @property

900 def is_linear(f):

901 """Returns ``True`` if ``f`` is linear in all its variables. """

902 return all(sum(monom) <= 1 for monom in dmp_to_dict(f.rep, f.lev, f.dom).keys())

903

904 @property

905 def is_quadratic(f):

906 """Returns ``True`` if ``f`` is quadratic in all its variables. """

907 return all(sum(monom) <= 2 for monom in dmp_to_dict(f.rep, f.lev, f.dom).keys())

908

909 @property

910 def is_monomial(f):

911 """Returns ``True`` if ``f`` is zero or has only one term. """

912 return len(f.to_dict()) <= 1

913

914 @property

915 def is_homogeneous(f):

916 """Returns ``True`` if ``f`` is a homogeneous polynomial. """

917 return f.homogeneous_order() is not None

918

919 @property

920 def is_irreducible(f):

921 """Returns ``True`` if ``f`` has no factors over its domain. """

922 return dmp_irreducible_p(f.rep, f.lev, f.dom)

923

924 @property

925 def is_cyclotomic(f):

926 """Returns ``True`` if ``f`` is a cyclotomic polynomial. """

927 if not f.lev:

928 return dup_cyclotomic_p(f.rep, f.dom)

929 else:

930 return False

931

932 def __abs__(f):

933 return f.abs()

934

935 def __neg__(f):

936 return f.neg()

937

938 def __add__(f, g):

939 if not isinstance(g, DMP):

940 try:

941 g = f.per(dmp_ground(f.dom.convert(g), f.lev))

942 except TypeError:

943 return NotImplemented

944 except (CoercionFailed, NotImplementedError):

945 if f.ring is not None:

946 try:

947 g = f.ring.convert(g)

948 except (CoercionFailed, NotImplementedError):

949 return NotImplemented

950

951 return f.add(g)

952

953 def __radd__(f, g):

954 return f.__add__(g)

955

956 def __sub__(f, g):

957 if not isinstance(g, DMP):

958 try:

959 g = f.per(dmp_ground(f.dom.convert(g), f.lev))

960 except TypeError:

961 return NotImplemented

962 except (CoercionFailed, NotImplementedError):

963 if f.ring is not None:

964 try:

965 g = f.ring.convert(g)

966 except (CoercionFailed, NotImplementedError):

967 return NotImplemented

968

969 return f.sub(g)

970

971 def __rsub__(f, g):

972 return (-f).__add__(g)

973

974 def __mul__(f, g):

975 if isinstance(g, DMP):

976 return f.mul(g)

977 else:

978 try:

979 return f.mul_ground(g)

980 except TypeError:

981 return NotImplemented

982 except (CoercionFailed, NotImplementedError):

983 if f.ring is not None:

984 try:

985 return f.mul(f.ring.convert(g))

986 except (CoercionFailed, NotImplementedError):

987 pass

988 return NotImplemented

989

990 def __truediv__(f, g):

991 if isinstance(g, DMP):

992 return f.exquo(g)

993 else:

994 try:

995 return f.mul_ground(g)

996 except TypeError:

997 return NotImplemented

998 except (CoercionFailed, NotImplementedError):

999 if f.ring is not None:

1000 try:

1001 return f.exquo(f.ring.convert(g))

1002 except (CoercionFailed, NotImplementedError):

1003 pass

1004 return NotImplemented

1005

1006 def __rtruediv__(f, g):

1007 if isinstance(g, DMP):

1008 return g.exquo(f)

1009 elif f.ring is not None:

1010 try:

1011 return f.ring.convert(g).exquo(f)

1012 except (CoercionFailed, NotImplementedError):

1013 pass

1014 return NotImplemented

1015

1016 def __rmul__(f, g):

1017 return f.__mul__(g)

1018

1019 def __pow__(f, n):

1020 return f.pow(n)

1021

1022 def __divmod__(f, g):

1023 return f.div(g)

1024

1025 def __mod__(f, g):

1026 return f.rem(g)

1027

1028 def __floordiv__(f, g):

1029 if isinstance(g, DMP):

1030 return f.quo(g)

1031 else:

1032 try:

1033 return f.quo_ground(g)

1034 except TypeError:

1035 return NotImplemented

1036

1037 def __eq__(f, g):

1038 try:

1039 _, _, _, F, G = f.unify(g)

1040

1041 if f.lev == g.lev:

1042 return F == G

1043 except UnificationFailed:

1044 pass

1045

1046 return False

1047

1048 def __ne__(f, g):

1049 return not f == g

1050

1051 def eq(f, g, strict=False):

1052 if not strict:

1053 return f == g

1054 else:

1055 return f._strict_eq(g)

1056

1057 def ne(f, g, strict=False):

1058 return not f.eq(g, strict=strict)

1059

1060 def _strict_eq(f, g):

1061 return isinstance(g, f.__class__) and f.lev == g.lev \

1062 and f.dom == g.dom \

1063 and f.rep == g.rep

1064

1065 def __lt__(f, g):

1066 _, _, _, F, G = f.unify(g)

1067 return F < G

1068

1069 def __le__(f, g):

1070 _, _, _, F, G = f.unify(g)

1071 return F <= G

1072

1073 def __gt__(f, g):

1074 _, _, _, F, G = f.unify(g)

1075 return F > G

1076

1077 def __ge__(f, g):

1078 _, _, _, F, G = f.unify(g)

1079 return F >= G

1080

1081 def __bool__(f):

1082 return not dmp_zero_p(f.rep, f.lev)

1083

1084

1085def init_normal_DMF(num, den, lev, dom):

1086 return DMF(dmp_normal(num, lev, dom),

1087 dmp_normal(den, lev, dom), dom, lev)

1088

1089

1090class DMF(PicklableWithSlots, CantSympify):

1091 """Dense Multivariate Fractions over `K`. """

1092

1093 __slots__ = ('num', 'den', 'lev', 'dom', 'ring')

1094

1095 def __init__(self, rep, dom, lev=None, ring=None):

1096 num, den, lev = self._parse(rep, dom, lev)

1097 num, den = dmp_cancel(num, den, lev, dom)

1098

1099 self.num = num

1100 self.den = den

1101 self.lev = lev

1102 self.dom = dom

1103 self.ring = ring

1104

1105 @classmethod

1106 def new(cls, rep, dom, lev=None, ring=None):

1107 num, den, lev = cls._parse(rep, dom, lev)

1108

1109 obj = object.__new__(cls)

1110

1111 obj.num = num

1112 obj.den = den

1113 obj.lev = lev

1114 obj.dom = dom

1115 obj.ring = ring

1116

1117 return obj

1118

1119 @classmethod

1120 def _parse(cls, rep, dom, lev=None):

1121 if isinstance(rep, tuple):

1122 num, den = rep

1123

1124 if lev is not None:

1125 if isinstance(num, dict):

1126 num = dmp_from_dict(num, lev, dom)

1127

1128 if isinstance(den, dict):

1129 den = dmp_from_dict(den, lev, dom)

1130 else:

1131 num, num_lev = dmp_validate(num)

1132 den, den_lev = dmp_validate(den)

1133

1134 if num_lev == den_lev:

1135 lev = num_lev

1136 else:

1137 raise ValueError('inconsistent number of levels')

1138

1139 if dmp_zero_p(den, lev):

1140 raise ZeroDivisionError('fraction denominator')

1141

1142 if dmp_zero_p(num, lev):

1143 den = dmp_one(lev, dom)

1144 else:

1145 if dmp_negative_p(den, lev, dom):

1146 num = dmp_neg(num, lev, dom)

1147 den = dmp_neg(den, lev, dom)

1148 else:

1149 num = rep

1150

1151 if lev is not None:

1152 if isinstance(num, dict):

1153 num = dmp_from_dict(num, lev, dom)

1154 elif not isinstance(num, list):

1155 num = dmp_ground(dom.convert(num), lev)

1156 else:

1157 num, lev = dmp_validate(num)

1158

1159 den = dmp_one(lev, dom)

1160

1161 return num, den, lev

1162

1163 def __repr__(f):

1164 return "%s((%s, %s), %s, %s)" % (f.__class__.__name__, f.num, f.den,

1165 f.dom, f.ring)

1166

1167 def __hash__(f):

1168 return hash((f.__class__.__name__, dmp_to_tuple(f.num, f.lev),

1169 dmp_to_tuple(f.den, f.lev), f.lev, f.dom, f.ring))

1170

1171 def poly_unify(f, g):

1172 """Unify a multivariate fraction and a polynomial. """

1173 if not isinstance(g, DMP) or f.lev != g.lev:

1174 raise UnificationFailed("Cannot unify %s with %s" % (f, g))

1175

1176 if f.dom == g.dom and f.ring == g.ring:

1177 return (f.lev, f.dom, f.per, (f.num, f.den), g.rep)

1178 else:

1179 lev, dom = f.lev, f.dom.unify(g.dom)

1180 ring = f.ring

1181 if g.ring is not None:

1182 if ring is not None:

1183 ring = ring.unify(g.ring)

1184 else:

1185 ring = g.ring

1186

1187 F = (dmp_convert(f.num, lev, f.dom, dom),

1188 dmp_convert(f.den, lev, f.dom, dom))

1189

1190 G = dmp_convert(g.rep, lev, g.dom, dom)

1191

1192 def per(num, den, cancel=True, kill=False, lev=lev):

1193 if kill:

1194 if not lev:

1195 return num/den

1196 else:

1197 lev = lev - 1

1198

1199 if cancel:

1200 num, den = dmp_cancel(num, den, lev, dom)

1201

1202 return f.__class__.new((num, den), dom, lev, ring=ring)

1203

1204 return lev, dom, per, F, G

1205

1206 def frac_unify(f, g):

1207 """Unify representations of two multivariate fractions. """

1208 if not isinstance(g, DMF) or f.lev != g.lev:

1209 raise UnificationFailed("Cannot unify %s with %s" % (f, g))

1210

1211 if f.dom == g.dom and f.ring == g.ring:

1212 return (f.lev, f.dom, f.per, (f.num, f.den),

1213 (g.num, g.den))

1214 else:

1215 lev, dom = f.lev, f.dom.unify(g.dom)

1216 ring = f.ring

1217 if g.ring is not None:

1218 if ring is not None:

1219 ring = ring.unify(g.ring)

1220 else:

1221 ring = g.ring

1222

1223 F = (dmp_convert(f.num, lev, f.dom, dom),

1224 dmp_convert(f.den, lev, f.dom, dom))

1225

1226 G = (dmp_convert(g.num, lev, g.dom, dom),

1227 dmp_convert(g.den, lev, g.dom, dom))

1228

1229 def per(num, den, cancel=True, kill=False, lev=lev):

1230 if kill:

1231 if not lev:

1232 return num/den

1233 else:

1234 lev = lev - 1

1235

1236 if cancel:

1237 num, den = dmp_cancel(num, den, lev, dom)

1238

1239 return f.__class__.new((num, den), dom, lev, ring=ring)

1240

1241 return lev, dom, per, F, G

1242

1243 def per(f, num, den, cancel=True, kill=False, ring=None):

1244 """Create a DMF out of the given representation. """

1245 lev, dom = f.lev, f.dom

1246

1247 if kill:

1248 if not lev:

1249 return num/den

1250 else:

1251 lev -= 1

1252

1253 if cancel:

1254 num, den = dmp_cancel(num, den, lev, dom)

1255

1256 if ring is None:

1257 ring = f.ring

1258

1259 return f.__class__.new((num, den), dom, lev, ring=ring)

1260

1261 def half_per(f, rep, kill=False):

1262 """Create a DMP out of the given representation. """

1263 lev = f.lev

1264

1265 if kill:

1266 if not lev:

1267 return rep

1268 else:

1269 lev -= 1

1270

1271 return DMP(rep, f.dom, lev)

1272

1273 @classmethod

1274 def zero(cls, lev, dom, ring=None):

1275 return cls.new(0, dom, lev, ring=ring)

1276

1277 @classmethod

1278 def one(cls, lev, dom, ring=None):

1279 return cls.new(1, dom, lev, ring=ring)

1280

1281 def numer(f):

1282 """Returns the numerator of ``f``. """

1283 return f.half_per(f.num)

1284

1285 def denom(f):

1286 """Returns the denominator of ``f``. """

1287 return f.half_per(f.den)

1288

1289 def cancel(f):

1290 """Remove common factors from ``f.num`` and ``f.den``. """

1291 return f.per(f.num, f.den)

1292

1293 def neg(f):

1294 """Negate all coefficients in ``f``. """

1295 return f.per(dmp_neg(f.num, f.lev, f.dom), f.den, cancel=False)

1296

1297 def add(f, g):

1298 """Add two multivariate fractions ``f`` and ``g``. """

1299 if isinstance(g, DMP):

1300 lev, dom, per, (F_num, F_den), G = f.poly_unify(g)

1301 num, den = dmp_add_mul(F_num, F_den, G, lev, dom), F_den

1302 else:

1303 lev, dom, per, F, G = f.frac_unify(g)

1304 (F_num, F_den), (G_num, G_den) = F, G

1305

1306 num = dmp_add(dmp_mul(F_num, G_den, lev, dom),

1307 dmp_mul(F_den, G_num, lev, dom), lev, dom)

1308 den = dmp_mul(F_den, G_den, lev, dom)

1309

1310 return per(num, den)

1311

1312 def sub(f, g):

1313 """Subtract two multivariate fractions ``f`` and ``g``. """

1314 if isinstance(g, DMP):

1315 lev, dom, per, (F_num, F_den), G = f.poly_unify(g)

1316 num, den = dmp_sub_mul(F_num, F_den, G, lev, dom), F_den

1317 else:

1318 lev, dom, per, F, G = f.frac_unify(g)

1319 (F_num, F_den), (G_num, G_den) = F, G

1320

1321 num = dmp_sub(dmp_mul(F_num, G_den, lev, dom),

1322 dmp_mul(F_den, G_num, lev, dom), lev, dom)

1323 den = dmp_mul(F_den, G_den, lev, dom)

1324

1325 return per(num, den)

1326

1327 def mul(f, g):

1328 """Multiply two multivariate fractions ``f`` and ``g``. """

1329 if isinstance(g, DMP):

1330 lev, dom, per, (F_num, F_den), G = f.poly_unify(g)

1331 num, den = dmp_mul(F_num, G, lev, dom), F_den

1332 else:

1333 lev, dom, per, F, G = f.frac_unify(g)

1334 (F_num, F_den), (G_num, G_den) = F, G

1335

1336 num = dmp_mul(F_num, G_num, lev, dom)

1337 den = dmp_mul(F_den, G_den, lev, dom)

1338

1339 return per(num, den)

1340

1341 def pow(f, n):

1342 """Raise ``f`` to a non-negative power ``n``. """

1343 if isinstance(n, int):

1344 num, den = f.num, f.den

1345 if n < 0:

1346 num, den, n = den, num, -n

1347 return f.per(dmp_pow(num, n, f.lev, f.dom),

1348 dmp_pow(den, n, f.lev, f.dom), cancel=False)

1349 else:

1350 raise TypeError("``int`` expected, got %s" % type(n))

1351

1352 def quo(f, g):

1353 """Computes quotient of fractions ``f`` and ``g``. """

1354 if isinstance(g, DMP):

1355 lev, dom, per, (F_num, F_den), G = f.poly_unify(g)

1356 num, den = F_num, dmp_mul(F_den, G, lev, dom)

1357 else:

1358 lev, dom, per, F, G = f.frac_unify(g)

1359 (F_num, F_den), (G_num, G_den) = F, G

1360

1361 num = dmp_mul(F_num, G_den, lev, dom)

1362 den = dmp_mul(F_den, G_num, lev, dom)

1363

1364 res = per(num, den)

1365 if f.ring is not None and res not in f.ring:

1366 from sympy.polys.polyerrors import ExactQuotientFailed

1367 raise ExactQuotientFailed(f, g, f.ring)

1368 return res

1369

1370 exquo = quo

1371

1372 def invert(f, check=True):

1373 """Computes inverse of a fraction ``f``. """

1374 if check and f.ring is not None and not f.ring.is_unit(f):

1375 raise NotReversible(f, f.ring)

1376 res = f.per(f.den, f.num, cancel=False)

1377 return res

1378

1379 @property

1380 def is_zero(f):

1381 """Returns ``True`` if ``f`` is a zero fraction. """

1382 return dmp_zero_p(f.num, f.lev)

1383

1384 @property

1385 def is_one(f):

1386 """Returns ``True`` if ``f`` is a unit fraction. """

1387 return dmp_one_p(f.num, f.lev, f.dom) and \

1388 dmp_one_p(f.den, f.lev, f.dom)

1389

1390 def __neg__(f):

1391 return f.neg()

1392

1393 def __add__(f, g):

1394 if isinstance(g, (DMP, DMF)):

1395 return f.add(g)

1396

1397 try:

1398 return f.add(f.half_per(g))

1399 except TypeError:

1400 return NotImplemented

1401 except (CoercionFailed, NotImplementedError):

1402 if f.ring is not None:

1403 try:

1404 return f.add(f.ring.convert(g))

1405 except (CoercionFailed, NotImplementedError):

1406 pass

1407 return NotImplemented

1408

1409 def __radd__(f, g):

1410 return f.__add__(g)

1411

1412 def __sub__(f, g):

1413 if isinstance(g, (DMP, DMF)):

1414 return f.sub(g)

1415

1416 try:

1417 return f.sub(f.half_per(g))

1418 except TypeError:

1419 return NotImplemented

1420 except (CoercionFailed, NotImplementedError):

1421 if f.ring is not None:

1422 try:

1423 return f.sub(f.ring.convert(g))

1424 except (CoercionFailed, NotImplementedError):

1425 pass

1426 return NotImplemented

1427

1428 def __rsub__(f, g):

1429 return (-f).__add__(g)

1430

1431 def __mul__(f, g):

1432 if isinstance(g, (DMP, DMF)):

1433 return f.mul(g)

1434

1435 try:

1436 return f.mul(f.half_per(g))

1437 except TypeError:

1438 return NotImplemented

1439 except (CoercionFailed, NotImplementedError):

1440 if f.ring is not None:

1441 try:

1442 return f.mul(f.ring.convert(g))

1443 except (CoercionFailed, NotImplementedError):

1444 pass

1445 return NotImplemented

1446

1447 def __rmul__(f, g):

1448 return f.__mul__(g)

1449

1450 def __pow__(f, n):

1451 return f.pow(n)

1452

1453 def __truediv__(f, g):

1454 if isinstance(g, (DMP, DMF)):

1455 return f.quo(g)

1456

1457 try:

1458 return f.quo(f.half_per(g))

1459 except TypeError:

1460 return NotImplemented

1461 except (CoercionFailed, NotImplementedError):

1462 if f.ring is not None:

1463 try:

1464 return f.quo(f.ring.convert(g))

1465 except (CoercionFailed, NotImplementedError):

1466 pass

1467 return NotImplemented

1468

1469 def __rtruediv__(self, g):

1470 r = self.invert(check=False)*g

1471 if self.ring and r not in self.ring:

1472 from sympy.polys.polyerrors import ExactQuotientFailed

1473 raise ExactQuotientFailed(g, self, self.ring)

1474 return r

1475

1476 def __eq__(f, g):

1477 try:

1478 if isinstance(g, DMP):

1479 _, _, _, (F_num, F_den), G = f.poly_unify(g)

1480

1481 if f.lev == g.lev:

1482 return dmp_one_p(F_den, f.lev, f.dom) and F_num == G

1483 else:

1484 _, _, _, F, G = f.frac_unify(g)

1485

1486 if f.lev == g.lev:

1487 return F == G

1488 except UnificationFailed:

1489 pass

1490

1491 return False

1492

1493 def __ne__(f, g):

1494 try:

1495 if isinstance(g, DMP):

1496 _, _, _, (F_num, F_den), G = f.poly_unify(g)

1497

1498 if f.lev == g.lev:

1499 return not (dmp_one_p(F_den, f.lev, f.dom) and F_num == G)

1500 else:

1501 _, _, _, F, G = f.frac_unify(g)

1502

1503 if f.lev == g.lev:

1504 return F != G

1505 except UnificationFailed:

1506 pass

1507

1508 return True

1509

1510 def __lt__(f, g):

1511 _, _, _, F, G = f.frac_unify(g)

1512 return F < G

1513

1514 def __le__(f, g):

1515 _, _, _, F, G = f.frac_unify(g)

1516 return F <= G

1517

1518 def __gt__(f, g):

1519 _, _, _, F, G = f.frac_unify(g)

1520 return F > G

1521

1522 def __ge__(f, g):

1523 _, _, _, F, G = f.frac_unify(g)

1524 return F >= G

1525

1526 def __bool__(f):

1527 return not dmp_zero_p(f.num, f.lev)

1528

1529

1530def init_normal_ANP(rep, mod, dom):

1531 return ANP(dup_normal(rep, dom),

1532 dup_normal(mod, dom), dom)

1533

1534

1535class ANP(PicklableWithSlots, CantSympify):

1536 """Dense Algebraic Number Polynomials over a field. """

1537

1538 __slots__ = ('rep', 'mod', 'dom')

1539

1540 def __init__(self, rep, mod, dom):

1541 # Not possible to check with isinstance

1542 if type(rep) is dict:

1543 self.rep = dup_from_dict(rep, dom)

1544 else:

1545 if isinstance(rep, list):

1546 rep = [dom.convert(a) for a in rep]

1547 else:

1548 rep = [dom.convert(rep)]

1549

1550 self.rep = dup_strip(rep)

1551

1552 if isinstance(mod, DMP):

1553 self.mod = mod.rep

1554 else:

1555 if isinstance(mod, dict):

1556 self.mod = dup_from_dict(mod, dom)

1557 else:

1558 self.mod = dup_strip(mod)

1559

1560 self.dom = dom

1561

1562 def __repr__(f):

1563 return "%s(%s, %s, %s)" % (f.__class__.__name__, f.rep, f.mod, f.dom)

1564

1565 def __hash__(f):

1566 return hash((f.__class__.__name__, f.to_tuple(), dmp_to_tuple(f.mod, 0), f.dom))

1567

1568 def unify(f, g):

1569 """Unify representations of two algebraic numbers. """

1570 if not isinstance(g, ANP) or f.mod != g.mod:

1571 raise UnificationFailed("Cannot unify %s with %s" % (f, g))

1572

1573 if f.dom == g.dom:

1574 return f.dom, f.per, f.rep, g.rep, f.mod

1575 else:

1576 dom = f.dom.unify(g.dom)

1577

1578 F = dup_convert(f.rep, f.dom, dom)

1579 G = dup_convert(g.rep, g.dom, dom)

1580

1581 if dom != f.dom and dom != g.dom:

1582 mod = dup_convert(f.mod, f.dom, dom)

1583 else:

1584 if dom == f.dom:

1585 mod = f.mod

1586 else:

1587 mod = g.mod

1588

1589 per = lambda rep: ANP(rep, mod, dom)

1590

1591 return dom, per, F, G, mod

1592

1593 def per(f, rep, mod=None, dom=None):

1594 return ANP(rep, mod or f.mod, dom or f.dom)

1595

1596 @classmethod

1597 def zero(cls, mod, dom):

1598 return ANP(0, mod, dom)

1599

1600 @classmethod

1601 def one(cls, mod, dom):

1602 return ANP(1, mod, dom)

1603

1604 def to_dict(f):

1605 """Convert ``f`` to a dict representation with native coefficients. """

1606 return dmp_to_dict(f.rep, 0, f.dom)

1607

1608 def to_sympy_dict(f):

1609 """Convert ``f`` to a dict representation with SymPy coefficients. """

1610 rep = dmp_to_dict(f.rep, 0, f.dom)

1611

1612 for k, v in rep.items():

1613 rep[k] = f.dom.to_sympy(v)

1614

1615 return rep

1616

1617 def to_list(f):

1618 """Convert ``f`` to a list representation with native coefficients. """

1619 return f.rep

1620

1621 def to_sympy_list(f):

1622 """Convert ``f`` to a list representation with SymPy coefficients. """

1623 return [ f.dom.to_sympy(c) for c in f.rep ]

1624

1625 def to_tuple(f):

1626 """

1627 Convert ``f`` to a tuple representation with native coefficients.

1628

1629 This is needed for hashing.

1630 """

1631 return dmp_to_tuple(f.rep, 0)

1632

1633 @classmethod

1634 def from_list(cls, rep, mod, dom):

1635 return ANP(dup_strip(list(map(dom.convert, rep))), mod, dom)

1636

1637 def neg(f):

1638 return f.per(dup_neg(f.rep, f.dom))

1639

1640 def add(f, g):

1641 dom, per, F, G, mod = f.unify(g)

1642 return per(dup_add(F, G, dom))

1643

1644 def sub(f, g):

1645 dom, per, F, G, mod = f.unify(g)

1646 return per(dup_sub(F, G, dom))

1647

1648 def mul(f, g):

1649 dom, per, F, G, mod = f.unify(g)

1650 return per(dup_rem(dup_mul(F, G, dom), mod, dom))

1651

1652 def pow(f, n):

1653 """Raise ``f`` to a non-negative power ``n``. """

1654 if isinstance(n, int):

1655 if n < 0:

1656 F, n = dup_invert(f.rep, f.mod, f.dom), -n

1657 else:

1658 F = f.rep

1659

1660 return f.per(dup_rem(dup_pow(F, n, f.dom), f.mod, f.dom))

1661 else:

1662 raise TypeError("``int`` expected, got %s" % type(n))

1663

1664 def div(f, g):

1665 dom, per, F, G, mod = f.unify(g)

1666 return (per(dup_rem(dup_mul(F, dup_invert(G, mod, dom), dom), mod, dom)), f.zero(mod, dom))

1667

1668 def rem(f, g):

1669 dom, _, _, G, mod = f.unify(g)

1670

1671 s, h = dup_half_gcdex(G, mod, dom)

1672

1673 if h == [dom.one]:

1674 return f.zero(mod, dom)

1675 else:

1676 raise NotInvertible("zero divisor")

1677

1678 def quo(f, g):

1679 dom, per, F, G, mod = f.unify(g)

1680 return per(dup_rem(dup_mul(F, dup_invert(G, mod, dom), dom), mod, dom))

1681

1682 exquo = quo

1683

1684 def LC(f):

1685 """Returns the leading coefficient of ``f``. """

1686 return dup_LC(f.rep, f.dom)

1687

1688 def TC(f):

1689 """Returns the trailing coefficient of ``f``. """

1690 return dup_TC(f.rep, f.dom)

1691

1692 @property

1693 def is_zero(f):

1694 """Returns ``True`` if ``f`` is a zero algebraic number. """

1695 return not f

1696

1697 @property

1698 def is_one(f):

1699 """Returns ``True`` if ``f`` is a unit algebraic number. """

1700 return f.rep == [f.dom.one]

1701

1702 @property

1703 def is_ground(f):

1704 """Returns ``True`` if ``f`` is an element of the ground domain. """

1705 return not f.rep or len(f.rep) == 1

1706

1707 def __pos__(f):

1708 return f

1709

1710 def __neg__(f):

1711 return f.neg()

1712

1713 def __add__(f, g):

1714 if isinstance(g, ANP):

1715 return f.add(g)

1716 else:

1717 try:

1718 return f.add(f.per(g))

1719 except (CoercionFailed, TypeError):

1720 return NotImplemented

1721

1722 def __radd__(f, g):

1723 return f.__add__(g)

1724

1725 def __sub__(f, g):

1726 if isinstance(g, ANP):

1727 return f.sub(g)

1728 else:

1729 try:

1730 return f.sub(f.per(g))

1731 except (CoercionFailed, TypeError):

1732 return NotImplemented

1733

1734 def __rsub__(f, g):

1735 return (-f).__add__(g)

1736

1737 def __mul__(f, g):

1738 if isinstance(g, ANP):

1739 return f.mul(g)

1740 else:

1741 try:

1742 return f.mul(f.per(g))

1743 except (CoercionFailed, TypeError):

1744 return NotImplemented

1745

1746 def __rmul__(f, g):

1747 return f.__mul__(g)

1748

1749 def __pow__(f, n):

1750 return f.pow(n)

1751

1752 def __divmod__(f, g):

1753 return f.div(g)

1754

1755 def __mod__(f, g):

1756 return f.rem(g)

1757

1758 def __truediv__(f, g):

1759 if isinstance(g, ANP):

1760 return f.quo(g)

1761 else:

1762 try:

1763 return f.quo(f.per(g))

1764 except (CoercionFailed, TypeError):

1765 return NotImplemented

1766

1767 def __eq__(f, g):

1768 try:

1769 _, _, F, G, _ = f.unify(g)

1770

1771 return F == G

1772 except UnificationFailed:

1773 return False

1774

1775 def __ne__(f, g):

1776 try:

1777 _, _, F, G, _ = f.unify(g)

1778

1779 return F != G

1780 except UnificationFailed:

1781 return True

1782

1783 def __lt__(f, g):

1784 _, _, F, G, _ = f.unify(g)

1785 return F < G

1786

1787 def __le__(f, g):

1788 _, _, F, G, _ = f.unify(g)

1789 return F <= G

1790

1791 def __gt__(f, g):

1792 _, _, F, G, _ = f.unify(g)

1793 return F > G

1794

1795 def __ge__(f, g):

1796 _, _, F, G, _ = f.unify(g)

1797 return F >= G

1798

1799 def __bool__(f):

1800 return bool(f.rep)