"""Tools and arithmetics for monomials of distributed polynomials. """
from __future__ import print_function, division
from itertools import combinations_with_replacement, product
from textwrap import dedent
from sympy.core import Mul, S, Tuple, sympify
from sympy.core.compatibility import exec_, iterable, range
from sympy.polys.polyerrors import ExactQuotientFailed
from sympy.polys.polyutils import PicklableWithSlots, dict_from_expr
from sympy.utilities import public
[docs]@public
def itermonomials(variables, max_degree, min_degree = 0):
r"""
Generate a set of monomials of the degree greater than or equal
to `min_degree` and less than or equal to `max_degree`.
Given a set of variables `V` and a min_degree `N` and a max_degree `M`
generate a set of monomials of degree less than or equal to `N` and greater
than or equal to `M`. The total number of monomials in commutative
variables is huge and is given by the following formula if `M = 0`:
.. math::
\frac{(\#V + N)!}{\#V! N!}
For example if we would like to generate a dense polynomial of
a total degree `N = 50` and `M = 0`, which is the worst case, in 5
variables, assuming that exponents and all of coefficients are 32-bit long
and stored in an array we would need almost 80 GiB of memory! Fortunately
most polynomials, that we will encounter, are sparse.
Examples
========
Consider monomials in commutative variables `x` and `y`
and non-commutative variables `a` and `b`::
>>> from sympy import symbols
>>> from sympy.polys.monomials import itermonomials
>>> from sympy.polys.orderings import monomial_key
>>> from sympy.abc import x, y
>>> sorted(itermonomials([x, y], 2), key=monomial_key('grlex', [y, x]))
[1, x, y, x**2, x*y, y**2]
>>> sorted(itermonomials([x, y], 3), key=monomial_key('grlex', [y, x]))
[1, x, y, x**2, x*y, y**2, x**3, x**2*y, x*y**2, y**3]
>>> a, b = symbols('a, b', commutative=False)
>>> itermonomials([a, b, x], 2)
{1, a, a**2, b, b**2, x, x**2, a*b, b*a, x*a, x*b}
>>> sorted(itermonomials([x, y], 2, 1), key=monomial_key('grlex', [y, x]))
[x, y, x**2, x*y, y**2]
"""
if max_degree < 0 or min_degree > max_degree:
return set()
if not variables or max_degree == 0:
return {S(1)}
# Force to list in case of passed tuple or other incompatible collection
variables = list(variables) + [S(1)]
if all(variable.is_commutative for variable in variables):
monomials_list_comm = []
for item in combinations_with_replacement(variables, max_degree):
powers = dict()
for variable in variables:
powers[variable] = 0
for variable in item:
if variable != 1:
powers[variable] += 1
if max(powers.values()) >= min_degree:
monomials_list_comm.append(Mul(*item))
return set(monomials_list_comm)
else:
monomials_list_non_comm = []
for item in product(variables, repeat=max_degree):
powers = dict()
for variable in variables:
powers[variable] = 0
for variable in item:
if variable != 1:
powers[variable] += 1
if max(powers.values()) >= min_degree:
monomials_list_non_comm.append(Mul(*item))
return set(monomials_list_non_comm)
[docs]def monomial_count(V, N):
r"""
Computes the number of monomials.
The number of monomials is given by the following formula:
.. math::
\frac{(\#V + N)!}{\#V! N!}
where `N` is a total degree and `V` is a set of variables.
Examples
========
>>> from sympy.polys.monomials import itermonomials, monomial_count
>>> from sympy.polys.orderings import monomial_key
>>> from sympy.abc import x, y
>>> monomial_count(2, 2)
6
>>> M = itermonomials([x, y], 2)
>>> sorted(M, key=monomial_key('grlex', [y, x]))
[1, x, y, x**2, x*y, y**2]
>>> len(M)
6
"""
from sympy import factorial
return factorial(V + N) / factorial(V) / factorial(N)
def monomial_mul(A, B):
"""
Multiplication of tuples representing monomials.
Examples
========
Lets multiply `x**3*y**4*z` with `x*y**2`::
>>> from sympy.polys.monomials import monomial_mul
>>> monomial_mul((3, 4, 1), (1, 2, 0))
(4, 6, 1)
which gives `x**4*y**5*z`.
"""
return tuple([ a + b for a, b in zip(A, B) ])
def monomial_div(A, B):
"""
Division of tuples representing monomials.
Examples
========
Lets divide `x**3*y**4*z` by `x*y**2`::
>>> from sympy.polys.monomials import monomial_div
>>> monomial_div((3, 4, 1), (1, 2, 0))
(2, 2, 1)
which gives `x**2*y**2*z`. However::
>>> monomial_div((3, 4, 1), (1, 2, 2)) is None
True
`x*y**2*z**2` does not divide `x**3*y**4*z`.
"""
C = monomial_ldiv(A, B)
if all(c >= 0 for c in C):
return tuple(C)
else:
return None
def monomial_ldiv(A, B):
"""
Division of tuples representing monomials.
Examples
========
Lets divide `x**3*y**4*z` by `x*y**2`::
>>> from sympy.polys.monomials import monomial_ldiv
>>> monomial_ldiv((3, 4, 1), (1, 2, 0))
(2, 2, 1)
which gives `x**2*y**2*z`.
>>> monomial_ldiv((3, 4, 1), (1, 2, 2))
(2, 2, -1)
which gives `x**2*y**2*z**-1`.
"""
return tuple([ a - b for a, b in zip(A, B) ])
def monomial_pow(A, n):
"""Return the n-th pow of the monomial. """
return tuple([ a*n for a in A ])
def monomial_gcd(A, B):
"""
Greatest common divisor of tuples representing monomials.
Examples
========
Lets compute GCD of `x*y**4*z` and `x**3*y**2`::
>>> from sympy.polys.monomials import monomial_gcd
>>> monomial_gcd((1, 4, 1), (3, 2, 0))
(1, 2, 0)
which gives `x*y**2`.
"""
return tuple([ min(a, b) for a, b in zip(A, B) ])
def monomial_lcm(A, B):
"""
Least common multiple of tuples representing monomials.
Examples
========
Lets compute LCM of `x*y**4*z` and `x**3*y**2`::
>>> from sympy.polys.monomials import monomial_lcm
>>> monomial_lcm((1, 4, 1), (3, 2, 0))
(3, 4, 1)
which gives `x**3*y**4*z`.
"""
return tuple([ max(a, b) for a, b in zip(A, B) ])
def monomial_divides(A, B):
"""
Does there exist a monomial X such that XA == B?
Examples
========
>>> from sympy.polys.monomials import monomial_divides
>>> monomial_divides((1, 2), (3, 4))
True
>>> monomial_divides((1, 2), (0, 2))
False
"""
return all(a <= b for a, b in zip(A, B))
def monomial_max(*monoms):
"""
Returns maximal degree for each variable in a set of monomials.
Examples
========
Consider monomials `x**3*y**4*z**5`, `y**5*z` and `x**6*y**3*z**9`.
We wish to find out what is the maximal degree for each of `x`, `y`
and `z` variables::
>>> from sympy.polys.monomials import monomial_max
>>> monomial_max((3,4,5), (0,5,1), (6,3,9))
(6, 5, 9)
"""
M = list(monoms[0])
for N in monoms[1:]:
for i, n in enumerate(N):
M[i] = max(M[i], n)
return tuple(M)
def monomial_min(*monoms):
"""
Returns minimal degree for each variable in a set of monomials.
Examples
========
Consider monomials `x**3*y**4*z**5`, `y**5*z` and `x**6*y**3*z**9`.
We wish to find out what is the minimal degree for each of `x`, `y`
and `z` variables::
>>> from sympy.polys.monomials import monomial_min
>>> monomial_min((3,4,5), (0,5,1), (6,3,9))
(0, 3, 1)
"""
M = list(monoms[0])
for N in monoms[1:]:
for i, n in enumerate(N):
M[i] = min(M[i], n)
return tuple(M)
def monomial_deg(M):
"""
Returns the total degree of a monomial.
Examples
========
The total degree of `xy^2` is 3:
>>> from sympy.polys.monomials import monomial_deg
>>> monomial_deg((1, 2))
3
"""
return sum(M)
def term_div(a, b, domain):
"""Division of two terms in over a ring/field. """
a_lm, a_lc = a
b_lm, b_lc = b
monom = monomial_div(a_lm, b_lm)
if domain.is_Field:
if monom is not None:
return monom, domain.quo(a_lc, b_lc)
else:
return None
else:
if not (monom is None or a_lc % b_lc):
return monom, domain.quo(a_lc, b_lc)
else:
return None
class MonomialOps(object):
"""Code generator of fast monomial arithmetic functions. """
def __init__(self, ngens):
self.ngens = ngens
def _build(self, code, name):
ns = {}
exec_(code, ns)
return ns[name]
def _vars(self, name):
return [ "%s%s" % (name, i) for i in range(self.ngens) ]
def mul(self):
name = "monomial_mul"
template = dedent("""\
def %(name)s(A, B):
(%(A)s,) = A
(%(B)s,) = B
return (%(AB)s,)
""")
A = self._vars("a")
B = self._vars("b")
AB = [ "%s + %s" % (a, b) for a, b in zip(A, B) ]
code = template % dict(name=name, A=", ".join(A), B=", ".join(B), AB=", ".join(AB))
return self._build(code, name)
def pow(self):
name = "monomial_pow"
template = dedent("""\
def %(name)s(A, k):
(%(A)s,) = A
return (%(Ak)s,)
""")
A = self._vars("a")
Ak = [ "%s*k" % a for a in A ]
code = template % dict(name=name, A=", ".join(A), Ak=", ".join(Ak))
return self._build(code, name)
def mulpow(self):
name = "monomial_mulpow"
template = dedent("""\
def %(name)s(A, B, k):
(%(A)s,) = A
(%(B)s,) = B
return (%(ABk)s,)
""")
A = self._vars("a")
B = self._vars("b")
ABk = [ "%s + %s*k" % (a, b) for a, b in zip(A, B) ]
code = template % dict(name=name, A=", ".join(A), B=", ".join(B), ABk=", ".join(ABk))
return self._build(code, name)
def ldiv(self):
name = "monomial_ldiv"
template = dedent("""\
def %(name)s(A, B):
(%(A)s,) = A
(%(B)s,) = B
return (%(AB)s,)
""")
A = self._vars("a")
B = self._vars("b")
AB = [ "%s - %s" % (a, b) for a, b in zip(A, B) ]
code = template % dict(name=name, A=", ".join(A), B=", ".join(B), AB=", ".join(AB))
return self._build(code, name)
def div(self):
name = "monomial_div"
template = dedent("""\
def %(name)s(A, B):
(%(A)s,) = A
(%(B)s,) = B
%(RAB)s
return (%(R)s,)
""")
A = self._vars("a")
B = self._vars("b")
RAB = [ "r%(i)s = a%(i)s - b%(i)s\n if r%(i)s < 0: return None" % dict(i=i) for i in range(self.ngens) ]
R = self._vars("r")
code = template % dict(name=name, A=", ".join(A), B=", ".join(B), RAB="\n ".join(RAB), R=", ".join(R))
return self._build(code, name)
def lcm(self):
name = "monomial_lcm"
template = dedent("""\
def %(name)s(A, B):
(%(A)s,) = A
(%(B)s,) = B
return (%(AB)s,)
""")
A = self._vars("a")
B = self._vars("b")
AB = [ "%s if %s >= %s else %s" % (a, a, b, b) for a, b in zip(A, B) ]
code = template % dict(name=name, A=", ".join(A), B=", ".join(B), AB=", ".join(AB))
return self._build(code, name)
def gcd(self):
name = "monomial_gcd"
template = dedent("""\
def %(name)s(A, B):
(%(A)s,) = A
(%(B)s,) = B
return (%(AB)s,)
""")
A = self._vars("a")
B = self._vars("b")
AB = [ "%s if %s <= %s else %s" % (a, a, b, b) for a, b in zip(A, B) ]
code = template % dict(name=name, A=", ".join(A), B=", ".join(B), AB=", ".join(AB))
return self._build(code, name)
[docs]@public
class Monomial(PicklableWithSlots):
"""Class representing a monomial, i.e. a product of powers. """
__slots__ = ['exponents', 'gens']
def __init__(self, monom, gens=None):
if not iterable(monom):
rep, gens = dict_from_expr(sympify(monom), gens=gens)
if len(rep) == 1 and list(rep.values())[0] == 1:
monom = list(rep.keys())[0]
else:
raise ValueError("Expected a monomial got %s" % monom)
self.exponents = tuple(map(int, monom))
self.gens = gens
def rebuild(self, exponents, gens=None):
return self.__class__(exponents, gens or self.gens)
def __len__(self):
return len(self.exponents)
def __iter__(self):
return iter(self.exponents)
def __getitem__(self, item):
return self.exponents[item]
def __hash__(self):
return hash((self.__class__.__name__, self.exponents, self.gens))
def __str__(self):
if self.gens:
return "*".join([ "%s**%s" % (gen, exp) for gen, exp in zip(self.gens, self.exponents) ])
else:
return "%s(%s)" % (self.__class__.__name__, self.exponents)
def as_expr(self, *gens):
"""Convert a monomial instance to a SymPy expression. """
gens = gens or self.gens
if not gens:
raise ValueError(
"can't convert %s to an expression without generators" % self)
return Mul(*[ gen**exp for gen, exp in zip(gens, self.exponents) ])
def __eq__(self, other):
if isinstance(other, Monomial):
exponents = other.exponents
elif isinstance(other, (tuple, Tuple)):
exponents = other
else:
return False
return self.exponents == exponents
def __ne__(self, other):
return not self == other
def __mul__(self, other):
if isinstance(other, Monomial):
exponents = other.exponents
elif isinstance(other, (tuple, Tuple)):
exponents = other
else:
return NotImplementedError
return self.rebuild(monomial_mul(self.exponents, exponents))
def __div__(self, other):
if isinstance(other, Monomial):
exponents = other.exponents
elif isinstance(other, (tuple, Tuple)):
exponents = other
else:
return NotImplementedError
result = monomial_div(self.exponents, exponents)
if result is not None:
return self.rebuild(result)
else:
raise ExactQuotientFailed(self, Monomial(other))
__floordiv__ = __truediv__ = __div__
def __pow__(self, other):
n = int(other)
if not n:
return self.rebuild([0]*len(self))
elif n > 0:
exponents = self.exponents
for i in range(1, n):
exponents = monomial_mul(exponents, self.exponents)
return self.rebuild(exponents)
else:
raise ValueError("a non-negative integer expected, got %s" % other)
def gcd(self, other):
"""Greatest common divisor of monomials. """
if isinstance(other, Monomial):
exponents = other.exponents
elif isinstance(other, (tuple, Tuple)):
exponents = other
else:
raise TypeError(
"an instance of Monomial class expected, got %s" % other)
return self.rebuild(monomial_gcd(self.exponents, exponents))
def lcm(self, other):
"""Least common multiple of monomials. """
if isinstance(other, Monomial):
exponents = other.exponents
elif isinstance(other, (tuple, Tuple)):
exponents = other
else:
raise TypeError(
"an instance of Monomial class expected, got %s" % other)
return self.rebuild(monomial_lcm(self.exponents, exponents))
