puiseux.py 26 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795
  1. """
  2. Puiseux rings. These are used by the ring_series module to represented
  3. truncated Puiseux series. Elements of a Puiseux ring are like polynomials
  4. except that the exponents can be negative or rational rather than just
  5. non-negative integers.
  6. """
  7. # Previously the ring_series module used PolyElement to represent Puiseux
  8. # series. This is problematic because it means that PolyElement has to support
  9. # negative and non-integer exponents which most polynomial representations do
  10. # not support. This module provides an implementation of a ring for Puiseux
  11. # series that can be used by ring_series without breaking the basic invariants
  12. # of polynomial rings.
  13. #
  14. # Ideally there would be more of a proper series type that can keep track of
  15. # not just the leading terms of a truncated series but also the precision
  16. # of the series. For now the rings here are just introduced to keep the
  17. # interface that ring_series was using before.
  18. from __future__ import annotations
  19. from sympy.polys.domains import QQ
  20. from sympy.polys.rings import PolyRing, PolyElement
  21. from sympy.core.add import Add
  22. from sympy.core.mul import Mul
  23. from sympy.external.gmpy import gcd, lcm
  24. from typing import TYPE_CHECKING
  25. if TYPE_CHECKING:
  26. from typing import Any, Unpack
  27. from sympy.core.expr import Expr
  28. from sympy.polys.domains import Domain
  29. from collections.abc import Iterable, Iterator
  30. def puiseux_ring(
  31. symbols: str | list[Expr], domain: Domain
  32. ) -> tuple[PuiseuxRing, Unpack[tuple[PuiseuxPoly, ...]]]:
  33. """Construct a Puiseux ring.
  34. This function constructs a Puiseux ring with the given symbols and domain.
  35. >>> from sympy.polys.domains import QQ
  36. >>> from sympy.polys.puiseux import puiseux_ring
  37. >>> R, x, y = puiseux_ring('x y', QQ)
  38. >>> R
  39. PuiseuxRing((x, y), QQ)
  40. >>> p = 5*x**QQ(1,2) + 7/y
  41. >>> p
  42. 7*y**(-1) + 5*x**(1/2)
  43. """
  44. ring = PuiseuxRing(symbols, domain)
  45. return (ring,) + ring.gens # type: ignore
  46. class PuiseuxRing:
  47. """Ring of Puiseux polynomials.
  48. A Puiseux polynomial is a truncated Puiseux series. The exponents of the
  49. monomials can be negative or rational numbers. This ring is used by the
  50. ring_series module:
  51. >>> from sympy.polys.domains import QQ
  52. >>> from sympy.polys.puiseux import puiseux_ring
  53. >>> from sympy.polys.ring_series import rs_exp, rs_nth_root
  54. >>> ring, x, y = puiseux_ring('x y', QQ)
  55. >>> f = x**2 + y**3
  56. >>> f
  57. y**3 + x**2
  58. >>> f.diff(x)
  59. 2*x
  60. >>> rs_exp(x, x, 5)
  61. 1 + x + 1/2*x**2 + 1/6*x**3 + 1/24*x**4
  62. Importantly the Puiseux ring can represent truncated series with negative
  63. and fractional exponents:
  64. >>> f = 1/x + 1/y**2
  65. >>> f
  66. x**(-1) + y**(-2)
  67. >>> f.diff(x)
  68. -1*x**(-2)
  69. >>> rs_nth_root(8*x + x**2 + x**3, 3, x, 5)
  70. 2*x**(1/3) + 1/12*x**(4/3) + 23/288*x**(7/3) + -139/20736*x**(10/3)
  71. See Also
  72. ========
  73. sympy.polys.ring_series.rs_series
  74. PuiseuxPoly
  75. """
  76. def __init__(self, symbols: str | list[Expr], domain: Domain):
  77. poly_ring = PolyRing(symbols, domain)
  78. domain = poly_ring.domain
  79. ngens = poly_ring.ngens
  80. self.poly_ring = poly_ring
  81. self.domain = domain
  82. self.symbols = poly_ring.symbols
  83. self.gens = tuple([self.from_poly(g) for g in poly_ring.gens])
  84. self.ngens = ngens
  85. self.zero = self.from_poly(poly_ring.zero)
  86. self.one = self.from_poly(poly_ring.one)
  87. self.zero_monom = poly_ring.zero_monom # type: ignore
  88. self.monomial_mul = poly_ring.monomial_mul # type: ignore
  89. def __repr__(self) -> str:
  90. return f"PuiseuxRing({self.symbols}, {self.domain})"
  91. def __eq__(self, other: Any) -> bool:
  92. if not isinstance(other, PuiseuxRing):
  93. return NotImplemented
  94. return self.symbols == other.symbols and self.domain == other.domain
  95. def from_poly(self, poly: PolyElement) -> PuiseuxPoly:
  96. """Create a Puiseux polynomial from a polynomial.
  97. >>> from sympy.polys.domains import QQ
  98. >>> from sympy.polys.rings import ring
  99. >>> from sympy.polys.puiseux import puiseux_ring
  100. >>> R1, x1 = ring('x', QQ)
  101. >>> R2, x2 = puiseux_ring('x', QQ)
  102. >>> R2.from_poly(x1**2)
  103. x**2
  104. """
  105. return PuiseuxPoly(poly, self)
  106. def from_dict(self, terms: dict[tuple[int, ...], Any]) -> PuiseuxPoly:
  107. """Create a Puiseux polynomial from a dictionary of terms.
  108. >>> from sympy.polys.domains import QQ
  109. >>> from sympy.polys.puiseux import puiseux_ring
  110. >>> R, x = puiseux_ring('x', QQ)
  111. >>> R.from_dict({(QQ(1,2),): QQ(3)})
  112. 3*x**(1/2)
  113. """
  114. return PuiseuxPoly.from_dict(terms, self)
  115. def from_int(self, n: int) -> PuiseuxPoly:
  116. """Create a Puiseux polynomial from an integer.
  117. >>> from sympy.polys.domains import QQ
  118. >>> from sympy.polys.puiseux import puiseux_ring
  119. >>> R, x = puiseux_ring('x', QQ)
  120. >>> R.from_int(3)
  121. 3
  122. """
  123. return self.from_poly(self.poly_ring(n))
  124. def domain_new(self, arg: Any) -> Any:
  125. """Create a new element of the domain.
  126. >>> from sympy.polys.domains import QQ
  127. >>> from sympy.polys.puiseux import puiseux_ring
  128. >>> R, x = puiseux_ring('x', QQ)
  129. >>> R.domain_new(3)
  130. 3
  131. >>> QQ.of_type(_)
  132. True
  133. """
  134. return self.poly_ring.domain_new(arg)
  135. def ground_new(self, arg: Any) -> PuiseuxPoly:
  136. """Create a new element from a ground element.
  137. >>> from sympy.polys.domains import QQ
  138. >>> from sympy.polys.puiseux import puiseux_ring, PuiseuxPoly
  139. >>> R, x = puiseux_ring('x', QQ)
  140. >>> R.ground_new(3)
  141. 3
  142. >>> isinstance(_, PuiseuxPoly)
  143. True
  144. """
  145. return self.from_poly(self.poly_ring.ground_new(arg))
  146. def __call__(self, arg: Any) -> PuiseuxPoly:
  147. """Coerce an element into the ring.
  148. >>> from sympy.polys.domains import QQ
  149. >>> from sympy.polys.puiseux import puiseux_ring
  150. >>> R, x = puiseux_ring('x', QQ)
  151. >>> R(3)
  152. 3
  153. >>> R({(QQ(1,2),): QQ(3)})
  154. 3*x**(1/2)
  155. """
  156. if isinstance(arg, dict):
  157. return self.from_dict(arg)
  158. else:
  159. return self.from_poly(self.poly_ring(arg))
  160. def index(self, x: PuiseuxPoly) -> int:
  161. """Return the index of a generator.
  162. >>> from sympy.polys.domains import QQ
  163. >>> from sympy.polys.puiseux import puiseux_ring
  164. >>> R, x, y = puiseux_ring('x y', QQ)
  165. >>> R.index(x)
  166. 0
  167. >>> R.index(y)
  168. 1
  169. """
  170. return self.gens.index(x)
  171. def _div_poly_monom(poly: PolyElement, monom: Iterable[int]) -> PolyElement:
  172. ring = poly.ring
  173. div = ring.monomial_div
  174. return ring.from_dict({div(m, monom): c for m, c in poly.terms()})
  175. def _mul_poly_monom(poly: PolyElement, monom: Iterable[int]) -> PolyElement:
  176. ring = poly.ring
  177. mul = ring.monomial_mul
  178. return ring.from_dict({mul(m, monom): c for m, c in poly.terms()})
  179. def _div_monom(monom: Iterable[int], div: Iterable[int]) -> tuple[int, ...]:
  180. return tuple(mi - di for mi, di in zip(monom, div))
  181. class PuiseuxPoly:
  182. """Puiseux polynomial. Represents a truncated Puiseux series.
  183. See the :class:`PuiseuxRing` class for more information.
  184. >>> from sympy import QQ
  185. >>> from sympy.polys.puiseux import puiseux_ring
  186. >>> R, x, y = puiseux_ring('x, y', QQ)
  187. >>> p = 5*x**2 + 7*y**3
  188. >>> p
  189. 7*y**3 + 5*x**2
  190. The internal representation of a Puiseux polynomial wraps a normal
  191. polynomial. To support negative powers the polynomial is considered to be
  192. divided by a monomial.
  193. >>> p2 = 1/x + 1/y**2
  194. >>> p2.monom # x*y**2
  195. (1, 2)
  196. >>> p2.poly
  197. x + y**2
  198. >>> (y**2 + x) / (x*y**2) == p2
  199. True
  200. To support fractional powers the polynomial is considered to be a function
  201. of ``x**(1/nx), y**(1/ny), ...``. The representation keeps track of a
  202. monomial and a list of exponent denominators so that the polynomial can be
  203. used to represent both negative and fractional powers.
  204. >>> p3 = x**QQ(1,2) + y**QQ(2,3)
  205. >>> p3.ns
  206. (2, 3)
  207. >>> p3.poly
  208. x + y**2
  209. See Also
  210. ========
  211. sympy.polys.puiseux.PuiseuxRing
  212. sympy.polys.rings.PolyElement
  213. """
  214. ring: PuiseuxRing
  215. poly: PolyElement
  216. monom: tuple[int, ...] | None
  217. ns: tuple[int, ...] | None
  218. def __new__(cls, poly: PolyElement, ring: PuiseuxRing) -> PuiseuxPoly:
  219. return cls._new(ring, poly, None, None)
  220. @classmethod
  221. def _new(
  222. cls,
  223. ring: PuiseuxRing,
  224. poly: PolyElement,
  225. monom: tuple[int, ...] | None,
  226. ns: tuple[int, ...] | None,
  227. ) -> PuiseuxPoly:
  228. poly, monom, ns = cls._normalize(poly, monom, ns)
  229. return cls._new_raw(ring, poly, monom, ns)
  230. @classmethod
  231. def _new_raw(
  232. cls,
  233. ring: PuiseuxRing,
  234. poly: PolyElement,
  235. monom: tuple[int, ...] | None,
  236. ns: tuple[int, ...] | None,
  237. ) -> PuiseuxPoly:
  238. obj = object.__new__(cls)
  239. obj.ring = ring
  240. obj.poly = poly
  241. obj.monom = monom
  242. obj.ns = ns
  243. return obj
  244. def __eq__(self, other: Any) -> bool:
  245. if isinstance(other, PuiseuxPoly):
  246. return (
  247. self.poly == other.poly
  248. and self.monom == other.monom
  249. and self.ns == other.ns
  250. )
  251. elif self.monom is None and self.ns is None:
  252. return self.poly.__eq__(other)
  253. else:
  254. return NotImplemented
  255. @classmethod
  256. def _normalize(
  257. cls,
  258. poly: PolyElement,
  259. monom: tuple[int, ...] | None,
  260. ns: tuple[int, ...] | None,
  261. ) -> tuple[PolyElement, tuple[int, ...] | None, tuple[int, ...] | None]:
  262. if monom is None and ns is None:
  263. return poly, None, None
  264. if monom is not None:
  265. degs = [max(d, 0) for d in poly.tail_degrees()]
  266. if all(di >= mi for di, mi in zip(degs, monom)):
  267. poly = _div_poly_monom(poly, monom)
  268. monom = None
  269. elif any(degs):
  270. poly = _div_poly_monom(poly, degs)
  271. monom = _div_monom(monom, degs)
  272. if ns is not None:
  273. factors_d, [poly_d] = poly.deflate()
  274. degrees = poly.degrees()
  275. monom_d = monom if monom is not None else [0] * len(degrees)
  276. ns_new = []
  277. monom_new = []
  278. inflations = []
  279. for fi, ni, di, mi in zip(factors_d, ns, degrees, monom_d):
  280. if di == 0:
  281. g = gcd(ni, mi)
  282. else:
  283. g = gcd(fi, ni, mi)
  284. ns_new.append(ni // g)
  285. monom_new.append(mi // g)
  286. inflations.append(fi // g)
  287. if any(infl > 1 for infl in inflations):
  288. poly_d = poly_d.inflate(inflations)
  289. poly = poly_d
  290. if monom is not None:
  291. monom = tuple(monom_new)
  292. if all(n == 1 for n in ns_new):
  293. ns = None
  294. else:
  295. ns = tuple(ns_new)
  296. return poly, monom, ns
  297. @classmethod
  298. def _monom_fromint(
  299. cls,
  300. monom: tuple[int, ...],
  301. dmonom: tuple[int, ...] | None,
  302. ns: tuple[int, ...] | None,
  303. ) -> tuple[Any, ...]:
  304. if dmonom is not None and ns is not None:
  305. return tuple(QQ(mi - di, ni) for mi, di, ni in zip(monom, dmonom, ns))
  306. elif dmonom is not None:
  307. return tuple(QQ(mi - di) for mi, di in zip(monom, dmonom))
  308. elif ns is not None:
  309. return tuple(QQ(mi, ni) for mi, ni in zip(monom, ns))
  310. else:
  311. return tuple(QQ(mi) for mi in monom)
  312. @classmethod
  313. def _monom_toint(
  314. cls,
  315. monom: tuple[Any, ...],
  316. dmonom: tuple[int, ...] | None,
  317. ns: tuple[int, ...] | None,
  318. ) -> tuple[int, ...]:
  319. if dmonom is not None and ns is not None:
  320. return tuple(
  321. int((mi * ni).numerator + di) for mi, di, ni in zip(monom, dmonom, ns)
  322. )
  323. elif dmonom is not None:
  324. return tuple(int(mi.numerator + di) for mi, di in zip(monom, dmonom))
  325. elif ns is not None:
  326. return tuple(int((mi * ni).numerator) for mi, ni in zip(monom, ns))
  327. else:
  328. return tuple(int(mi.numerator) for mi in monom)
  329. def itermonoms(self) -> Iterator[tuple[Any, ...]]:
  330. """Iterate over the monomials of a Puiseux polynomial.
  331. >>> from sympy import QQ
  332. >>> from sympy.polys.puiseux import puiseux_ring
  333. >>> R, x, y = puiseux_ring('x, y', QQ)
  334. >>> p = 5*x**2 + 7*y**3
  335. >>> list(p.itermonoms())
  336. [(2, 0), (0, 3)]
  337. >>> p[(2, 0)]
  338. 5
  339. """
  340. monom, ns = self.monom, self.ns
  341. for m in self.poly.itermonoms():
  342. yield self._monom_fromint(m, monom, ns)
  343. def monoms(self) -> list[tuple[Any, ...]]:
  344. """Return a list of the monomials of a Puiseux polynomial."""
  345. return list(self.itermonoms())
  346. def __iter__(self) -> Iterator[tuple[tuple[Any, ...], Any]]:
  347. return self.itermonoms()
  348. def __getitem__(self, monom: tuple[int, ...]) -> Any:
  349. monom = self._monom_toint(monom, self.monom, self.ns)
  350. return self.poly[monom]
  351. def __len__(self) -> int:
  352. return len(self.poly)
  353. def iterterms(self) -> Iterator[tuple[tuple[Any, ...], Any]]:
  354. """Iterate over the terms of a Puiseux polynomial.
  355. >>> from sympy import QQ
  356. >>> from sympy.polys.puiseux import puiseux_ring
  357. >>> R, x, y = puiseux_ring('x, y', QQ)
  358. >>> p = 5*x**2 + 7*y**3
  359. >>> list(p.iterterms())
  360. [((2, 0), 5), ((0, 3), 7)]
  361. """
  362. monom, ns = self.monom, self.ns
  363. for m, coeff in self.poly.iterterms():
  364. mq = self._monom_fromint(m, monom, ns)
  365. yield mq, coeff
  366. def terms(self) -> list[tuple[tuple[Any, ...], Any]]:
  367. """Return a list of the terms of a Puiseux polynomial."""
  368. return list(self.iterterms())
  369. @property
  370. def is_term(self) -> bool:
  371. """Return True if the Puiseux polynomial is a single term."""
  372. return self.poly.is_term
  373. def to_dict(self) -> dict[tuple[int, ...], Any]:
  374. """Return a dictionary representation of a Puiseux polynomial."""
  375. return dict(self.iterterms())
  376. @classmethod
  377. def from_dict(
  378. cls, terms: dict[tuple[Any, ...], Any], ring: PuiseuxRing
  379. ) -> PuiseuxPoly:
  380. """Create a Puiseux polynomial from a dictionary of terms.
  381. >>> from sympy import QQ
  382. >>> from sympy.polys.puiseux import puiseux_ring, PuiseuxPoly
  383. >>> R, x = puiseux_ring('x', QQ)
  384. >>> PuiseuxPoly.from_dict({(QQ(1,2),): QQ(3)}, R)
  385. 3*x**(1/2)
  386. >>> R.from_dict({(QQ(1,2),): QQ(3)})
  387. 3*x**(1/2)
  388. """
  389. ns = [1] * ring.ngens
  390. mon = [0] * ring.ngens
  391. for mo in terms:
  392. ns = [lcm(n, m.denominator) for n, m in zip(ns, mo)]
  393. mon = [min(m, n) for m, n in zip(mo, mon)]
  394. if not any(mon):
  395. monom = None
  396. else:
  397. monom = tuple(-int((m * n).numerator) for m, n in zip(mon, ns))
  398. if all(n == 1 for n in ns):
  399. ns_final = None
  400. else:
  401. ns_final = tuple(ns)
  402. terms_p = {cls._monom_toint(m, monom, ns_final): coeff for m, coeff in terms.items()}
  403. poly = ring.poly_ring.from_dict(terms_p)
  404. return cls._new(ring, poly, monom, ns_final)
  405. def as_expr(self) -> Expr:
  406. """Convert a Puiseux polynomial to :class:`~sympy.core.expr.Expr`.
  407. >>> from sympy import QQ, Expr
  408. >>> from sympy.polys.puiseux import puiseux_ring
  409. >>> R, x = puiseux_ring('x', QQ)
  410. >>> p = 5*x**2 + 7*x**3
  411. >>> p.as_expr()
  412. 7*x**3 + 5*x**2
  413. >>> isinstance(_, Expr)
  414. True
  415. """
  416. ring = self.ring
  417. dom = ring.domain
  418. symbols = ring.symbols
  419. terms = []
  420. for monom, coeff in self.iterterms():
  421. coeff_expr = dom.to_sympy(coeff)
  422. monoms_expr = []
  423. for i, m in enumerate(monom):
  424. monoms_expr.append(symbols[i] ** m)
  425. terms.append(Mul(coeff_expr, *monoms_expr))
  426. return Add(*terms)
  427. def __repr__(self) -> str:
  428. def format_power(base: str, exp: int) -> str:
  429. if exp == 1:
  430. return base
  431. elif exp >= 0 and int(exp) == exp:
  432. return f"{base}**{exp}"
  433. else:
  434. return f"{base}**({exp})"
  435. ring = self.ring
  436. dom = ring.domain
  437. syms = [str(s) for s in ring.symbols]
  438. terms_str = []
  439. for monom, coeff in sorted(self.terms()):
  440. monom_str = "*".join(format_power(s, e) for s, e in zip(syms, monom) if e)
  441. if coeff == dom.one:
  442. if monom_str:
  443. terms_str.append(monom_str)
  444. else:
  445. terms_str.append("1")
  446. elif not monom_str:
  447. terms_str.append(str(coeff))
  448. else:
  449. terms_str.append(f"{coeff}*{monom_str}")
  450. return " + ".join(terms_str)
  451. def _unify(
  452. self, other: PuiseuxPoly
  453. ) -> tuple[
  454. PolyElement, PolyElement, tuple[int, ...] | None, tuple[int, ...] | None
  455. ]:
  456. """Bring two Puiseux polynomials to a common monom and ns."""
  457. poly1, monom1, ns1 = self.poly, self.monom, self.ns
  458. poly2, monom2, ns2 = other.poly, other.monom, other.ns
  459. if monom1 == monom2 and ns1 == ns2:
  460. return poly1, poly2, monom1, ns1
  461. if ns1 == ns2:
  462. ns = ns1
  463. elif ns1 is not None and ns2 is not None:
  464. ns = tuple(lcm(n1, n2) for n1, n2 in zip(ns1, ns2))
  465. f1 = [n // n1 for n, n1 in zip(ns, ns1)]
  466. f2 = [n // n2 for n, n2 in zip(ns, ns2)]
  467. poly1 = poly1.inflate(f1)
  468. poly2 = poly2.inflate(f2)
  469. if monom1 is not None:
  470. monom1 = tuple(m * f for m, f in zip(monom1, f1))
  471. if monom2 is not None:
  472. monom2 = tuple(m * f for m, f in zip(monom2, f2))
  473. elif ns2 is not None:
  474. ns = ns2
  475. poly1 = poly1.inflate(ns)
  476. if monom1 is not None:
  477. monom1 = tuple(m * n for m, n in zip(monom1, ns))
  478. elif ns1 is not None:
  479. ns = ns1
  480. poly2 = poly2.inflate(ns)
  481. if monom2 is not None:
  482. monom2 = tuple(m * n for m, n in zip(monom2, ns))
  483. else:
  484. assert False
  485. if monom1 == monom2:
  486. monom = monom1
  487. elif monom1 is not None and monom2 is not None:
  488. monom = tuple(max(m1, m2) for m1, m2 in zip(monom1, monom2))
  489. poly1 = _mul_poly_monom(poly1, _div_monom(monom, monom1))
  490. poly2 = _mul_poly_monom(poly2, _div_monom(monom, monom2))
  491. elif monom2 is not None:
  492. monom = monom2
  493. poly1 = _mul_poly_monom(poly1, monom2)
  494. elif monom1 is not None:
  495. monom = monom1
  496. poly2 = _mul_poly_monom(poly2, monom1)
  497. else:
  498. assert False
  499. return poly1, poly2, monom, ns
  500. def __pos__(self) -> PuiseuxPoly:
  501. return self
  502. def __neg__(self) -> PuiseuxPoly:
  503. return self._new_raw(self.ring, -self.poly, self.monom, self.ns)
  504. def __add__(self, other: Any) -> PuiseuxPoly:
  505. if isinstance(other, PuiseuxPoly):
  506. if self.ring != other.ring:
  507. raise ValueError("Cannot add Puiseux polynomials from different rings")
  508. return self._add(other)
  509. domain = self.ring.domain
  510. if isinstance(other, int):
  511. return self._add_ground(domain.convert_from(QQ(other), QQ))
  512. elif domain.of_type(other):
  513. return self._add_ground(other)
  514. else:
  515. return NotImplemented
  516. def __radd__(self, other: Any) -> PuiseuxPoly:
  517. domain = self.ring.domain
  518. if isinstance(other, int):
  519. return self._add_ground(domain.convert_from(QQ(other), QQ))
  520. elif domain.of_type(other):
  521. return self._add_ground(other)
  522. else:
  523. return NotImplemented
  524. def __sub__(self, other: Any) -> PuiseuxPoly:
  525. if isinstance(other, PuiseuxPoly):
  526. if self.ring != other.ring:
  527. raise ValueError(
  528. "Cannot subtract Puiseux polynomials from different rings"
  529. )
  530. return self._sub(other)
  531. domain = self.ring.domain
  532. if isinstance(other, int):
  533. return self._sub_ground(domain.convert_from(QQ(other), QQ))
  534. elif domain.of_type(other):
  535. return self._sub_ground(other)
  536. else:
  537. return NotImplemented
  538. def __rsub__(self, other: Any) -> PuiseuxPoly:
  539. domain = self.ring.domain
  540. if isinstance(other, int):
  541. return self._rsub_ground(domain.convert_from(QQ(other), QQ))
  542. elif domain.of_type(other):
  543. return self._rsub_ground(other)
  544. else:
  545. return NotImplemented
  546. def __mul__(self, other: Any) -> PuiseuxPoly:
  547. if isinstance(other, PuiseuxPoly):
  548. if self.ring != other.ring:
  549. raise ValueError(
  550. "Cannot multiply Puiseux polynomials from different rings"
  551. )
  552. return self._mul(other)
  553. domain = self.ring.domain
  554. if isinstance(other, int):
  555. return self._mul_ground(domain.convert_from(QQ(other), QQ))
  556. elif domain.of_type(other):
  557. return self._mul_ground(other)
  558. else:
  559. return NotImplemented
  560. def __rmul__(self, other: Any) -> PuiseuxPoly:
  561. domain = self.ring.domain
  562. if isinstance(other, int):
  563. return self._mul_ground(domain.convert_from(QQ(other), QQ))
  564. elif domain.of_type(other):
  565. return self._mul_ground(other)
  566. else:
  567. return NotImplemented
  568. def __pow__(self, other: Any) -> PuiseuxPoly:
  569. if isinstance(other, int):
  570. if other >= 0:
  571. return self._pow_pint(other)
  572. else:
  573. return self._pow_nint(-other)
  574. elif QQ.of_type(other):
  575. return self._pow_rational(other)
  576. else:
  577. return NotImplemented
  578. def __truediv__(self, other: Any) -> PuiseuxPoly:
  579. if isinstance(other, PuiseuxPoly):
  580. if self.ring != other.ring:
  581. raise ValueError(
  582. "Cannot divide Puiseux polynomials from different rings"
  583. )
  584. return self._mul(other._inv())
  585. domain = self.ring.domain
  586. if isinstance(other, int):
  587. return self._mul_ground(domain.convert_from(QQ(1, other), QQ))
  588. elif domain.of_type(other):
  589. return self._div_ground(other)
  590. else:
  591. return NotImplemented
  592. def __rtruediv__(self, other: Any) -> PuiseuxPoly:
  593. if isinstance(other, int):
  594. return self._inv()._mul_ground(self.ring.domain.convert_from(QQ(other), QQ))
  595. elif self.ring.domain.of_type(other):
  596. return self._inv()._mul_ground(other)
  597. else:
  598. return NotImplemented
  599. def _add(self, other: PuiseuxPoly) -> PuiseuxPoly:
  600. poly1, poly2, monom, ns = self._unify(other)
  601. return self._new(self.ring, poly1 + poly2, monom, ns)
  602. def _add_ground(self, ground: Any) -> PuiseuxPoly:
  603. return self._add(self.ring.ground_new(ground))
  604. def _sub(self, other: PuiseuxPoly) -> PuiseuxPoly:
  605. poly1, poly2, monom, ns = self._unify(other)
  606. return self._new(self.ring, poly1 - poly2, monom, ns)
  607. def _sub_ground(self, ground: Any) -> PuiseuxPoly:
  608. return self._sub(self.ring.ground_new(ground))
  609. def _rsub_ground(self, ground: Any) -> PuiseuxPoly:
  610. return self.ring.ground_new(ground)._sub(self)
  611. def _mul(self, other: PuiseuxPoly) -> PuiseuxPoly:
  612. poly1, poly2, monom, ns = self._unify(other)
  613. if monom is not None:
  614. monom = tuple(2 * e for e in monom)
  615. return self._new(self.ring, poly1 * poly2, monom, ns)
  616. def _mul_ground(self, ground: Any) -> PuiseuxPoly:
  617. return self._new_raw(self.ring, self.poly * ground, self.monom, self.ns)
  618. def _div_ground(self, ground: Any) -> PuiseuxPoly:
  619. return self._new_raw(self.ring, self.poly / ground, self.monom, self.ns)
  620. def _pow_pint(self, n: int) -> PuiseuxPoly:
  621. assert n >= 0
  622. monom = self.monom
  623. if monom is not None:
  624. monom = tuple(m * n for m in monom)
  625. return self._new(self.ring, self.poly**n, monom, self.ns)
  626. def _pow_nint(self, n: int) -> PuiseuxPoly:
  627. return self._inv()._pow_pint(n)
  628. def _pow_rational(self, n: Any) -> PuiseuxPoly:
  629. if not self.is_term:
  630. raise ValueError("Only monomials can be raised to a rational power")
  631. [(monom, coeff)] = self.terms()
  632. domain = self.ring.domain
  633. if not domain.is_one(coeff):
  634. raise ValueError("Only monomials can be raised to a rational power")
  635. monom = tuple(m * n for m in monom)
  636. return self.ring.from_dict({monom: domain.one})
  637. def _inv(self) -> PuiseuxPoly:
  638. if not self.is_term:
  639. raise ValueError("Only terms can be inverted")
  640. [(monom, coeff)] = self.terms()
  641. domain = self.ring.domain
  642. if not domain.is_Field and not domain.is_one(coeff):
  643. raise ValueError("Cannot invert non-unit coefficient")
  644. monom = tuple(-m for m in monom)
  645. coeff = 1 / coeff
  646. return self.ring.from_dict({monom: coeff})
  647. def diff(self, x: PuiseuxPoly) -> PuiseuxPoly:
  648. """Differentiate a Puiseux polynomial with respect to a variable.
  649. >>> from sympy import QQ
  650. >>> from sympy.polys.puiseux import puiseux_ring
  651. >>> R, x, y = puiseux_ring('x, y', QQ)
  652. >>> p = 5*x**2 + 7*y**3
  653. >>> p.diff(x)
  654. 10*x
  655. >>> p.diff(y)
  656. 21*y**2
  657. """
  658. ring = self.ring
  659. i = ring.index(x)
  660. g = {}
  661. for expv, coeff in self.iterterms():
  662. n = expv[i]
  663. if n:
  664. e = list(expv)
  665. e[i] -= 1
  666. g[tuple(e)] = coeff * n
  667. return ring(g)