""" This module contains functions for two multivariate resultants. These are: - Dixon's resultant. - Macaulay's resultant. Multivariate resultants are used to identify whether a multivariate system has common roots. That is when the resultant is equal to zero. """ from sympy.core.mul import (Mul, prod) from sympy.matrices.dense import (Matrix, diag) from sympy.polys.polytools import (Poly, degree_list, rem) from sympy.simplify.simplify import simplify from sympy.tensor.indexed import IndexedBase from sympy.polys.monomials import itermonomials, monomial_deg from sympy.polys.orderings import monomial_key from sympy.polys.polytools import poly_from_expr, total_degree from sympy.functions.combinatorial.factorials import binomial from itertools import combinations_with_replacement from sympy.utilities.exceptions import sympy_deprecation_warning class DixonResultant(): """ A class for retrieving the Dixon's resultant of a multivariate system. Examples ======== >>> from sympy import symbols >>> from sympy.polys.multivariate_resultants import DixonResultant >>> x, y = symbols('x, y') >>> p = x + y >>> q = x ** 2 + y ** 3 >>> h = x ** 2 + y >>> dixon = DixonResultant(variables=[x, y], polynomials=[p, q, h]) >>> poly = dixon.get_dixon_polynomial() >>> matrix = dixon.get_dixon_matrix(polynomial=poly) >>> matrix Matrix([ [ 0, 0, -1, 0, -1], [ 0, -1, 0, -1, 0], [-1, 0, 1, 0, 0], [ 0, -1, 0, 0, 1], [-1, 0, 0, 1, 0]]) >>> matrix.det() 0 See Also ======== Notebook in examples: sympy/example/notebooks. References ========== .. [1] [Kapur1994]_ .. [2] [Palancz08]_ """ def __init__(self, polynomials, variables): """ A class that takes two lists, a list of polynomials and list of variables. Returns the Dixon matrix of the multivariate system. Parameters ---------- polynomials : list of polynomials A list of m n-degree polynomials variables: list A list of all n variables """ self.polynomials = polynomials self.variables = variables self.n = len(self.variables) self.m = len(self.polynomials) a = IndexedBase("alpha") # A list of n alpha variables (the replacing variables) self.dummy_variables = [a[i] for i in range(self.n)] # A list of the d_max of each variable. self._max_degrees = [max(degree_list(poly)[i] for poly in self.polynomials) for i in range(self.n)] @property def max_degrees(self): sympy_deprecation_warning( """ The max_degrees property of DixonResultant is deprecated. """, deprecated_since_version="1.5", active_deprecations_target="deprecated-dixonresultant-properties", ) return self._max_degrees def get_dixon_polynomial(self): r""" Returns ======= dixon_polynomial: polynomial Dixon's polynomial is calculated as: delta = Delta(A) / ((x_1 - a_1) ... (x_n - a_n)) where, A = |p_1(x_1,... x_n), ..., p_n(x_1,... x_n)| |p_1(a_1,... x_n), ..., p_n(a_1,... x_n)| |... , ..., ...| |p_1(a_1,... a_n), ..., p_n(a_1,... a_n)| """ if self.m != (self.n + 1): raise ValueError('Method invalid for given combination.') # First row rows = [self.polynomials] temp = list(self.variables) for idx in range(self.n): temp[idx] = self.dummy_variables[idx] substitution = {var: t for var, t in zip(self.variables, temp)} rows.append([f.subs(substitution) for f in self.polynomials]) A = Matrix(rows) terms = zip(self.variables, self.dummy_variables) product_of_differences = Mul(*[a - b for a, b in terms]) dixon_polynomial = (A.det() / product_of_differences).factor() return poly_from_expr(dixon_polynomial, self.dummy_variables)[0] def get_upper_degree(self): sympy_deprecation_warning( """ The get_upper_degree() method of DixonResultant is deprecated. Use get_max_degrees() instead. """, deprecated_since_version="1.5", active_deprecations_target="deprecated-dixonresultant-properties" ) list_of_products = [self.variables[i] ** self._max_degrees[i] for i in range(self.n)] product = prod(list_of_products) product = Poly(product).monoms() return monomial_deg(*product) def get_max_degrees(self, polynomial): r""" Returns a list of the maximum degree of each variable appearing in the coefficients of the Dixon polynomial. The coefficients are viewed as polys in $x_1, x_2, \dots, x_n$. """ deg_lists = [degree_list(Poly(poly, self.variables)) for poly in polynomial.coeffs()] max_degrees = [max(degs) for degs in zip(*deg_lists)] return max_degrees def get_dixon_matrix(self, polynomial): r""" Construct the Dixon matrix from the coefficients of polynomial \alpha. Each coefficient is viewed as a polynomial of x_1, ..., x_n. """ max_degrees = self.get_max_degrees(polynomial) # list of column headers of the Dixon matrix. monomials = itermonomials(self.variables, max_degrees) monomials = sorted(monomials, reverse=True, key=monomial_key('lex', self.variables)) dixon_matrix = Matrix([[Poly(c, *self.variables).coeff_monomial(m) for m in monomials] for c in polynomial.coeffs()]) # remove columns if needed if dixon_matrix.shape[0] != dixon_matrix.shape[1]: keep = [column for column in range(dixon_matrix.shape[-1]) if any(element != 0 for element in dixon_matrix[:, column])] dixon_matrix = dixon_matrix[:, keep] return dixon_matrix def KSY_precondition(self, matrix): """ Test for the validity of the Kapur-Saxena-Yang precondition. The precondition requires that the column corresponding to the monomial 1 = x_1 ^ 0 * x_2 ^ 0 * ... * x_n ^ 0 is not a linear combination of the remaining ones. In SymPy notation this is the last column. For the precondition to hold the last non-zero row of the rref matrix should be of the form [0, 0, ..., 1]. """ if matrix.is_zero_matrix: return False m, n = matrix.shape # simplify the matrix and keep only its non-zero rows matrix = simplify(matrix.rref()[0]) rows = [i for i in range(m) if any(matrix[i, j] != 0 for j in range(n))] matrix = matrix[rows,:] condition = Matrix([[0]*(n-1) + [1]]) if matrix[-1,:] == condition: return True else: return False def delete_zero_rows_and_columns(self, matrix): """Remove the zero rows and columns of the matrix.""" rows = [ i for i in range(matrix.rows) if not matrix.row(i).is_zero_matrix] cols = [ j for j in range(matrix.cols) if not matrix.col(j).is_zero_matrix] return matrix[rows, cols] def product_leading_entries(self, matrix): """Calculate the product of the leading entries of the matrix.""" res = 1 for row in range(matrix.rows): for el in matrix.row(row): if el != 0: res = res * el break return res def get_KSY_Dixon_resultant(self, matrix): """Calculate the Kapur-Saxena-Yang approach to the Dixon Resultant.""" matrix = self.delete_zero_rows_and_columns(matrix) _, U, _ = matrix.LUdecomposition() matrix = self.delete_zero_rows_and_columns(simplify(U)) return self.product_leading_entries(matrix) class MacaulayResultant(): """ A class for calculating the Macaulay resultant. Note that the polynomials must be homogenized and their coefficients must be given as symbols. Examples ======== >>> from sympy import symbols >>> from sympy.polys.multivariate_resultants import MacaulayResultant >>> x, y, z = symbols('x, y, z') >>> a_0, a_1, a_2 = symbols('a_0, a_1, a_2') >>> b_0, b_1, b_2 = symbols('b_0, b_1, b_2') >>> c_0, c_1, c_2,c_3, c_4 = symbols('c_0, c_1, c_2, c_3, c_4') >>> f = a_0 * y - a_1 * x + a_2 * z >>> g = b_1 * x ** 2 + b_0 * y ** 2 - b_2 * z ** 2 >>> h = c_0 * y * z ** 2 - c_1 * x ** 3 + c_2 * x ** 2 * z - c_3 * x * z ** 2 + c_4 * z ** 3 >>> mac = MacaulayResultant(polynomials=[f, g, h], variables=[x, y, z]) >>> mac.monomial_set [x**4, x**3*y, x**3*z, x**2*y**2, x**2*y*z, x**2*z**2, x*y**3, x*y**2*z, x*y*z**2, x*z**3, y**4, y**3*z, y**2*z**2, y*z**3, z**4] >>> matrix = mac.get_matrix() >>> submatrix = mac.get_submatrix(matrix) >>> submatrix Matrix([ [-a_1, a_0, a_2, 0], [ 0, -a_1, 0, 0], [ 0, 0, -a_1, 0], [ 0, 0, 0, -a_1]]) See Also ======== Notebook in examples: sympy/example/notebooks. References ========== .. [1] [Bruce97]_ .. [2] [Stiller96]_ """ def __init__(self, polynomials, variables): """ Parameters ========== variables: list A list of all n variables polynomials : list of SymPy polynomials A list of m n-degree polynomials """ self.polynomials = polynomials self.variables = variables self.n = len(variables) # A list of the d_max of each polynomial. self.degrees = [total_degree(poly, *self.variables) for poly in self.polynomials] self.degree_m = self._get_degree_m() self.monomials_size = self.get_size() # The set T of all possible monomials of degree degree_m self.monomial_set = self.get_monomials_of_certain_degree(self.degree_m) def _get_degree_m(self): r""" Returns ======= degree_m: int The degree_m is calculated as 1 + \sum_1 ^ n (d_i - 1), where d_i is the degree of the i polynomial """ return 1 + sum(d - 1 for d in self.degrees) def get_size(self): r""" Returns ======= size: int The size of set T. Set T is the set of all possible monomials of the n variables for degree equal to the degree_m """ return binomial(self.degree_m + self.n - 1, self.n - 1) def get_monomials_of_certain_degree(self, degree): """ Returns ======= monomials: list A list of monomials of a certain degree. """ monomials = [Mul(*monomial) for monomial in combinations_with_replacement(self.variables, degree)] return sorted(monomials, reverse=True, key=monomial_key('lex', self.variables)) def get_row_coefficients(self): """ Returns ======= row_coefficients: list The row coefficients of Macaulay's matrix """ row_coefficients = [] divisible = [] for i in range(self.n): if i == 0: degree = self.degree_m - self.degrees[i] monomial = self.get_monomials_of_certain_degree(degree) row_coefficients.append(monomial) else: divisible.append(self.variables[i - 1] ** self.degrees[i - 1]) degree = self.degree_m - self.degrees[i] poss_rows = self.get_monomials_of_certain_degree(degree) for div in divisible: for p in poss_rows: if rem(p, div) == 0: poss_rows = [item for item in poss_rows if item != p] row_coefficients.append(poss_rows) return row_coefficients def get_matrix(self): """ Returns ======= macaulay_matrix: Matrix The Macaulay numerator matrix """ rows = [] row_coefficients = self.get_row_coefficients() for i in range(self.n): for multiplier in row_coefficients[i]: coefficients = [] poly = Poly(self.polynomials[i] * multiplier, *self.variables) for mono in self.monomial_set: coefficients.append(poly.coeff_monomial(mono)) rows.append(coefficients) macaulay_matrix = Matrix(rows) return macaulay_matrix def get_reduced_nonreduced(self): r""" Returns ======= reduced: list A list of the reduced monomials non_reduced: list A list of the monomials that are not reduced Definition ========== A polynomial is said to be reduced in x_i, if its degree (the maximum degree of its monomials) in x_i is less than d_i. A polynomial that is reduced in all variables but one is said simply to be reduced. """ divisible = [] for m in self.monomial_set: temp = [] for i, v in enumerate(self.variables): temp.append(bool(total_degree(m, v) >= self.degrees[i])) divisible.append(temp) reduced = [i for i, r in enumerate(divisible) if sum(r) < self.n - 1] non_reduced = [i for i, r in enumerate(divisible) if sum(r) >= self.n -1] return reduced, non_reduced def get_submatrix(self, matrix): r""" Returns ======= macaulay_submatrix: Matrix The Macaulay denominator matrix. Columns that are non reduced are kept. The row which contains one of the a_{i}s is dropped. a_{i}s are the coefficients of x_i ^ {d_i}. """ reduced, non_reduced = self.get_reduced_nonreduced() # if reduced == [], then det(matrix) should be 1 if reduced == []: return diag([1]) # reduced != [] reduction_set = [v ** self.degrees[i] for i, v in enumerate(self.variables)] ais = list([self.polynomials[i].coeff(reduction_set[i]) for i in range(self.n)]) reduced_matrix = matrix[:, reduced] keep = [] for row in range(reduced_matrix.rows): check = [ai in reduced_matrix[row, :] for ai in ais] if True not in check: keep.append(row) return matrix[keep, non_reduced]