Source code for diofant.polys.monomials

"""Tools and arithmetics for monomials of distributed polynomials. """

from ..core import Mul, S, Tuple, sympify
from ..core.compatibility import iterable
from .polyerrors import ExactQuotientFailed
from .polyutils import dict_from_expr


__all__ = 'Monomial', 'itermonomials'


[docs]def itermonomials(variables, degree): r""" Generate a set of monomials of the given total degree or less. Given a set of variables `V` and a total degree `N` generate a set of monomials of degree at most `N`. The total number of monomials is huge and is given by the following formula: .. math:: \frac{(\#V + N)!}{\#V! N!} For example if we would like to generate a dense polynomial of a total degree `N = 50` 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 variables `x` and `y`:: >>> from diofant.polys.orderings import monomial_key >>> 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] """ if not variables: return {S.One} else: x, tail = variables[0], variables[1:] monoms = itermonomials(tail, degree) for i in range(1, degree + 1): monoms |= { x**i * m for m in itermonomials(tail, degree - i) } return monoms
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 diofant.polys.orderings import monomial_key >>> 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 ..functions import factorial return factorial(V + N) / factorial(V) / factorial(N) def monomial_mul(A, B): """ Multiplication of tuples representing monomials. Lets multiply `x**3*y**4*z` with `x*y**2`:: >>> 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. Lets divide `x**3*y**4*z` by `x*y**2`:: >>> 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 def monomial_ldiv(A, B): """ Division of tuples representing monomials. Lets divide `x**3*y**4*z` by `x*y**2`:: >>> 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. Lets compute GCD of `x*y**4*z` and `x**3*y**2`:: >>> 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. Lets compute LCM of `x*y**4*z` and `x**3*y**2`:: >>> 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? >>> 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. 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:: >>> 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. 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:: >>> 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. For example, the total degree of `xy^2` is 3: >>> 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: if not (monom is None or a_lc % b_lc): return monom, domain.quo(a_lc, b_lc)
[docs]class Monomial: """Class representing a monomial, i.e. a product of powers. """ 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)[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)
[docs] def as_expr(self, *gens): """Convert a monomial instance to a Diofant 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 __mul__(self, other): if isinstance(other, Monomial): exponents = other.exponents elif isinstance(other, (tuple, Tuple)): exponents = other else: return NotImplemented return self.rebuild(monomial_mul(self.exponents, exponents)) def __truediv__(self, other): if isinstance(other, Monomial): exponents = other.exponents elif isinstance(other, (tuple, Tuple)): exponents = other else: return NotImplemented result = monomial_div(self.exponents, exponents) if result is not None: return self.rebuild(result) else: raise ExactQuotientFailed(self, Monomial(other)) __floordiv__ = __truediv__ 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)
[docs] 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))
[docs] 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))