fields.py 21 KB

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