474 lines
15 KiB
Python
474 lines
15 KiB
Python
|
"""
|
||
|
|
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 math import prod
|
||
|
|
|
||
|
|
from sympy.core.mul import Mul
|
||
|
|
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 = [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]
|