polyutils.py 17 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584
  1. """Useful utilities for higher level polynomial classes. """
  2. from __future__ import annotations
  3. from sympy.external.gmpy import GROUND_TYPES
  4. from sympy.core import (S, Add, Mul, Pow, Eq, Expr,
  5. expand_mul, expand_multinomial)
  6. from sympy.core.exprtools import decompose_power, decompose_power_rat
  7. from sympy.core.numbers import _illegal
  8. from sympy.polys.polyerrors import PolynomialError, GeneratorsError
  9. from sympy.polys.polyoptions import build_options
  10. import re
  11. _gens_order = {
  12. 'a': 301, 'b': 302, 'c': 303, 'd': 304,
  13. 'e': 305, 'f': 306, 'g': 307, 'h': 308,
  14. 'i': 309, 'j': 310, 'k': 311, 'l': 312,
  15. 'm': 313, 'n': 314, 'o': 315, 'p': 216,
  16. 'q': 217, 'r': 218, 's': 219, 't': 220,
  17. 'u': 221, 'v': 222, 'w': 223, 'x': 124,
  18. 'y': 125, 'z': 126,
  19. }
  20. _max_order = 1000
  21. _re_gen = re.compile(r"^(.*?)(\d*)$", re.MULTILINE)
  22. def _nsort(roots, separated=False):
  23. """Sort the numerical roots putting the real roots first, then sorting
  24. according to real and imaginary parts. If ``separated`` is True, then
  25. the real and imaginary roots will be returned in two lists, respectively.
  26. This routine tries to avoid issue 6137 by separating the roots into real
  27. and imaginary parts before evaluation. In addition, the sorting will raise
  28. an error if any computation cannot be done with precision.
  29. """
  30. if not all(r.is_number for r in roots):
  31. raise NotImplementedError
  32. if not len(roots):
  33. return [] if not separated else ([], [])
  34. # see issue 6137:
  35. # get the real part of the evaluated real and imaginary parts of each root
  36. key = [[i.n(2).as_real_imag()[0] for i in r.as_real_imag()] for r in roots]
  37. # make sure the parts were computed with precision
  38. if len(roots) > 1 and any(i._prec == 1 for k in key for i in k):
  39. raise NotImplementedError("could not compute root with precision")
  40. # insert a key to indicate if the root has an imaginary part
  41. key = [(1 if i else 0, r, i) for r, i in key]
  42. key = sorted(zip(key, roots))
  43. # return the real and imaginary roots separately if desired
  44. if separated:
  45. r = []
  46. i = []
  47. for (im, _, _), v in key:
  48. if im:
  49. i.append(v)
  50. else:
  51. r.append(v)
  52. return r, i
  53. _, roots = zip(*key)
  54. return list(roots)
  55. def _sort_gens(gens, **args):
  56. """Sort generators in a reasonably intelligent way. """
  57. opt = build_options(args)
  58. gens_order, wrt = {}, None
  59. if opt is not None:
  60. gens_order, wrt = {}, opt.wrt
  61. for i, gen in enumerate(opt.sort):
  62. gens_order[gen] = i + 1
  63. def order_key(gen):
  64. gen = str(gen)
  65. if wrt is not None:
  66. try:
  67. return (-len(wrt) + wrt.index(gen), gen, 0)
  68. except ValueError:
  69. pass
  70. name, index = _re_gen.match(gen).groups()
  71. if index:
  72. index = int(index)
  73. else:
  74. index = 0
  75. try:
  76. return ( gens_order[name], name, index)
  77. except KeyError:
  78. pass
  79. try:
  80. return (_gens_order[name], name, index)
  81. except KeyError:
  82. pass
  83. return (_max_order, name, index)
  84. try:
  85. gens = sorted(gens, key=order_key)
  86. except TypeError: # pragma: no cover
  87. pass
  88. return tuple(gens)
  89. def _unify_gens(f_gens, g_gens):
  90. """Unify generators in a reasonably intelligent way. """
  91. f_gens = list(f_gens)
  92. g_gens = list(g_gens)
  93. if f_gens == g_gens:
  94. return tuple(f_gens)
  95. gens, common, k = [], [], 0
  96. for gen in f_gens:
  97. if gen in g_gens:
  98. common.append(gen)
  99. for i, gen in enumerate(g_gens):
  100. if gen in common:
  101. g_gens[i], k = common[k], k + 1
  102. for gen in common:
  103. i = f_gens.index(gen)
  104. gens.extend(f_gens[:i])
  105. f_gens = f_gens[i + 1:]
  106. i = g_gens.index(gen)
  107. gens.extend(g_gens[:i])
  108. g_gens = g_gens[i + 1:]
  109. gens.append(gen)
  110. gens.extend(f_gens)
  111. gens.extend(g_gens)
  112. return tuple(gens)
  113. def _analyze_gens(gens):
  114. """Support for passing generators as `*gens` and `[gens]`. """
  115. if len(gens) == 1 and hasattr(gens[0], '__iter__'):
  116. return tuple(gens[0])
  117. else:
  118. return tuple(gens)
  119. def _sort_factors(factors, **args):
  120. """Sort low-level factors in increasing 'complexity' order. """
  121. # XXX: GF(p) does not support comparisons so we need a key function to sort
  122. # the factors if python-flint is being used. A better solution might be to
  123. # add a sort key method to each domain.
  124. def order_key(factor):
  125. if isinstance(factor, _GF_types):
  126. return int(factor)
  127. elif isinstance(factor, list):
  128. return [order_key(f) for f in factor]
  129. else:
  130. return factor
  131. def order_if_multiple_key(factor):
  132. (f, n) = factor
  133. return (len(f), n, order_key(f))
  134. def order_no_multiple_key(f):
  135. return (len(f), order_key(f))
  136. if args.get('multiple', True):
  137. return sorted(factors, key=order_if_multiple_key)
  138. else:
  139. return sorted(factors, key=order_no_multiple_key)
  140. illegal_types = [type(obj) for obj in _illegal]
  141. finf = [float(i) for i in _illegal[1:3]]
  142. def _not_a_coeff(expr):
  143. """Do not treat NaN and infinities as valid polynomial coefficients. """
  144. if type(expr) in illegal_types or expr in finf:
  145. return True
  146. if isinstance(expr, float) and float(expr) != expr:
  147. return True # nan
  148. return # could be
  149. def _parallel_dict_from_expr_if_gens(exprs, opt):
  150. """Transform expressions into a multinomial form given generators. """
  151. k, indices = len(opt.gens), {}
  152. for i, g in enumerate(opt.gens):
  153. indices[g] = i
  154. polys = []
  155. for expr in exprs:
  156. poly = {}
  157. if expr.is_Equality:
  158. expr = expr.lhs - expr.rhs
  159. for term in Add.make_args(expr):
  160. coeff, monom = [], [0]*k
  161. for factor in Mul.make_args(term):
  162. if not _not_a_coeff(factor) and factor.is_Number:
  163. coeff.append(factor)
  164. else:
  165. try:
  166. if opt.series is False:
  167. base, exp = decompose_power(factor)
  168. if exp < 0:
  169. exp, base = -exp, Pow(base, -S.One)
  170. else:
  171. base, exp = decompose_power_rat(factor)
  172. monom[indices[base]] = exp
  173. except KeyError:
  174. if not factor.has_free(*opt.gens):
  175. coeff.append(factor)
  176. else:
  177. raise PolynomialError("%s contains an element of "
  178. "the set of generators." % factor)
  179. monom = tuple(monom)
  180. if monom in poly:
  181. poly[monom] += Mul(*coeff)
  182. else:
  183. poly[monom] = Mul(*coeff)
  184. polys.append(poly)
  185. return polys, opt.gens
  186. def _parallel_dict_from_expr_no_gens(exprs, opt):
  187. """Transform expressions into a multinomial form and figure out generators. """
  188. if opt.domain is not None:
  189. def _is_coeff(factor):
  190. return factor in opt.domain
  191. elif opt.extension is True:
  192. def _is_coeff(factor):
  193. return factor.is_algebraic
  194. elif opt.greedy is not False:
  195. def _is_coeff(factor):
  196. return factor is S.ImaginaryUnit
  197. else:
  198. def _is_coeff(factor):
  199. return factor.is_number
  200. gens, reprs = set(), []
  201. for expr in exprs:
  202. terms = []
  203. if expr.is_Equality:
  204. expr = expr.lhs - expr.rhs
  205. for term in Add.make_args(expr):
  206. coeff, elements = [], {}
  207. for factor in Mul.make_args(term):
  208. if not _not_a_coeff(factor) and (factor.is_Number or _is_coeff(factor)):
  209. coeff.append(factor)
  210. else:
  211. if opt.series is False:
  212. base, exp = decompose_power(factor)
  213. if exp < 0:
  214. exp, base = -exp, Pow(base, -S.One)
  215. else:
  216. base, exp = decompose_power_rat(factor)
  217. elements[base] = elements.setdefault(base, 0) + exp
  218. gens.add(base)
  219. terms.append((coeff, elements))
  220. reprs.append(terms)
  221. gens = _sort_gens(gens, opt=opt)
  222. k, indices = len(gens), {}
  223. for i, g in enumerate(gens):
  224. indices[g] = i
  225. polys = []
  226. for terms in reprs:
  227. poly = {}
  228. for coeff, term in terms:
  229. monom = [0]*k
  230. for base, exp in term.items():
  231. monom[indices[base]] = exp
  232. monom = tuple(monom)
  233. if monom in poly:
  234. poly[monom] += Mul(*coeff)
  235. else:
  236. poly[monom] = Mul(*coeff)
  237. polys.append(poly)
  238. return polys, tuple(gens)
  239. def _dict_from_expr_if_gens(expr, opt):
  240. """Transform an expression into a multinomial form given generators. """
  241. (poly,), gens = _parallel_dict_from_expr_if_gens((expr,), opt)
  242. return poly, gens
  243. def _dict_from_expr_no_gens(expr, opt):
  244. """Transform an expression into a multinomial form and figure out generators. """
  245. (poly,), gens = _parallel_dict_from_expr_no_gens((expr,), opt)
  246. return poly, gens
  247. def parallel_dict_from_expr(exprs, **args):
  248. """Transform expressions into a multinomial form. """
  249. reps, opt = _parallel_dict_from_expr(exprs, build_options(args))
  250. return reps, opt.gens
  251. def _parallel_dict_from_expr(exprs, opt):
  252. """Transform expressions into a multinomial form. """
  253. if opt.expand is not False:
  254. exprs = [ expr.expand() for expr in exprs ]
  255. if any(expr.is_commutative is False for expr in exprs):
  256. raise PolynomialError('non-commutative expressions are not supported')
  257. if opt.gens:
  258. reps, gens = _parallel_dict_from_expr_if_gens(exprs, opt)
  259. else:
  260. reps, gens = _parallel_dict_from_expr_no_gens(exprs, opt)
  261. return reps, opt.clone({'gens': gens})
  262. def dict_from_expr(expr, **args):
  263. """Transform an expression into a multinomial form. """
  264. rep, opt = _dict_from_expr(expr, build_options(args))
  265. return rep, opt.gens
  266. def _dict_from_expr(expr, opt):
  267. """Transform an expression into a multinomial form. """
  268. if expr.is_commutative is False:
  269. raise PolynomialError('non-commutative expressions are not supported')
  270. def _is_expandable_pow(expr):
  271. return (expr.is_Pow and expr.exp.is_positive and expr.exp.is_Integer
  272. and expr.base.is_Add)
  273. if opt.expand is not False:
  274. if not isinstance(expr, (Expr, Eq)):
  275. raise PolynomialError('expression must be of type Expr')
  276. expr = expr.expand()
  277. # TODO: Integrate this into expand() itself
  278. while any(_is_expandable_pow(i) or i.is_Mul and
  279. any(_is_expandable_pow(j) for j in i.args) for i in
  280. Add.make_args(expr)):
  281. expr = expand_multinomial(expr)
  282. while any(i.is_Mul and any(j.is_Add for j in i.args) for i in Add.make_args(expr)):
  283. expr = expand_mul(expr)
  284. if opt.gens:
  285. rep, gens = _dict_from_expr_if_gens(expr, opt)
  286. else:
  287. rep, gens = _dict_from_expr_no_gens(expr, opt)
  288. return rep, opt.clone({'gens': gens})
  289. def expr_from_dict(rep, *gens):
  290. """Convert a multinomial form into an expression. """
  291. result = []
  292. for monom, coeff in rep.items():
  293. term = [coeff]
  294. for g, m in zip(gens, monom):
  295. if m:
  296. term.append(Pow(g, m))
  297. result.append(Mul(*term))
  298. return Add(*result)
  299. parallel_dict_from_basic = parallel_dict_from_expr
  300. dict_from_basic = dict_from_expr
  301. basic_from_dict = expr_from_dict
  302. def _dict_reorder(rep, gens, new_gens):
  303. """Reorder levels using dict representation. """
  304. gens = list(gens)
  305. monoms = rep.keys()
  306. coeffs = rep.values()
  307. new_monoms = [ [] for _ in range(len(rep)) ]
  308. used_indices = set()
  309. for gen in new_gens:
  310. try:
  311. j = gens.index(gen)
  312. used_indices.add(j)
  313. for M, new_M in zip(monoms, new_monoms):
  314. new_M.append(M[j])
  315. except ValueError:
  316. for new_M in new_monoms:
  317. new_M.append(0)
  318. for i, _ in enumerate(gens):
  319. if i not in used_indices:
  320. for monom in monoms:
  321. if monom[i]:
  322. raise GeneratorsError("unable to drop generators")
  323. return map(tuple, new_monoms), coeffs
  324. class PicklableWithSlots:
  325. """
  326. Mixin class that allows to pickle objects with ``__slots__``.
  327. Examples
  328. ========
  329. First define a class that mixes :class:`PicklableWithSlots` in::
  330. >>> from sympy.polys.polyutils import PicklableWithSlots
  331. >>> class Some(PicklableWithSlots):
  332. ... __slots__ = ('foo', 'bar')
  333. ...
  334. ... def __init__(self, foo, bar):
  335. ... self.foo = foo
  336. ... self.bar = bar
  337. To make :mod:`pickle` happy in doctest we have to use these hacks::
  338. >>> import builtins
  339. >>> builtins.Some = Some
  340. >>> from sympy.polys import polyutils
  341. >>> polyutils.Some = Some
  342. Next lets see if we can create an instance, pickle it and unpickle::
  343. >>> some = Some('abc', 10)
  344. >>> some.foo, some.bar
  345. ('abc', 10)
  346. >>> from pickle import dumps, loads
  347. >>> some2 = loads(dumps(some))
  348. >>> some2.foo, some2.bar
  349. ('abc', 10)
  350. """
  351. __slots__ = ()
  352. def __getstate__(self, cls=None):
  353. if cls is None:
  354. # This is the case for the instance that gets pickled
  355. cls = self.__class__
  356. d = {}
  357. # Get all data that should be stored from super classes
  358. for c in cls.__bases__:
  359. # XXX: Python 3.11 defines object.__getstate__ and it does not
  360. # accept any arguments so we need to make sure not to call it with
  361. # an argument here. To be compatible with Python < 3.11 we need to
  362. # be careful not to assume that c or object has a __getstate__
  363. # method though.
  364. getstate = getattr(c, "__getstate__", None)
  365. objstate = getattr(object, "__getstate__", None)
  366. if getstate is not None and getstate is not objstate:
  367. d.update(getstate(self, c))
  368. # Get all information that should be stored from cls and return the dict
  369. for name in cls.__slots__:
  370. if hasattr(self, name):
  371. d[name] = getattr(self, name)
  372. return d
  373. def __setstate__(self, d):
  374. # All values that were pickled are now assigned to a fresh instance
  375. for name, value in d.items():
  376. setattr(self, name, value)
  377. class IntegerPowerable:
  378. r"""
  379. Mixin class for classes that define a `__mul__` method, and want to be
  380. raised to integer powers in the natural way that follows. Implements
  381. powering via binary expansion, for efficiency.
  382. By default, only integer powers $\geq 2$ are supported. To support the
  383. first, zeroth, or negative powers, override the corresponding methods,
  384. `_first_power`, `_zeroth_power`, `_negative_power`, below.
  385. """
  386. def __pow__(self, e, modulo=None):
  387. if e < 2:
  388. try:
  389. if e == 1:
  390. return self._first_power()
  391. elif e == 0:
  392. return self._zeroth_power()
  393. else:
  394. return self._negative_power(e, modulo=modulo)
  395. except NotImplementedError:
  396. return NotImplemented
  397. else:
  398. bits = [int(d) for d in reversed(bin(e)[2:])]
  399. n = len(bits)
  400. p = self
  401. first = True
  402. for i in range(n):
  403. if bits[i]:
  404. if first:
  405. r = p
  406. first = False
  407. else:
  408. r *= p
  409. if modulo is not None:
  410. r %= modulo
  411. if i < n - 1:
  412. p *= p
  413. if modulo is not None:
  414. p %= modulo
  415. return r
  416. def _negative_power(self, e, modulo=None):
  417. """
  418. Compute inverse of self, then raise that to the abs(e) power.
  419. For example, if the class has an `inv()` method,
  420. return self.inv() ** abs(e) % modulo
  421. """
  422. raise NotImplementedError
  423. def _zeroth_power(self):
  424. """Return unity element of algebraic struct to which self belongs."""
  425. raise NotImplementedError
  426. def _first_power(self):
  427. """Return a copy of self."""
  428. raise NotImplementedError
  429. _GF_types: tuple[type, ...]
  430. if GROUND_TYPES == 'flint':
  431. import flint
  432. _GF_types = (flint.nmod, flint.fmpz_mod)
  433. else:
  434. from sympy.polys.domains.modularinteger import ModularInteger
  435. flint = None
  436. _GF_types = (ModularInteger,)