monomials.py 18 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628
  1. """Tools and arithmetics for monomials of distributed polynomials. """
  2. from itertools import combinations_with_replacement, product
  3. from textwrap import dedent
  4. from sympy.core.cache import cacheit
  5. from sympy.core import Mul, S, Tuple, sympify
  6. from sympy.polys.polyerrors import ExactQuotientFailed
  7. from sympy.polys.polyutils import PicklableWithSlots, dict_from_expr
  8. from sympy.utilities import public
  9. from sympy.utilities.iterables import is_sequence, iterable
  10. @public
  11. def itermonomials(variables, max_degrees, min_degrees=None):
  12. r"""
  13. ``max_degrees`` and ``min_degrees`` are either both integers or both lists.
  14. Unless otherwise specified, ``min_degrees`` is either ``0`` or
  15. ``[0, ..., 0]``.
  16. A generator of all monomials ``monom`` is returned, such that
  17. either
  18. ``min_degree <= total_degree(monom) <= max_degree``,
  19. or
  20. ``min_degrees[i] <= degree_list(monom)[i] <= max_degrees[i]``,
  21. for all ``i``.
  22. Case I. ``max_degrees`` and ``min_degrees`` are both integers
  23. =============================================================
  24. Given a set of variables $V$ and a min_degree $N$ and a max_degree $M$
  25. generate a set of monomials of degree less than or equal to $N$ and greater
  26. than or equal to $M$. The total number of monomials in commutative
  27. variables is huge and is given by the following formula if $M = 0$:
  28. .. math::
  29. \frac{(\#V + N)!}{\#V! N!}
  30. For example if we would like to generate a dense polynomial of
  31. a total degree $N = 50$ and $M = 0$, which is the worst case, in 5
  32. variables, assuming that exponents and all of coefficients are 32-bit long
  33. and stored in an array we would need almost 80 GiB of memory! Fortunately
  34. most polynomials, that we will encounter, are sparse.
  35. Consider monomials in commutative variables $x$ and $y$
  36. and non-commutative variables $a$ and $b$::
  37. >>> from sympy import symbols
  38. >>> from sympy.polys.monomials import itermonomials
  39. >>> from sympy.polys.orderings import monomial_key
  40. >>> from sympy.abc import x, y
  41. >>> sorted(itermonomials([x, y], 2), key=monomial_key('grlex', [y, x]))
  42. [1, x, y, x**2, x*y, y**2]
  43. >>> sorted(itermonomials([x, y], 3), key=monomial_key('grlex', [y, x]))
  44. [1, x, y, x**2, x*y, y**2, x**3, x**2*y, x*y**2, y**3]
  45. >>> a, b = symbols('a, b', commutative=False)
  46. >>> set(itermonomials([a, b, x], 2))
  47. {1, a, a**2, b, b**2, x, x**2, a*b, b*a, x*a, x*b}
  48. >>> sorted(itermonomials([x, y], 2, 1), key=monomial_key('grlex', [y, x]))
  49. [x, y, x**2, x*y, y**2]
  50. Case II. ``max_degrees`` and ``min_degrees`` are both lists
  51. ===========================================================
  52. If ``max_degrees = [d_1, ..., d_n]`` and
  53. ``min_degrees = [e_1, ..., e_n]``, the number of monomials generated
  54. is:
  55. .. math::
  56. (d_1 - e_1 + 1) (d_2 - e_2 + 1) \cdots (d_n - e_n + 1)
  57. Let us generate all monomials ``monom`` in variables $x$ and $y$
  58. such that ``[1, 2][i] <= degree_list(monom)[i] <= [2, 4][i]``,
  59. ``i = 0, 1`` ::
  60. >>> from sympy import symbols
  61. >>> from sympy.polys.monomials import itermonomials
  62. >>> from sympy.polys.orderings import monomial_key
  63. >>> from sympy.abc import x, y
  64. >>> sorted(itermonomials([x, y], [2, 4], [1, 2]), reverse=True, key=monomial_key('lex', [x, y]))
  65. [x**2*y**4, x**2*y**3, x**2*y**2, x*y**4, x*y**3, x*y**2]
  66. """
  67. if is_sequence(max_degrees):
  68. n = len(variables)
  69. if len(max_degrees) != n:
  70. raise ValueError('Argument sizes do not match')
  71. if min_degrees is None:
  72. min_degrees = [0]*n
  73. elif not is_sequence(min_degrees):
  74. raise ValueError('min_degrees is not a list')
  75. else:
  76. if len(min_degrees) != n:
  77. raise ValueError('Argument sizes do not match')
  78. if any(i < 0 for i in min_degrees):
  79. raise ValueError("min_degrees cannot contain negative numbers")
  80. if any(min_degrees[i] > max_degrees[i] for i in range(n)):
  81. raise ValueError('min_degrees[i] must be <= max_degrees[i] for all i')
  82. power_lists = []
  83. for var, min_d, max_d in zip(variables, min_degrees, max_degrees):
  84. power_lists.append([var**i for i in range(min_d, max_d + 1)])
  85. for powers in product(*power_lists):
  86. yield Mul(*powers)
  87. else:
  88. max_degree = max_degrees
  89. if max_degree < 0:
  90. raise ValueError("max_degrees cannot be negative")
  91. if min_degrees is None:
  92. min_degree = 0
  93. else:
  94. if min_degrees < 0:
  95. raise ValueError("min_degrees cannot be negative")
  96. min_degree = min_degrees
  97. if min_degree > max_degree:
  98. return
  99. if not variables or max_degree == 0:
  100. yield S.One
  101. return
  102. # Force to list in case of passed tuple or other incompatible collection
  103. variables = list(variables) + [S.One]
  104. if all(variable.is_commutative for variable in variables):
  105. it = combinations_with_replacement(variables, max_degree)
  106. else:
  107. it = product(variables, repeat=max_degree)
  108. monomials_set = set()
  109. d = max_degree - min_degree
  110. for item in it:
  111. count = 0
  112. for variable in item:
  113. if variable == 1:
  114. count += 1
  115. if d < count:
  116. break
  117. else:
  118. monomials_set.add(Mul(*item))
  119. yield from monomials_set
  120. def monomial_count(V, N):
  121. r"""
  122. Computes the number of monomials.
  123. The number of monomials is given by the following formula:
  124. .. math::
  125. \frac{(\#V + N)!}{\#V! N!}
  126. where `N` is a total degree and `V` is a set of variables.
  127. Examples
  128. ========
  129. >>> from sympy.polys.monomials import itermonomials, monomial_count
  130. >>> from sympy.polys.orderings import monomial_key
  131. >>> from sympy.abc import x, y
  132. >>> monomial_count(2, 2)
  133. 6
  134. >>> M = list(itermonomials([x, y], 2))
  135. >>> sorted(M, key=monomial_key('grlex', [y, x]))
  136. [1, x, y, x**2, x*y, y**2]
  137. >>> len(M)
  138. 6
  139. """
  140. from sympy.functions.combinatorial.factorials import factorial
  141. return factorial(V + N) / factorial(V) / factorial(N)
  142. def monomial_mul(A, B):
  143. """
  144. Multiplication of tuples representing monomials.
  145. Examples
  146. ========
  147. Lets multiply `x**3*y**4*z` with `x*y**2`::
  148. >>> from sympy.polys.monomials import monomial_mul
  149. >>> monomial_mul((3, 4, 1), (1, 2, 0))
  150. (4, 6, 1)
  151. which gives `x**4*y**5*z`.
  152. """
  153. return tuple([ a + b for a, b in zip(A, B) ])
  154. def monomial_div(A, B):
  155. """
  156. Division of tuples representing monomials.
  157. Examples
  158. ========
  159. Lets divide `x**3*y**4*z` by `x*y**2`::
  160. >>> from sympy.polys.monomials import monomial_div
  161. >>> monomial_div((3, 4, 1), (1, 2, 0))
  162. (2, 2, 1)
  163. which gives `x**2*y**2*z`. However::
  164. >>> monomial_div((3, 4, 1), (1, 2, 2)) is None
  165. True
  166. `x*y**2*z**2` does not divide `x**3*y**4*z`.
  167. """
  168. C = monomial_ldiv(A, B)
  169. if all(c >= 0 for c in C):
  170. return tuple(C)
  171. else:
  172. return None
  173. def monomial_ldiv(A, B):
  174. """
  175. Division of tuples representing monomials.
  176. Examples
  177. ========
  178. Lets divide `x**3*y**4*z` by `x*y**2`::
  179. >>> from sympy.polys.monomials import monomial_ldiv
  180. >>> monomial_ldiv((3, 4, 1), (1, 2, 0))
  181. (2, 2, 1)
  182. which gives `x**2*y**2*z`.
  183. >>> monomial_ldiv((3, 4, 1), (1, 2, 2))
  184. (2, 2, -1)
  185. which gives `x**2*y**2*z**-1`.
  186. """
  187. return tuple([ a - b for a, b in zip(A, B) ])
  188. def monomial_pow(A, n):
  189. """Return the n-th pow of the monomial. """
  190. return tuple([ a*n for a in A ])
  191. def monomial_gcd(A, B):
  192. """
  193. Greatest common divisor of tuples representing monomials.
  194. Examples
  195. ========
  196. Lets compute GCD of `x*y**4*z` and `x**3*y**2`::
  197. >>> from sympy.polys.monomials import monomial_gcd
  198. >>> monomial_gcd((1, 4, 1), (3, 2, 0))
  199. (1, 2, 0)
  200. which gives `x*y**2`.
  201. """
  202. return tuple([ min(a, b) for a, b in zip(A, B) ])
  203. def monomial_lcm(A, B):
  204. """
  205. Least common multiple of tuples representing monomials.
  206. Examples
  207. ========
  208. Lets compute LCM of `x*y**4*z` and `x**3*y**2`::
  209. >>> from sympy.polys.monomials import monomial_lcm
  210. >>> monomial_lcm((1, 4, 1), (3, 2, 0))
  211. (3, 4, 1)
  212. which gives `x**3*y**4*z`.
  213. """
  214. return tuple([ max(a, b) for a, b in zip(A, B) ])
  215. def monomial_divides(A, B):
  216. """
  217. Does there exist a monomial X such that XA == B?
  218. Examples
  219. ========
  220. >>> from sympy.polys.monomials import monomial_divides
  221. >>> monomial_divides((1, 2), (3, 4))
  222. True
  223. >>> monomial_divides((1, 2), (0, 2))
  224. False
  225. """
  226. return all(a <= b for a, b in zip(A, B))
  227. def monomial_max(*monoms):
  228. """
  229. Returns maximal degree for each variable in a set of monomials.
  230. Examples
  231. ========
  232. Consider monomials `x**3*y**4*z**5`, `y**5*z` and `x**6*y**3*z**9`.
  233. We wish to find out what is the maximal degree for each of `x`, `y`
  234. and `z` variables::
  235. >>> from sympy.polys.monomials import monomial_max
  236. >>> monomial_max((3,4,5), (0,5,1), (6,3,9))
  237. (6, 5, 9)
  238. """
  239. M = list(monoms[0])
  240. for N in monoms[1:]:
  241. for i, n in enumerate(N):
  242. M[i] = max(M[i], n)
  243. return tuple(M)
  244. def monomial_min(*monoms):
  245. """
  246. Returns minimal degree for each variable in a set of monomials.
  247. Examples
  248. ========
  249. Consider monomials `x**3*y**4*z**5`, `y**5*z` and `x**6*y**3*z**9`.
  250. We wish to find out what is the minimal degree for each of `x`, `y`
  251. and `z` variables::
  252. >>> from sympy.polys.monomials import monomial_min
  253. >>> monomial_min((3,4,5), (0,5,1), (6,3,9))
  254. (0, 3, 1)
  255. """
  256. M = list(monoms[0])
  257. for N in monoms[1:]:
  258. for i, n in enumerate(N):
  259. M[i] = min(M[i], n)
  260. return tuple(M)
  261. def monomial_deg(M):
  262. """
  263. Returns the total degree of a monomial.
  264. Examples
  265. ========
  266. The total degree of `xy^2` is 3:
  267. >>> from sympy.polys.monomials import monomial_deg
  268. >>> monomial_deg((1, 2))
  269. 3
  270. """
  271. return sum(M)
  272. def term_div(a, b, domain):
  273. """Division of two terms in over a ring/field. """
  274. a_lm, a_lc = a
  275. b_lm, b_lc = b
  276. monom = monomial_div(a_lm, b_lm)
  277. if domain.is_Field:
  278. if monom is not None:
  279. return monom, domain.quo(a_lc, b_lc)
  280. else:
  281. return None
  282. else:
  283. if not (monom is None or a_lc % b_lc):
  284. return monom, domain.quo(a_lc, b_lc)
  285. else:
  286. return None
  287. class MonomialOps:
  288. """Code generator of fast monomial arithmetic functions. """
  289. @cacheit
  290. def __new__(cls, ngens):
  291. obj = super().__new__(cls)
  292. obj.ngens = ngens
  293. return obj
  294. def __getnewargs__(self):
  295. return (self.ngens,)
  296. def _build(self, code, name):
  297. ns = {}
  298. exec(code, ns)
  299. return ns[name]
  300. def _vars(self, name):
  301. return [ "%s%s" % (name, i) for i in range(self.ngens) ]
  302. @cacheit
  303. def mul(self):
  304. name = "monomial_mul"
  305. template = dedent("""\
  306. def %(name)s(A, B):
  307. (%(A)s,) = A
  308. (%(B)s,) = B
  309. return (%(AB)s,)
  310. """)
  311. A = self._vars("a")
  312. B = self._vars("b")
  313. AB = [ "%s + %s" % (a, b) for a, b in zip(A, B) ]
  314. code = template % {"name": name, "A": ", ".join(A), "B": ", ".join(B), "AB": ", ".join(AB)}
  315. return self._build(code, name)
  316. @cacheit
  317. def pow(self):
  318. name = "monomial_pow"
  319. template = dedent("""\
  320. def %(name)s(A, k):
  321. (%(A)s,) = A
  322. return (%(Ak)s,)
  323. """)
  324. A = self._vars("a")
  325. Ak = [ "%s*k" % a for a in A ]
  326. code = template % {"name": name, "A": ", ".join(A), "Ak": ", ".join(Ak)}
  327. return self._build(code, name)
  328. @cacheit
  329. def mulpow(self):
  330. name = "monomial_mulpow"
  331. template = dedent("""\
  332. def %(name)s(A, B, k):
  333. (%(A)s,) = A
  334. (%(B)s,) = B
  335. return (%(ABk)s,)
  336. """)
  337. A = self._vars("a")
  338. B = self._vars("b")
  339. ABk = [ "%s + %s*k" % (a, b) for a, b in zip(A, B) ]
  340. code = template % {"name": name, "A": ", ".join(A), "B": ", ".join(B), "ABk": ", ".join(ABk)}
  341. return self._build(code, name)
  342. @cacheit
  343. def ldiv(self):
  344. name = "monomial_ldiv"
  345. template = dedent("""\
  346. def %(name)s(A, B):
  347. (%(A)s,) = A
  348. (%(B)s,) = B
  349. return (%(AB)s,)
  350. """)
  351. A = self._vars("a")
  352. B = self._vars("b")
  353. AB = [ "%s - %s" % (a, b) for a, b in zip(A, B) ]
  354. code = template % {"name": name, "A": ", ".join(A), "B": ", ".join(B), "AB": ", ".join(AB)}
  355. return self._build(code, name)
  356. @cacheit
  357. def div(self):
  358. name = "monomial_div"
  359. template = dedent("""\
  360. def %(name)s(A, B):
  361. (%(A)s,) = A
  362. (%(B)s,) = B
  363. %(RAB)s
  364. return (%(R)s,)
  365. """)
  366. A = self._vars("a")
  367. B = self._vars("b")
  368. RAB = [ "r%(i)s = a%(i)s - b%(i)s\n if r%(i)s < 0: return None" % {"i": i} for i in range(self.ngens) ]
  369. R = self._vars("r")
  370. code = template % {"name": name, "A": ", ".join(A), "B": ", ".join(B), "RAB": "\n ".join(RAB), "R": ", ".join(R)}
  371. return self._build(code, name)
  372. @cacheit
  373. def lcm(self):
  374. name = "monomial_lcm"
  375. template = dedent("""\
  376. def %(name)s(A, B):
  377. (%(A)s,) = A
  378. (%(B)s,) = B
  379. return (%(AB)s,)
  380. """)
  381. A = self._vars("a")
  382. B = self._vars("b")
  383. AB = [ "%s if %s >= %s else %s" % (a, a, b, b) for a, b in zip(A, B) ]
  384. code = template % {"name": name, "A": ", ".join(A), "B": ", ".join(B), "AB": ", ".join(AB)}
  385. return self._build(code, name)
  386. @cacheit
  387. def gcd(self):
  388. name = "monomial_gcd"
  389. template = dedent("""\
  390. def %(name)s(A, B):
  391. (%(A)s,) = A
  392. (%(B)s,) = B
  393. return (%(AB)s,)
  394. """)
  395. A = self._vars("a")
  396. B = self._vars("b")
  397. AB = [ "%s if %s <= %s else %s" % (a, a, b, b) for a, b in zip(A, B) ]
  398. code = template % {"name": name, "A": ", ".join(A), "B": ", ".join(B), "AB": ", ".join(AB)}
  399. return self._build(code, name)
  400. @public
  401. class Monomial(PicklableWithSlots):
  402. """Class representing a monomial, i.e. a product of powers. """
  403. __slots__ = ('exponents', 'gens')
  404. def __init__(self, monom, gens=None):
  405. if not iterable(monom):
  406. rep, gens = dict_from_expr(sympify(monom), gens=gens)
  407. if len(rep) == 1 and list(rep.values())[0] == 1:
  408. monom = list(rep.keys())[0]
  409. else:
  410. raise ValueError("Expected a monomial got {}".format(monom))
  411. self.exponents = tuple(map(int, monom))
  412. self.gens = gens
  413. def rebuild(self, exponents, gens=None):
  414. return self.__class__(exponents, gens or self.gens)
  415. def __len__(self):
  416. return len(self.exponents)
  417. def __iter__(self):
  418. return iter(self.exponents)
  419. def __getitem__(self, item):
  420. return self.exponents[item]
  421. def __hash__(self):
  422. return hash((self.__class__.__name__, self.exponents, self.gens))
  423. def __str__(self):
  424. if self.gens:
  425. return "*".join([ "%s**%s" % (gen, exp) for gen, exp in zip(self.gens, self.exponents) ])
  426. else:
  427. return "%s(%s)" % (self.__class__.__name__, self.exponents)
  428. def as_expr(self, *gens):
  429. """Convert a monomial instance to a SymPy expression. """
  430. gens = gens or self.gens
  431. if not gens:
  432. raise ValueError(
  433. "Cannot convert %s to an expression without generators" % self)
  434. return Mul(*[ gen**exp for gen, exp in zip(gens, self.exponents) ])
  435. def __eq__(self, other):
  436. if isinstance(other, Monomial):
  437. exponents = other.exponents
  438. elif isinstance(other, (tuple, Tuple)):
  439. exponents = other
  440. else:
  441. return False
  442. return self.exponents == exponents
  443. def __ne__(self, other):
  444. return not self == other
  445. def __mul__(self, other):
  446. if isinstance(other, Monomial):
  447. exponents = other.exponents
  448. elif isinstance(other, (tuple, Tuple)):
  449. exponents = other
  450. else:
  451. raise NotImplementedError
  452. return self.rebuild(monomial_mul(self.exponents, exponents))
  453. def __truediv__(self, other):
  454. if isinstance(other, Monomial):
  455. exponents = other.exponents
  456. elif isinstance(other, (tuple, Tuple)):
  457. exponents = other
  458. else:
  459. raise NotImplementedError
  460. result = monomial_div(self.exponents, exponents)
  461. if result is not None:
  462. return self.rebuild(result)
  463. else:
  464. raise ExactQuotientFailed(self, Monomial(other))
  465. __floordiv__ = __truediv__
  466. def __pow__(self, other):
  467. n = int(other)
  468. if n < 0:
  469. raise ValueError("a non-negative integer expected, got %s" % other)
  470. return self.rebuild(monomial_pow(self.exponents, n))
  471. def gcd(self, other):
  472. """Greatest common divisor of monomials. """
  473. if isinstance(other, Monomial):
  474. exponents = other.exponents
  475. elif isinstance(other, (tuple, Tuple)):
  476. exponents = other
  477. else:
  478. raise TypeError(
  479. "an instance of Monomial class expected, got %s" % other)
  480. return self.rebuild(monomial_gcd(self.exponents, exponents))
  481. def lcm(self, other):
  482. """Least common multiple of monomials. """
  483. if isinstance(other, Monomial):
  484. exponents = other.exponents
  485. elif isinstance(other, (tuple, Tuple)):
  486. exponents = other
  487. else:
  488. raise TypeError(
  489. "an instance of Monomial class expected, got %s" % other)
  490. return self.rebuild(monomial_lcm(self.exponents, exponents))