fields.py 21 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639
  1. """Sparse rational function fields. """
  2. from __future__ import annotations
  3. from functools import reduce
  4. from operator import add, mul, lt, le, gt, ge
  5. from sympy.core.expr import Expr
  6. from sympy.core.mod import Mod
  7. from sympy.core.numbers import Exp1
  8. from sympy.core.singleton import S
  9. from sympy.core.symbol import Symbol
  10. from sympy.core.sympify import CantSympify, sympify
  11. from sympy.functions.elementary.exponential import ExpBase
  12. from sympy.polys.domains.domain import Domain
  13. from sympy.polys.domains.domainelement import DomainElement
  14. from sympy.polys.domains.fractionfield import FractionField
  15. from sympy.polys.domains.polynomialring import PolynomialRing
  16. from sympy.polys.constructor import construct_domain
  17. from sympy.polys.orderings import lex, MonomialOrder
  18. from sympy.polys.polyerrors import CoercionFailed
  19. from sympy.polys.polyoptions import build_options
  20. from sympy.polys.polyutils import _parallel_dict_from_expr
  21. from sympy.polys.rings import PolyRing, PolyElement
  22. from sympy.printing.defaults import DefaultPrinting
  23. from sympy.utilities import public
  24. from sympy.utilities.iterables import is_sequence
  25. from sympy.utilities.magic import pollute
  26. @public
  27. def field(symbols, domain, order=lex):
  28. """Construct new rational function field returning (field, x1, ..., xn). """
  29. _field = FracField(symbols, domain, order)
  30. return (_field,) + _field.gens
  31. @public
  32. def xfield(symbols, domain, order=lex):
  33. """Construct new rational function field returning (field, (x1, ..., xn)). """
  34. _field = FracField(symbols, domain, order)
  35. return (_field, _field.gens)
  36. @public
  37. def vfield(symbols, domain, order=lex):
  38. """Construct new rational function field and inject generators into global namespace. """
  39. _field = FracField(symbols, domain, order)
  40. pollute([ sym.name for sym in _field.symbols ], _field.gens)
  41. return _field
  42. @public
  43. def sfield(exprs, *symbols, **options):
  44. """Construct a field deriving generators and domain
  45. from options and input expressions.
  46. Parameters
  47. ==========
  48. exprs : py:class:`~.Expr` or sequence of :py:class:`~.Expr` (sympifiable)
  49. symbols : sequence of :py:class:`~.Symbol`/:py:class:`~.Expr`
  50. options : keyword arguments understood by :py:class:`~.Options`
  51. Examples
  52. ========
  53. >>> from sympy import exp, log, symbols, sfield
  54. >>> x = symbols("x")
  55. >>> K, f = sfield((x*log(x) + 4*x**2)*exp(1/x + log(x)/3)/x**2)
  56. >>> K
  57. Rational function field in x, exp(1/x), log(x), x**(1/3) over ZZ with lex order
  58. >>> f
  59. (4*x**2*(exp(1/x)) + x*(exp(1/x))*(log(x)))/((x**(1/3))**5)
  60. """
  61. single = False
  62. if not is_sequence(exprs):
  63. exprs, single = [exprs], True
  64. exprs = list(map(sympify, exprs))
  65. opt = build_options(symbols, options)
  66. numdens = []
  67. for expr in exprs:
  68. numdens.extend(expr.as_numer_denom())
  69. reps, opt = _parallel_dict_from_expr(numdens, opt)
  70. if opt.domain is None:
  71. # NOTE: this is inefficient because construct_domain() automatically
  72. # performs conversion to the target domain. It shouldn't do this.
  73. coeffs = sum([list(rep.values()) for rep in reps], [])
  74. opt.domain, _ = construct_domain(coeffs, opt=opt)
  75. _field = FracField(opt.gens, opt.domain, opt.order)
  76. fracs = []
  77. for i in range(0, len(reps), 2):
  78. fracs.append(_field(tuple(reps[i:i+2])))
  79. if single:
  80. return (_field, fracs[0])
  81. else:
  82. return (_field, fracs)
  83. class FracField(DefaultPrinting):
  84. """Multivariate distributed rational function field. """
  85. ring: PolyRing
  86. gens: tuple[FracElement, ...]
  87. symbols: tuple[Expr, ...]
  88. ngens: int
  89. domain: Domain
  90. order: MonomialOrder
  91. def __new__(cls, symbols, domain, order=lex):
  92. ring = PolyRing(symbols, domain, order)
  93. symbols = ring.symbols
  94. ngens = ring.ngens
  95. domain = ring.domain
  96. order = ring.order
  97. _hash_tuple = (cls.__name__, symbols, ngens, domain, order)
  98. obj = object.__new__(cls)
  99. obj._hash_tuple = _hash_tuple
  100. obj._hash = hash(_hash_tuple)
  101. obj.ring = ring
  102. obj.symbols = symbols
  103. obj.ngens = ngens
  104. obj.domain = domain
  105. obj.order = order
  106. obj.dtype = FracElement(obj, ring.zero).raw_new
  107. obj.zero = obj.dtype(ring.zero)
  108. obj.one = obj.dtype(ring.one)
  109. obj.gens = obj._gens()
  110. for symbol, generator in zip(obj.symbols, obj.gens):
  111. if isinstance(symbol, Symbol):
  112. name = symbol.name
  113. if not hasattr(obj, name):
  114. setattr(obj, name, generator)
  115. return obj
  116. def _gens(self):
  117. """Return a list of polynomial generators. """
  118. return tuple([ self.dtype(gen) for gen in self.ring.gens ])
  119. def __getnewargs__(self):
  120. return (self.symbols, self.domain, self.order)
  121. def __hash__(self):
  122. return self._hash
  123. def index(self, gen):
  124. if self.is_element(gen):
  125. return self.ring.index(gen.to_poly())
  126. else:
  127. raise ValueError("expected a %s, got %s instead" % (self.dtype,gen))
  128. def __eq__(self, other):
  129. return isinstance(other, FracField) and \
  130. (self.symbols, self.ngens, self.domain, self.order) == \
  131. (other.symbols, other.ngens, other.domain, other.order)
  132. def __ne__(self, other):
  133. return not self == other
  134. def is_element(self, element):
  135. """True if ``element`` is an element of this field. False otherwise. """
  136. return isinstance(element, FracElement) and element.field == self
  137. def raw_new(self, numer, denom=None):
  138. return self.dtype(numer, denom)
  139. def new(self, numer, denom=None):
  140. if denom is None: denom = self.ring.one
  141. numer, denom = numer.cancel(denom)
  142. return self.raw_new(numer, denom)
  143. def domain_new(self, element):
  144. return self.domain.convert(element)
  145. def ground_new(self, element):
  146. try:
  147. return self.new(self.ring.ground_new(element))
  148. except CoercionFailed:
  149. domain = self.domain
  150. if not domain.is_Field and domain.has_assoc_Field:
  151. ring = self.ring
  152. ground_field = domain.get_field()
  153. element = ground_field.convert(element)
  154. numer = ring.ground_new(ground_field.numer(element))
  155. denom = ring.ground_new(ground_field.denom(element))
  156. return self.raw_new(numer, denom)
  157. else:
  158. raise
  159. def field_new(self, element):
  160. if isinstance(element, FracElement):
  161. if self == element.field:
  162. return element
  163. if isinstance(self.domain, FractionField) and \
  164. self.domain.field == element.field:
  165. return self.ground_new(element)
  166. elif isinstance(self.domain, PolynomialRing) and \
  167. self.domain.ring.to_field() == element.field:
  168. return self.ground_new(element)
  169. else:
  170. raise NotImplementedError("conversion")
  171. elif isinstance(element, PolyElement):
  172. denom, numer = element.clear_denoms()
  173. if isinstance(self.domain, PolynomialRing) and \
  174. numer.ring == self.domain.ring:
  175. numer = self.ring.ground_new(numer)
  176. elif isinstance(self.domain, FractionField) and \
  177. numer.ring == self.domain.field.to_ring():
  178. numer = self.ring.ground_new(numer)
  179. else:
  180. numer = numer.set_ring(self.ring)
  181. denom = self.ring.ground_new(denom)
  182. return self.raw_new(numer, denom)
  183. elif isinstance(element, tuple) and len(element) == 2:
  184. numer, denom = list(map(self.ring.ring_new, element))
  185. return self.new(numer, denom)
  186. elif isinstance(element, str):
  187. raise NotImplementedError("parsing")
  188. elif isinstance(element, Expr):
  189. return self.from_expr(element)
  190. else:
  191. return self.ground_new(element)
  192. __call__ = field_new
  193. def _rebuild_expr(self, expr, mapping):
  194. domain = self.domain
  195. powers = tuple((gen, gen.as_base_exp()) for gen in mapping.keys()
  196. if gen.is_Pow or isinstance(gen, ExpBase))
  197. def _rebuild(expr):
  198. generator = mapping.get(expr)
  199. if generator is not None:
  200. return generator
  201. elif expr.is_Add:
  202. return reduce(add, list(map(_rebuild, expr.args)))
  203. elif expr.is_Mul:
  204. return reduce(mul, list(map(_rebuild, expr.args)))
  205. elif expr.is_Pow or isinstance(expr, (ExpBase, Exp1)):
  206. b, e = expr.as_base_exp()
  207. # look for bg**eg whose integer power may be b**e
  208. for gen, (bg, eg) in powers:
  209. if bg == b and Mod(e, eg) == 0:
  210. return mapping.get(gen)**int(e/eg)
  211. if e.is_Integer and e is not S.One:
  212. return _rebuild(b)**int(e)
  213. elif mapping.get(1/expr) is not None:
  214. return 1/mapping.get(1/expr)
  215. try:
  216. return domain.convert(expr)
  217. except CoercionFailed:
  218. if not domain.is_Field and domain.has_assoc_Field:
  219. return domain.get_field().convert(expr)
  220. else:
  221. raise
  222. return _rebuild(expr)
  223. def from_expr(self, expr):
  224. mapping = dict(list(zip(self.symbols, self.gens)))
  225. try:
  226. frac = self._rebuild_expr(sympify(expr), mapping)
  227. except CoercionFailed:
  228. raise ValueError("expected an expression convertible to a rational function in %s, got %s" % (self, expr))
  229. else:
  230. return self.field_new(frac)
  231. def to_domain(self):
  232. return FractionField(self)
  233. def to_ring(self):
  234. return PolyRing(self.symbols, self.domain, self.order)
  235. class FracElement(DomainElement, DefaultPrinting, CantSympify):
  236. """Element of multivariate distributed rational function field. """
  237. def __init__(self, field, numer, denom=None):
  238. if denom is None:
  239. denom = field.ring.one
  240. elif not denom:
  241. raise ZeroDivisionError("zero denominator")
  242. self.field = field
  243. self.numer = numer
  244. self.denom = denom
  245. def raw_new(f, numer, denom=None):
  246. return f.__class__(f.field, numer, denom)
  247. def new(f, numer, denom):
  248. return f.raw_new(*numer.cancel(denom))
  249. def to_poly(f):
  250. if f.denom != 1:
  251. raise ValueError("f.denom should be 1")
  252. return f.numer
  253. def parent(self):
  254. return self.field.to_domain()
  255. def __getnewargs__(self):
  256. return (self.field, self.numer, self.denom)
  257. _hash = None
  258. def __hash__(self):
  259. _hash = self._hash
  260. if _hash is None:
  261. self._hash = _hash = hash((self.field, self.numer, self.denom))
  262. return _hash
  263. def copy(self):
  264. return self.raw_new(self.numer.copy(), self.denom.copy())
  265. def set_field(self, new_field):
  266. if self.field == new_field:
  267. return self
  268. else:
  269. new_ring = new_field.ring
  270. numer = self.numer.set_ring(new_ring)
  271. denom = self.denom.set_ring(new_ring)
  272. return new_field.new(numer, denom)
  273. def as_expr(self, *symbols):
  274. return self.numer.as_expr(*symbols)/self.denom.as_expr(*symbols)
  275. def __eq__(f, g):
  276. if isinstance(g, FracElement) and f.field == g.field:
  277. return f.numer == g.numer and f.denom == g.denom
  278. else:
  279. return f.numer == g and f.denom == f.field.ring.one
  280. def __ne__(f, g):
  281. return not f == g
  282. def __bool__(f):
  283. return bool(f.numer)
  284. def sort_key(self):
  285. return (self.denom.sort_key(), self.numer.sort_key())
  286. def _cmp(f1, f2, op):
  287. if f1.field.is_element(f2):
  288. return op(f1.sort_key(), f2.sort_key())
  289. else:
  290. return NotImplemented
  291. def __lt__(f1, f2):
  292. return f1._cmp(f2, lt)
  293. def __le__(f1, f2):
  294. return f1._cmp(f2, le)
  295. def __gt__(f1, f2):
  296. return f1._cmp(f2, gt)
  297. def __ge__(f1, f2):
  298. return f1._cmp(f2, ge)
  299. def __pos__(f):
  300. """Negate all coefficients in ``f``. """
  301. return f.raw_new(f.numer, f.denom)
  302. def __neg__(f):
  303. """Negate all coefficients in ``f``. """
  304. return f.raw_new(-f.numer, f.denom)
  305. def _extract_ground(self, element):
  306. domain = self.field.domain
  307. try:
  308. element = domain.convert(element)
  309. except CoercionFailed:
  310. if not domain.is_Field and domain.has_assoc_Field:
  311. ground_field = domain.get_field()
  312. try:
  313. element = ground_field.convert(element)
  314. except CoercionFailed:
  315. pass
  316. else:
  317. return -1, ground_field.numer(element), ground_field.denom(element)
  318. return 0, None, None
  319. else:
  320. return 1, element, None
  321. def __add__(f, g):
  322. """Add rational functions ``f`` and ``g``. """
  323. field = f.field
  324. if not g:
  325. return f
  326. elif not f:
  327. return g
  328. elif field.is_element(g):
  329. if f.denom == g.denom:
  330. return f.new(f.numer + g.numer, f.denom)
  331. else:
  332. return f.new(f.numer*g.denom + f.denom*g.numer, f.denom*g.denom)
  333. elif field.ring.is_element(g):
  334. return f.new(f.numer + f.denom*g, f.denom)
  335. else:
  336. if isinstance(g, FracElement):
  337. if isinstance(field.domain, FractionField) and field.domain.field == g.field:
  338. pass
  339. elif isinstance(g.field.domain, FractionField) and g.field.domain.field == field:
  340. return g.__radd__(f)
  341. else:
  342. return NotImplemented
  343. elif isinstance(g, PolyElement):
  344. if isinstance(field.domain, PolynomialRing) and field.domain.ring == g.ring:
  345. pass
  346. else:
  347. return g.__radd__(f)
  348. return f.__radd__(g)
  349. def __radd__(f, c):
  350. if f.field.ring.is_element(c):
  351. return f.new(f.numer + f.denom*c, f.denom)
  352. op, g_numer, g_denom = f._extract_ground(c)
  353. if op == 1:
  354. return f.new(f.numer + f.denom*g_numer, f.denom)
  355. elif not op:
  356. return NotImplemented
  357. else:
  358. return f.new(f.numer*g_denom + f.denom*g_numer, f.denom*g_denom)
  359. def __sub__(f, g):
  360. """Subtract rational functions ``f`` and ``g``. """
  361. field = f.field
  362. if not g:
  363. return f
  364. elif not f:
  365. return -g
  366. elif field.is_element(g):
  367. if f.denom == g.denom:
  368. return f.new(f.numer - g.numer, f.denom)
  369. else:
  370. return f.new(f.numer*g.denom - f.denom*g.numer, f.denom*g.denom)
  371. elif field.ring.is_element(g):
  372. return f.new(f.numer - f.denom*g, f.denom)
  373. else:
  374. if isinstance(g, FracElement):
  375. if isinstance(field.domain, FractionField) and field.domain.field == g.field:
  376. pass
  377. elif isinstance(g.field.domain, FractionField) and g.field.domain.field == field:
  378. return g.__rsub__(f)
  379. else:
  380. return NotImplemented
  381. elif isinstance(g, PolyElement):
  382. if isinstance(field.domain, PolynomialRing) and field.domain.ring == g.ring:
  383. pass
  384. else:
  385. return g.__rsub__(f)
  386. op, g_numer, g_denom = f._extract_ground(g)
  387. if op == 1:
  388. return f.new(f.numer - f.denom*g_numer, f.denom)
  389. elif not op:
  390. return NotImplemented
  391. else:
  392. return f.new(f.numer*g_denom - f.denom*g_numer, f.denom*g_denom)
  393. def __rsub__(f, c):
  394. if f.field.ring.is_element(c):
  395. return f.new(-f.numer + f.denom*c, f.denom)
  396. op, g_numer, g_denom = f._extract_ground(c)
  397. if op == 1:
  398. return f.new(-f.numer + f.denom*g_numer, f.denom)
  399. elif not op:
  400. return NotImplemented
  401. else:
  402. return f.new(-f.numer*g_denom + f.denom*g_numer, f.denom*g_denom)
  403. def __mul__(f, g):
  404. """Multiply rational functions ``f`` and ``g``. """
  405. field = f.field
  406. if not f or not g:
  407. return field.zero
  408. elif field.is_element(g):
  409. return f.new(f.numer*g.numer, f.denom*g.denom)
  410. elif field.ring.is_element(g):
  411. return f.new(f.numer*g, f.denom)
  412. else:
  413. if isinstance(g, FracElement):
  414. if isinstance(field.domain, FractionField) and field.domain.field == g.field:
  415. pass
  416. elif isinstance(g.field.domain, FractionField) and g.field.domain.field == field:
  417. return g.__rmul__(f)
  418. else:
  419. return NotImplemented
  420. elif isinstance(g, PolyElement):
  421. if isinstance(field.domain, PolynomialRing) and field.domain.ring == g.ring:
  422. pass
  423. else:
  424. return g.__rmul__(f)
  425. return f.__rmul__(g)
  426. def __rmul__(f, c):
  427. if f.field.ring.is_element(c):
  428. return f.new(f.numer*c, f.denom)
  429. op, g_numer, g_denom = f._extract_ground(c)
  430. if op == 1:
  431. return f.new(f.numer*g_numer, f.denom)
  432. elif not op:
  433. return NotImplemented
  434. else:
  435. return f.new(f.numer*g_numer, f.denom*g_denom)
  436. def __truediv__(f, g):
  437. """Computes quotient of fractions ``f`` and ``g``. """
  438. field = f.field
  439. if not g:
  440. raise ZeroDivisionError
  441. elif field.is_element(g):
  442. return f.new(f.numer*g.denom, f.denom*g.numer)
  443. elif field.ring.is_element(g):
  444. return f.new(f.numer, f.denom*g)
  445. else:
  446. if isinstance(g, FracElement):
  447. if isinstance(field.domain, FractionField) and field.domain.field == g.field:
  448. pass
  449. elif isinstance(g.field.domain, FractionField) and g.field.domain.field == field:
  450. return g.__rtruediv__(f)
  451. else:
  452. return NotImplemented
  453. elif isinstance(g, PolyElement):
  454. if isinstance(field.domain, PolynomialRing) and field.domain.ring == g.ring:
  455. pass
  456. else:
  457. return g.__rtruediv__(f)
  458. op, g_numer, g_denom = f._extract_ground(g)
  459. if op == 1:
  460. return f.new(f.numer, f.denom*g_numer)
  461. elif not op:
  462. return NotImplemented
  463. else:
  464. return f.new(f.numer*g_denom, f.denom*g_numer)
  465. def __rtruediv__(f, c):
  466. if not f:
  467. raise ZeroDivisionError
  468. elif f.field.ring.is_element(c):
  469. return f.new(f.denom*c, f.numer)
  470. op, g_numer, g_denom = f._extract_ground(c)
  471. if op == 1:
  472. return f.new(f.denom*g_numer, f.numer)
  473. elif not op:
  474. return NotImplemented
  475. else:
  476. return f.new(f.denom*g_numer, f.numer*g_denom)
  477. def __pow__(f, n):
  478. """Raise ``f`` to a non-negative power ``n``. """
  479. if n >= 0:
  480. return f.raw_new(f.numer**n, f.denom**n)
  481. elif not f:
  482. raise ZeroDivisionError
  483. else:
  484. return f.raw_new(f.denom**-n, f.numer**-n)
  485. def diff(f, x):
  486. """Computes partial derivative in ``x``.
  487. Examples
  488. ========
  489. >>> from sympy.polys.fields import field
  490. >>> from sympy.polys.domains import ZZ
  491. >>> _, x, y, z = field("x,y,z", ZZ)
  492. >>> ((x**2 + y)/(z + 1)).diff(x)
  493. 2*x/(z + 1)
  494. """
  495. x = x.to_poly()
  496. return f.new(f.numer.diff(x)*f.denom - f.numer*f.denom.diff(x), f.denom**2)
  497. def __call__(f, *values):
  498. if 0 < len(values) <= f.field.ngens:
  499. return f.evaluate(list(zip(f.field.gens, values)))
  500. else:
  501. raise ValueError("expected at least 1 and at most %s values, got %s" % (f.field.ngens, len(values)))
  502. def evaluate(f, x, a=None):
  503. if isinstance(x, list) and a is None:
  504. x = [ (X.to_poly(), a) for X, a in x ]
  505. numer, denom = f.numer.evaluate(x), f.denom.evaluate(x)
  506. else:
  507. x = x.to_poly()
  508. numer, denom = f.numer.evaluate(x, a), f.denom.evaluate(x, a)
  509. field = numer.ring.to_field()
  510. return field.new(numer, denom)
  511. def subs(f, x, a=None):
  512. if isinstance(x, list) and a is None:
  513. x = [ (X.to_poly(), a) for X, a in x ]
  514. numer, denom = f.numer.subs(x), f.denom.subs(x)
  515. else:
  516. x = x.to_poly()
  517. numer, denom = f.numer.subs(x, a), f.denom.subs(x, a)
  518. return f.new(numer, denom)
  519. def compose(f, x, a=None):
  520. raise NotImplementedError