prjycq8up
/
image-handler
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
Branches Tags
${ item.name }
Create tag ${ searchTerm }
Create branch ${ searchTerm }
from 'master'
${ noResults }
image-handler/image-handler10.0/venv/lib/site-packages/sympy/polys/subresultants_qq_zz.py

2559 lines
86 KiB
Raw Permalink Normal View History Unescape

first commit
5 months ago
"""
This module contains functions for the computation
of Euclidean, (generalized) Sturmian, (modified) subresultant
polynomial remainder sequences (prs's) of two polynomials;
included are also three functions for the computation of the
resultant of two polynomials.
Except for the function res_z(), which computes the resultant
of two polynomials, the pseudo-remainder function prem()
of sympy is _not_ used by any of the functions in the module.
Instead of prem() we use the function
rem_z().
Included is also the function quo_z().
An explanation of why we avoid prem() can be found in the
references stated in the docstring of rem_z().
1. Theoretical background:
==========================
Consider the polynomials f, g in Z[x] of degrees deg(f) = n and
deg(g) = m with n >= m.
Definition 1:
=============
The sign sequence of a polynomial remainder sequence (prs) is the
sequence of signs of the leading coefficients of its polynomials.
Sign sequences can be computed with the function:
sign_seq(poly_seq, x)
Definition 2:
=============
A polynomial remainder sequence (prs) is called complete if the
degree difference between any two consecutive polynomials is 1;
otherwise, it called incomplete.
It is understood that f, g belong to the sequences mentioned in
the two definitions above.
1A. Euclidean and subresultant prs's:
=====================================
The subresultant prs of f, g is a sequence of polynomials in Z[x]
analogous to the Euclidean prs, the sequence obtained by applying
on f, g Euclid's algorithm for polynomial greatest common divisors
(gcd) in Q[x].
The subresultant prs differs from the Euclidean prs in that the
coefficients of each polynomial in the former sequence are determinants
--- also referred to as subresultants --- of appropriately selected
sub-matrices of sylvester1(f, g, x), Sylvester's matrix of 1840 of
dimensions (n + m) * (n + m).
Recall that the determinant of sylvester1(f, g, x) itself is
called the resultant of f, g and serves as a criterion of whether
the two polynomials have common roots or not.
In SymPy the resultant is computed with the function
resultant(f, g, x). This function does _not_ evaluate the
determinant of sylvester(f, g, x, 1); instead, it returns
the last member of the subresultant prs of f, g, multiplied
(if needed) by an appropriate power of -1; see the caveat below.
In this module we use three functions to compute the
resultant of f, g:
a) res(f, g, x) computes the resultant by evaluating
the determinant of sylvester(f, g, x, 1);
b) res_q(f, g, x) computes the resultant recursively, by
performing polynomial divisions in Q[x] with the function rem();
c) res_z(f, g, x) computes the resultant recursively, by
performing polynomial divisions in Z[x] with the function prem().
Caveat: If Df = degree(f, x) and Dg = degree(g, x), then:
resultant(f, g, x) = (-1)**(Df*Dg) * resultant(g, f, x).
For complete prs's the sign sequence of the Euclidean prs of f, g
is identical to the sign sequence of the subresultant prs of f, g
and the coefficients of one sequence are easily computed from the
coefficients of the other.
For incomplete prs's the polynomials in the subresultant prs, generally
differ in sign from those of the Euclidean prs, and --- unlike the
case of complete prs's --- it is not at all obvious how to compute
the coefficients of one sequence from the coefficients of the other.
1B. Sturmian and modified subresultant prs's:
=============================================
For the same polynomials f, g in Z[x] mentioned above, their ``modified''
subresultant prs is a sequence of polynomials similar to the Sturmian
prs, the sequence obtained by applying in Q[x] Sturm's algorithm on f, g.
The two sequences differ in that the coefficients of each polynomial
in the modified subresultant prs are the determinants --- also referred
to as modified subresultants --- of appropriately selected sub-matrices
of sylvester2(f, g, x), Sylvester's matrix of 1853 of dimensions 2n x 2n.
The determinant of sylvester2 itself is called the modified resultant
of f, g and it also can serve as a criterion of whether the two
polynomials have common roots or not.
For complete prs's the sign sequence of the Sturmian prs of f, g is
identical to the sign sequence of the modified subresultant prs of
f, g and the coefficients of one sequence are easily computed from
the coefficients of the other.
For incomplete prs's the polynomials in the modified subresultant prs,
generally differ in sign from those of the Sturmian prs, and --- unlike
the case of complete prs's --- it is not at all obvious how to compute
the coefficients of one sequence from the coefficients of the other.
As Sylvester pointed out, the coefficients of the polynomial remainders
obtained as (modified) subresultants are the smallest possible without
introducing rationals and without computing (integer) greatest common
divisors.
1C. On terminology:
===================
Whence the terminology? Well generalized Sturmian prs's are
``modifications'' of Euclidean prs's; the hint came from the title
of the Pell-Gordon paper of 1917.
In the literature one also encounters the name ``non signed'' and
``signed'' prs for Euclidean and Sturmian prs respectively.
Likewise ``non signed'' and ``signed'' subresultant prs for
subresultant and modified subresultant prs respectively.
2. Functions in the module:
===========================
No function utilizes SymPy's function prem().
2A. Matrices:
=============
The functions sylvester(f, g, x, method=1) and
sylvester(f, g, x, method=2) compute either Sylvester matrix.
They can be used to compute (modified) subresultant prs's by
direct determinant evaluation.
The function bezout(f, g, x, method='prs') provides a matrix of
smaller dimensions than either Sylvester matrix. It is the function
of choice for computing (modified) subresultant prs's by direct
determinant evaluation.
sylvester(f, g, x, method=1)
sylvester(f, g, x, method=2)
bezout(f, g, x, method='prs')
The following identity holds:
bezout(f, g, x, method='prs') =
backward_eye(deg(f))*bezout(f, g, x, method='bz')*backward_eye(deg(f))
2B. Subresultant and modified subresultant prs's by
===================================================
determinant evaluations:
=======================
We use the Sylvester matrices of 1840 and 1853 to
compute, respectively, subresultant and modified
subresultant polynomial remainder sequences. However,
for large matrices this approach takes a lot of time.
Instead of utilizing the Sylvester matrices, we can
employ the Bezout matrix which is of smaller dimensions.
subresultants_sylv(f, g, x)
modified_subresultants_sylv(f, g, x)
subresultants_bezout(f, g, x)
modified_subresultants_bezout(f, g, x)
2C. Subresultant prs's by ONE determinant evaluation:
=====================================================
All three functions in this section evaluate one determinant
per remainder polynomial; this is the determinant of an
appropriately selected sub-matrix of sylvester1(f, g, x),
Sylvester's matrix of 1840.
To compute the remainder polynomials the function
subresultants_rem(f, g, x) employs rem(f, g, x).
By contrast, the other two functions implement Van Vleck's ideas
of 1900 and compute the remainder polynomials by trinagularizing
sylvester2(f, g, x), Sylvester's matrix of 1853.
subresultants_rem(f, g, x)
subresultants_vv(f, g, x)
subresultants_vv_2(f, g, x).
2E. Euclidean, Sturmian prs's in Q[x]:
======================================
euclid_q(f, g, x)
sturm_q(f, g, x)
2F. Euclidean, Sturmian and (modified) subresultant prs's P-G:
==============================================================
All functions in this section are based on the Pell-Gordon (P-G)
theorem of 1917.
Computations are done in Q[x], employing the function rem(f, g, x)
for the computation of the remainder polynomials.
euclid_pg(f, g, x)
sturm pg(f, g, x)
subresultants_pg(f, g, x)
modified_subresultants_pg(f, g, x)
2G. Euclidean, Sturmian and (modified) subresultant prs's A-M-V:
================================================================
All functions in this section are based on the Akritas-Malaschonok-
Vigklas (A-M-V) theorem of 2015.
Computations are done in Z[x], employing the function rem_z(f, g, x)
for the computation of the remainder polynomials.
euclid_amv(f, g, x)
sturm_amv(f, g, x)
subresultants_amv(f, g, x)
modified_subresultants_amv(f, g, x)
2Ga. Exception:
===============
subresultants_amv_q(f, g, x)
This function employs rem(f, g, x) for the computation of
the remainder polynomials, despite the fact that it implements
the A-M-V Theorem.
It is included in our module in order to show that theorems P-G
and A-M-V can be implemented utilizing either the function
rem(f, g, x) or the function rem_z(f, g, x).
For clearly historical reasons --- since the Collins-Brown-Traub
coefficients-reduction factor beta_i was not available in 1917 ---
we have implemented the Pell-Gordon theorem with the function
rem(f, g, x) and the A-M-V Theorem with the function rem_z(f, g, x).
2H. Resultants:
===============
res(f, g, x)
res_q(f, g, x)
res_z(f, g, x)
"""
from sympy.concrete.summations import summation
from sympy.core.function import expand
from sympy.core.numbers import nan
from sympy.core.singleton import S
from sympy.core.symbol import Dummy as var
from sympy.functions.elementary.complexes import Abs, sign
from sympy.functions.elementary.integers import floor
from sympy.matrices.dense import eye, Matrix, zeros
from sympy.printing.pretty.pretty import pretty_print as pprint
from sympy.simplify.simplify import simplify
from sympy.polys.domains import QQ
from sympy.polys.polytools import degree, LC, Poly, pquo, quo, prem, rem
from sympy.polys.polyerrors import PolynomialError
def sylvester(f, g, x, method = 1):
'''
The input polynomials f, g are in Z[x] or in Q[x]. Let m = degree(f, x),
n = degree(g, x) and mx = max(m, n).
a. If method = 1 (default), computes sylvester1, Sylvester's matrix of 1840
of dimension (m + n) x (m + n). The determinants of properly chosen
submatrices of this matrix (a.k.a. subresultants) can be
used to compute the coefficients of the Euclidean PRS of f, g.
b. If method = 2, computes sylvester2, Sylvester's matrix of 1853
of dimension (2*mx) x (2*mx). The determinants of properly chosen
submatrices of this matrix (a.k.a. ``modified'' subresultants) can be
used to compute the coefficients of the Sturmian PRS of f, g.
Applications of these Matrices can be found in the references below.
Especially, for applications of sylvester2, see the first reference!!
References
==========
1. Akritas, A. G., G.I. Malaschonok and P.S. Vigklas: ``On a Theorem
by Van Vleck Regarding Sturm Sequences. Serdica Journal of Computing,
Vol. 7, No 4, 101-134, 2013.
2. Akritas, A. G., G.I. Malaschonok and P.S. Vigklas: ``Sturm Sequences
and Modified Subresultant Polynomial Remainder Sequences.''
Serdica Journal of Computing, Vol. 8, No 1, 29-46, 2014.
'''
# obtain degrees of polys
m, n = degree( Poly(f, x), x), degree( Poly(g, x), x)
# Special cases:
# A:: case m = n < 0 (i.e. both polys are 0)
if m == n and n < 0:
return Matrix([])
# B:: case m = n = 0 (i.e. both polys are constants)
if m == n and n == 0:
return Matrix([])
# C:: m == 0 and n < 0 or m < 0 and n == 0
# (i.e. one poly is constant and the other is 0)
if m == 0 and n < 0:
return Matrix([])
elif m < 0 and n == 0:
return Matrix([])
# D:: m >= 1 and n < 0 or m < 0 and n >=1
# (i.e. one poly is of degree >=1 and the other is 0)
if m >= 1 and n < 0:
return Matrix([0])
elif m < 0 and n >= 1:
return Matrix([0])
fp = Poly(f, x).all_coeffs()
gp = Poly(g, x).all_coeffs()
# Sylvester's matrix of 1840 (default; a.k.a. sylvester1)
if method <= 1:
M = zeros(m + n)
k = 0
for i in range(n):
j = k
for coeff in fp:
M[i, j] = coeff
j = j + 1
k = k + 1
k = 0
for i in range(n, m + n):
j = k
for coeff in gp:
M[i, j] = coeff
j = j + 1
k = k + 1
return M
# Sylvester's matrix of 1853 (a.k.a sylvester2)
if method >= 2:
if len(fp) < len(gp):
h = []
for i in range(len(gp) - len(fp)):
h.append(0)
fp[ : 0] = h
else:
h = []
for i in range(len(fp) - len(gp)):
h.append(0)
gp[ : 0] = h
mx = max(m, n)
dim = 2*mx
M = zeros( dim )
k = 0
for i in range( mx ):
j = k
for coeff in fp:
M[2*i, j] = coeff
j = j + 1
j = k
for coeff in gp:
M[2*i + 1, j] = coeff
j = j + 1
k = k + 1
return M
def process_matrix_output(poly_seq, x):
"""
poly_seq is a polynomial remainder sequence computed either by
(modified_)subresultants_bezout or by (modified_)subresultants_sylv.
This function removes from poly_seq all zero polynomials as well
as all those whose degree is equal to the degree of a preceding
polynomial in poly_seq, as we scan it from left to right.
"""
L = poly_seq[:] # get a copy of the input sequence
d = degree(L[1], x)
i = 2
while i < len(L):
d_i = degree(L[i], x)
if d_i < 0: # zero poly
L.remove(L[i])
i = i - 1
if d == d_i: # poly degree equals degree of previous poly
L.remove(L[i])
i = i - 1
if d_i >= 0:
d = d_i
i = i + 1
return L
def subresultants_sylv(f, g, x):
"""
The input polynomials f, g are in Z[x] or in Q[x]. It is assumed
that deg(f) >= deg(g).
Computes the subresultant polynomial remainder sequence (prs)
of f, g by evaluating determinants of appropriately selected
submatrices of sylvester(f, g, x, 1). The dimensions of the
latter are (deg(f) + deg(g)) x (deg(f) + deg(g)).
Each coefficient is computed by evaluating the determinant of the
corresponding submatrix of sylvester(f, g, x, 1).
If the subresultant prs is complete, then the output coincides
with the Euclidean sequence of the polynomials f, g.
References:
===========
1. G.M.Diaz-Toca,L.Gonzalez-Vega: Various New Expressions for Subresultants
and Their Applications. Appl. Algebra in Engin., Communic. and Comp.,
Vol. 15, 233-266, 2004.
"""
# make sure neither f nor g is 0
if f == 0 or g == 0:
return [f, g]
n = degF = degree(f, x)
m = degG = degree(g, x)
# make sure proper degrees
if n == 0 and m == 0:
return [f, g]
if n < m:
n, m, degF, degG, f, g = m, n, degG, degF, g, f
if n > 0 and m == 0:
return [f, g]
SR_L = [f, g] # subresultant list
# form matrix sylvester(f, g, x, 1)
S = sylvester(f, g, x, 1)
# pick appropriate submatrices of S
# and form subresultant polys
j = m - 1
while j > 0:
Sp = S[:, :] # copy of S
# delete last j rows of coeffs of g
for ind in range(m + n - j, m + n):
Sp.row_del(m + n - j)
# delete last j rows of coeffs of f
for ind in range(m - j, m):
Sp.row_del(m - j)
# evaluate determinants and form coefficients list
coeff_L, k, l = [], Sp.rows, 0
while l <= j:
coeff_L.append(Sp[:, 0:k].det())
Sp.col_swap(k - 1, k + l)
l += 1
# form poly and append to SP_L
SR_L.append(Poly(coeff_L, x).as_expr())
j -= 1
# j = 0
SR_L.append(S.det())
return process_matrix_output(SR_L, x)
def modified_subresultants_sylv(f, g, x):
"""
The input polynomials f, g are in Z[x] or in Q[x]. It is assumed
that deg(f) >= deg(g).
Computes the modified subresultant polynomial remainder sequence (prs)
of f, g by evaluating determinants of appropriately selected
submatrices of sylvester(f, g, x, 2). The dimensions of the
latter are (2*deg(f)) x (2*deg(f)).
Each coefficient is computed by evaluating the determinant of the
corresponding submatrix of sylvester(f, g, x, 2).
If the modified subresultant prs is complete, then the output coincides
with the Sturmian sequence of the polynomials f, g.
References:
===========
1. A. G. Akritas,G.I. Malaschonok and P.S. Vigklas:
Sturm Sequences and Modified Subresultant Polynomial Remainder
Sequences. Serdica Journal of Computing, Vol. 8, No 1, 29--46, 2014.
"""
# make sure neither f nor g is 0
if f == 0 or g == 0:
return [f, g]
n = degF = degree(f, x)
m = degG = degree(g, x)
# make sure proper degrees
if n == 0 and m == 0:
return [f, g]
if n < m:
n, m, degF, degG, f, g = m, n, degG, degF, g, f
if n > 0 and m == 0:
return [f, g]
SR_L = [f, g] # modified subresultant list
# form matrix sylvester(f, g, x, 2)
S = sylvester(f, g, x, 2)
# pick appropriate submatrices of S
# and form modified subresultant polys
j = m - 1
while j > 0:
# delete last 2*j rows of pairs of coeffs of f, g
Sp = S[0:2*n - 2*j, :] # copy of first 2*n - 2*j rows of S
# evaluate determinants and form coefficients list
coeff_L, k, l = [], Sp.rows, 0
while l <= j:
coeff_L.append(Sp[:, 0:k].det())
Sp.col_swap(k - 1, k + l)
l += 1
# form poly and append to SP_L
SR_L.append(Poly(coeff_L, x).as_expr())
j -= 1
# j = 0
SR_L.append(S.det())
return process_matrix_output(SR_L, x)
def res(f, g, x):
"""
The input polynomials f, g are in Z[x] or in Q[x].
The output is the resultant of f, g computed by evaluating
the determinant of the matrix sylvester(f, g, x, 1).
References:
===========
1. J. S. Cohen: Computer Algebra and Symbolic Computation
- Mathematical Methods. A. K. Peters, 2003.
"""
if f == 0 or g == 0:
raise PolynomialError("The resultant of %s and %s is not defined" % (f, g))
else:
return sylvester(f, g, x, 1).det()
def res_q(f, g, x):
"""
The input polynomials f, g are in Z[x] or in Q[x].
The output is the resultant of f, g computed recursively
by polynomial divisions in Q[x], using the function rem.
See Cohen's book p. 281.
References:
===========
1. J. S. Cohen: Computer Algebra and Symbolic Computation
- Mathematical Methods. A. K. Peters, 2003.
"""
m = degree(f, x)
n = degree(g, x)
if m < n:
return (-1)**(m*n) * res_q(g, f, x)
elif n == 0: # g is a constant
return g**m
else:
r = rem(f, g, x)
if r == 0:
return 0
else:
s = degree(r, x)
l = LC(g, x)
return (-1)**(m*n) * l**(m-s)*res_q(g, r, x)
def res_z(f, g, x):
"""
The input polynomials f, g are in Z[x] or in Q[x].
The output is the resultant of f, g computed recursively
by polynomial divisions in Z[x], using the function prem().
See Cohen's book p. 283.
References:
===========
1. J. S. Cohen: Computer Algebra and Symbolic Computation
- Mathematical Methods. A. K. Peters, 2003.
"""
m = degree(f, x)
n = degree(g, x)
if m < n:
return (-1)**(m*n) * res_z(g, f, x)
elif n == 0: # g is a constant
return g**m
else:
r = prem(f, g, x)
if r == 0:
return 0
else:
delta = m - n + 1
w = (-1)**(m*n) * res_z(g, r, x)
s = degree(r, x)
l = LC(g, x)
k = delta * n - m + s
return quo(w, l**k, x)
def sign_seq(poly_seq, x):
"""
Given a sequence of polynomials poly_seq, it returns
the sequence of signs of the leading coefficients of
the polynomials in poly_seq.
"""
return [sign(LC(poly_seq[i], x)) for i in range(len(poly_seq))]
def bezout(p, q, x, method='bz'):
"""
The input polynomials p, q are in Z[x] or in Q[x]. Let
mx = max(degree(p, x), degree(q, x)).
The default option bezout(p, q, x, method='bz') returns Bezout's
symmetric matrix of p and q, of dimensions (mx) x (mx). The
determinant of this matrix is equal to the determinant of sylvester2,
Sylvester's matrix of 1853, whose dimensions are (2*mx) x (2*mx);
however the subresultants of these two matrices may differ.
The other option, bezout(p, q, x, 'prs'), is of interest to us
in this module because it returns a matrix equivalent to sylvester2.
In this case all subresultants of the two matrices are identical.
Both the subresultant polynomial remainder sequence (prs) and
the modified subresultant prs of p and q can be computed by
evaluating determinants of appropriately selected submatrices of
bezout(p, q, x, 'prs') --- one determinant per coefficient of the
remainder polynomials.
The matrices bezout(p, q, x, 'bz') and bezout(p, q, x, 'prs')
are related by the formula
bezout(p, q, x, 'prs') =
backward_eye(deg(p)) * bezout(p, q, x, 'bz') * backward_eye(deg(p)),
where backward_eye() is the backward identity function.
References
==========
1. G.M.Diaz-Toca,L.Gonzalez-Vega: Various New Expressions for Subresultants
and Their Applications. Appl. Algebra in Engin., Communic. and Comp.,
Vol. 15, 233-266, 2004.
"""
# obtain degrees of polys
m, n = degree( Poly(p, x), x), degree( Poly(q, x), x)
# Special cases:
# A:: case m = n < 0 (i.e. both polys are 0)
if m == n and n < 0:
return Matrix([])
# B:: case m = n = 0 (i.e. both polys are constants)
if m == n and n == 0:
return Matrix([])
# C:: m == 0 and n < 0 or m < 0 and n == 0
# (i.e. one poly is constant and the other is 0)
if m == 0 and n < 0:
return Matrix([])
elif m < 0 and n == 0:
return Matrix([])
# D:: m >= 1 and n < 0 or m < 0 and n >=1
# (i.e. one poly is of degree >=1 and the other is 0)
if m >= 1 and n < 0:
return Matrix([0])
elif m < 0 and n >= 1:
return Matrix([0])
y = var('y')
# expr is 0 when x = y
expr = p * q.subs({x:y}) - p.subs({x:y}) * q
# hence expr is exactly divisible by x - y
poly = Poly( quo(expr, x-y), x, y)
# form Bezout matrix and store them in B as indicated to get
# the LC coefficient of each poly either in the first position
# of each row (method='prs') or in the last (method='bz').
mx = max(m, n)
B = zeros(mx)
for i in range(mx):
for j in range(mx):
if method == 'prs':
B[mx - 1 - i, mx - 1 - j] = poly.nth(i, j)
else:
B[i, j] = poly.nth(i, j)
return B
def backward_eye(n):
'''
Returns the backward identity matrix of dimensions n x n.
Needed to "turn" the Bezout matrices
so that the leading coefficients are first.
See docstring of the function bezout(p, q, x, method='bz').
'''
M = eye(n) # identity matrix of order n
for i in range(int(M.rows / 2)):
M.row_swap(0 + i, M.rows - 1 - i)
return M
def subresultants_bezout(p, q, x):
"""
The input polynomials p, q are in Z[x] or in Q[x]. It is assumed
that degree(p, x) >= degree(q, x).
Computes the subresultant polynomial remainder sequence
of p, q by evaluating determinants of appropriately selected
submatrices of bezout(p, q, x, 'prs'). The dimensions of the
latter are deg(p) x deg(p).
Each coefficient is computed by evaluating the determinant of the
corresponding submatrix of bezout(p, q, x, 'prs').
bezout(p, q, x, 'prs) is used instead of sylvester(p, q, x, 1),
Sylvester's matrix of 1840, because the dimensions of the latter
are (deg(p) + deg(q)) x (deg(p) + deg(q)).
If the subresultant prs is complete, then the output coincides
with the Euclidean sequence of the polynomials p, q.
References
==========
1. G.M.Diaz-Toca,L.Gonzalez-Vega: Various New Expressions for Subresultants
and Their Applications. Appl. Algebra in Engin., Communic. and Comp.,
Vol. 15, 233-266, 2004.
"""
# make sure neither p nor q is 0
if p == 0 or q == 0:
return [p, q]
f, g = p, q
n = degF = degree(f, x)
m = degG = degree(g, x)
# make sure proper degrees
if n == 0 and m == 0:
return [f, g]
if n < m:
n, m, degF, degG, f, g = m, n, degG, degF, g, f
if n > 0 and m == 0:
return [f, g]
SR_L = [f, g] # subresultant list
F = LC(f, x)**(degF - degG)
# form the bezout matrix
B = bezout(f, g, x, 'prs')
# pick appropriate submatrices of B
# and form subresultant polys
if degF > degG:
j = 2
if degF == degG:
j = 1
while j <= degF:
M = B[0:j, :]
k, coeff_L = j - 1, []
while k <= degF - 1:
coeff_L.append(M[:, 0:j].det())
if k < degF - 1:
M.col_swap(j - 1, k + 1)
k = k + 1
# apply Theorem 2.1 in the paper by Toca & Vega 2004
# to get correct signs
SR_L.append(int((-1)**(j*(j-1)/2)) * (Poly(coeff_L, x) / F).as_expr())
j = j + 1
return process_matrix_output(SR_L, x)
def modified_subresultants_bezout(p, q, x):
"""
The input polynomials p, q are in Z[x] or in Q[x]. It is assumed
that degree(p, x) >= degree(q, x).
Computes the modified subresultant polynomial remainder sequence
of p, q by evaluating determinants of appropriately selected
submatrices of bezout(p, q, x, 'prs'). The dimensions of the
latter are deg(p) x deg(p).
Each coefficient is computed by evaluating the determinant of the
corresponding submatrix of bezout(p, q, x, 'prs').
bezout(p, q, x, 'prs') is used instead of sylvester(p, q, x, 2),
Sylvester's matrix of 1853, because the dimensions of the latter
are 2*deg(p) x 2*deg(p).
If the modified subresultant prs is complete, and LC( p ) > 0, the output
coincides with the (generalized) Sturm's sequence of the polynomials p, q.
References
==========
1. Akritas, A. G., G.I. Malaschonok and P.S. Vigklas: ``Sturm Sequences
and Modified Subresultant Polynomial Remainder Sequences.''
Serdica Journal of Computing, Vol. 8, No 1, 29-46, 2014.
2. G.M.Diaz-Toca,L.Gonzalez-Vega: Various New Expressions for Subresultants
and Their Applications. Appl. Algebra in Engin., Communic. and Comp.,
Vol. 15, 233-266, 2004.
"""
# make sure neither p nor q is 0
if p == 0 or q == 0:
return [p, q]
f, g = p, q
n = degF = degree(f, x)
m = degG = degree(g, x)
# make sure proper degrees
if n == 0 and m == 0:
return [f, g]
if n < m:
n, m, degF, degG, f, g = m, n, degG, degF, g, f
if n > 0 and m == 0:
return [f, g]
SR_L = [f, g] # subresultant list
# form the bezout matrix
B = bezout(f, g, x, 'prs')
# pick appropriate submatrices of B
# and form subresultant polys
if degF > degG:
j = 2
if degF == degG:
j = 1
while j <= degF:
M = B[0:j, :]
k, coeff_L = j - 1, []
while k <= degF - 1:
coeff_L.append(M[:, 0:j].det())
if k < degF - 1:
M.col_swap(j - 1, k + 1)
k = k + 1
## Theorem 2.1 in the paper by Toca & Vega 2004 is _not needed_
## in this case since
## the bezout matrix is equivalent to sylvester2
SR_L.append(( Poly(coeff_L, x)).as_expr())
j = j + 1
return process_matrix_output(SR_L, x)
def sturm_pg(p, q, x, method=0):
"""
p, q are polynomials in Z[x] or Q[x]. It is assumed
that degree(p, x) >= degree(q, x).
Computes the (generalized) Sturm sequence of p and q in Z[x] or Q[x].
If q = diff(p, x, 1) it is the usual Sturm sequence.
A. If method == 0, default, the remainder coefficients of the sequence
are (in absolute value) ``modified'' subresultants, which for non-monic
polynomials are greater than the coefficients of the corresponding
subresultants by the factor Abs(LC(p)**( deg(p)- deg(q))).
B. If method == 1, the remainder coefficients of the sequence are (in
absolute value) subresultants, which for non-monic polynomials are
smaller than the coefficients of the corresponding ``modified''
subresultants by the factor Abs(LC(p)**( deg(p)- deg(q))).
If the Sturm sequence is complete, method=0 and LC( p ) > 0, the coefficients
of the polynomials in the sequence are ``modified'' subresultants.
That is, they are determinants of appropriately selected submatrices of
sylvester2, Sylvester's matrix of 1853. In this case the Sturm sequence
coincides with the ``modified'' subresultant prs, of the polynomials
p, q.
If the Sturm sequence is incomplete and method=0 then the signs of the
coefficients of the polynomials in the sequence may differ from the signs
of the coefficients of the corresponding polynomials in the ``modified''
subresultant prs; however, the absolute values are the same.
To compute the coefficients, no determinant evaluation takes place. Instead,
polynomial divisions in Q[x] are performed, using the function rem(p, q, x);
the coefficients of the remainders computed this way become (``modified'')
subresultants with the help of the Pell-Gordon Theorem of 1917.
See also the function euclid_pg(p, q, x).
References
==========
1. Pell A. J., R. L. Gordon. The Modified Remainders Obtained in Finding
the Highest Common Factor of Two Polynomials. Annals of MatheMatics,
Second Series, 18 (1917), No. 4, 188-193.
2. Akritas, A. G., G.I. Malaschonok and P.S. Vigklas: ``Sturm Sequences
and Modified Subresultant Polynomial Remainder Sequences.''
Serdica Journal of Computing, Vol. 8, No 1, 29-46, 2014.
"""
# make sure neither p nor q is 0
if p == 0 or q == 0:
return [p, q]
# make sure proper degrees
d0 = degree(p, x)
d1 = degree(q, x)
if d0 == 0 and d1 == 0:
return [p, q]
if d1 > d0:
d0, d1 = d1, d0
p, q = q, p
if d0 > 0 and d1 == 0:
return [p,q]
# make sure LC(p) > 0
flag = 0
if LC(p,x) < 0:
flag = 1
p = -p
q = -q
# initialize
lcf = LC(p, x)**(d0 - d1) # lcf * subr = modified subr
a0, a1 = p, q # the input polys
sturm_seq = [a0, a1] # the output list
del0 = d0 - d1 # degree difference
rho1 = LC(a1, x) # leading coeff of a1
exp_deg = d1 - 1 # expected degree of a2
a2 = - rem(a0, a1, domain=QQ) # first remainder
rho2 = LC(a2,x) # leading coeff of a2
d2 = degree(a2, x) # actual degree of a2
deg_diff_new = exp_deg - d2 # expected - actual degree
del1 = d1 - d2 # degree difference
# mul_fac is the factor by which a2 is multiplied to
# get integer coefficients
mul_fac_old = rho1**(del0 + del1 - deg_diff_new)
# append accordingly
if method == 0:
sturm_seq.append( simplify(lcf * a2 * Abs(mul_fac_old)))
else:
sturm_seq.append( simplify( a2 * Abs(mul_fac_old)))
# main loop
deg_diff_old = deg_diff_new
while d2 > 0:
a0, a1, d0, d1 = a1, a2, d1, d2 # update polys and degrees
del0 = del1 # update degree difference
exp_deg = d1 - 1 # new expected degree
a2 = - rem(a0, a1, domain=QQ) # new remainder
rho3 = LC(a2, x) # leading coeff of a2
d2 = degree(a2, x) # actual degree of a2
deg_diff_new = exp_deg - d2 # expected - actual degree
del1 = d1 - d2 # degree difference
# take into consideration the power
# rho1**deg_diff_old that was "left out"
expo_old = deg_diff_old # rho1 raised to this power
expo_new = del0 + del1 - deg_diff_new # rho2 raised to this power
# update variables and append
mul_fac_new = rho2**(expo_new) * rho1**(expo_old) * mul_fac_old
deg_diff_old, mul_fac_old = deg_diff_new, mul_fac_new
rho1, rho2 = rho2, rho3
if method == 0:
sturm_seq.append( simplify(lcf * a2 * Abs(mul_fac_old)))
else:
sturm_seq.append( simplify( a2 * Abs(mul_fac_old)))
if flag: # change the sign of the sequence
sturm_seq = [-i for i in sturm_seq]
# gcd is of degree > 0 ?
m = len(sturm_seq)
if sturm_seq[m - 1] == nan or sturm_seq[m - 1] == 0:
sturm_seq.pop(m - 1)
return sturm_seq
def sturm_q(p, q, x):
"""
p, q are polynomials in Z[x] or Q[x]. It is assumed
that degree(p, x) >= degree(q, x).
Computes the (generalized) Sturm sequence of p and q in Q[x].
Polynomial divisions in Q[x] are performed, using the function rem(p, q, x).
The coefficients of the polynomials in the Sturm sequence can be uniquely
determined from the corresponding coefficients of the polynomials found
either in:
(a) the ``modified'' subresultant prs, (references 1, 2)
or in
(b) the subresultant prs (reference 3).
References
==========
1. Pell A. J., R. L. Gordon. The Modified Remainders Obtained in Finding
the Highest Common Factor of Two Polynomials. Annals of MatheMatics,
Second Series, 18 (1917), No. 4, 188-193.
2 Akritas, A. G., G.I. Malaschonok and P.S. Vigklas: ``Sturm Sequences
and Modified Subresultant Polynomial Remainder Sequences.''
Serdica Journal of Computing, Vol. 8, No 1, 29-46, 2014.
3. Akritas, A. G., G.I. Malaschonok and P.S. Vigklas: ``A Basic Result
on the Theory of Subresultants.'' Serdica Journal of Computing 10 (2016), No.1, 31-48.
"""
# make sure neither p nor q is 0
if p == 0 or q == 0:
return [p, q]
# make sure proper degrees
d0 = degree(p, x)
d1 = degree(q, x)
if d0 == 0 and d1 == 0:
return [p, q]
if d1 > d0:
d0, d1 = d1, d0
p, q = q, p
if d0 > 0 and d1 == 0:
return [p,q]
# make sure LC(p) > 0
flag = 0
if LC(p,x) < 0:
flag = 1
p = -p
q = -q
# initialize
a0, a1 = p, q # the input polys
sturm_seq = [a0, a1] # the output list
a2 = -rem(a0, a1, domain=QQ) # first remainder
d2 = degree(a2, x) # degree of a2
sturm_seq.append( a2 )
# main loop
while d2 > 0:
a0, a1, d0, d1 = a1, a2, d1, d2 # update polys and degrees
a2 = -rem(a0, a1, domain=QQ) # new remainder
d2 = degree(a2, x) # actual degree of a2
sturm_seq.append( a2 )
if flag: # change the sign of the sequence
sturm_seq = [-i for i in sturm_seq]
# gcd is of degree > 0 ?
m = len(sturm_seq)
if sturm_seq[m - 1] == nan or sturm_seq[m - 1] == 0:
sturm_seq.pop(m - 1)
return sturm_seq
def sturm_amv(p, q, x, method=0):
"""
p, q are polynomials in Z[x] or Q[x]. It is assumed
that degree(p, x) >= degree(q, x).
Computes the (generalized) Sturm sequence of p and q in Z[x] or Q[x].
If q = diff(p, x, 1) it is the usual Sturm sequence.
A. If method == 0, default, the remainder coefficients of the
sequence are (in absolute value) ``modified'' subresultants, which
for non-monic polynomials are greater than the coefficients of the
corresponding subresultants by the factor Abs(LC(p)**( deg(p)- deg(q))).
B. If method == 1, the remainder coefficients of the sequence are (in
absolute value) subresultants, which for non-monic polynomials are
smaller than the coefficients of the corresponding ``modified''
subresultants by the factor Abs( LC(p)**( deg(p)- deg(q)) ).
If the Sturm sequence is complete, method=0 and LC( p ) > 0, then the
coefficients of the polynomials in the sequence are ``modified'' subresultants.
That is, they are determinants of appropriately selected submatrices of
sylvester2, Sylvester's matrix of 1853. In this case the Sturm sequence
coincides with the ``modified'' subresultant prs, of the polynomials
p, q.
If the Sturm sequence is incomplete and method=0 then the signs of the
coefficients of the polynomials in the sequence may differ from the signs
of the coefficients of the corresponding polynomials in the ``modified''
subresultant prs; however, the absolute values are the same.
To compute the coefficients, no determinant evaluation takes place.
Instead, we first compute the euclidean sequence of p and q using
euclid_amv(p, q, x) and then: (a) change the signs of the remainders in the
Euclidean sequence according to the pattern "-, -, +, +, -, -, +, +,..."
(see Lemma 1 in the 1st reference or Theorem 3 in the 2nd reference)
and (b) if method=0, assuming deg(p) > deg(q), we multiply the remainder
coefficients of the Euclidean sequence times the factor
Abs( LC(p)**( deg(p)- deg(q)) ) to make them modified subresultants.
See also the function sturm_pg(p, q, x).
References
==========
1. Akritas, A. G., G.I. Malaschonok and P.S. Vigklas: ``A Basic Result
on the Theory of Subresultants.'' Serdica Journal of Computing 10 (2016), No.1, 31-48.
2. Akritas, A. G., G.I. Malaschonok and P.S. Vigklas: ``On the Remainders
Obtained in Finding the Greatest Common Divisor of Two Polynomials.'' Serdica
Journal of Computing 9(2) (2015), 123-138.
3. Akritas, A. G., G.I. Malaschonok and P.S. Vigklas: ``Subresultant Polynomial
Remainder Sequences Obtained by Polynomial Divisions in Q[x] or in Z[x].''
Serdica Journal of Computing 10 (2016), No.3-4, 197-217.
"""
# compute the euclidean sequence
prs = euclid_amv(p, q, x)
# defensive
if prs == [] or len(prs) == 2:
return prs
# the coefficients in prs are subresultants and hence are smaller
# than the corresponding subresultants by the factor
# Abs( LC(prs[0])**( deg(prs[0]) - deg(prs[1])) ); Theorem 2, 2nd reference.
lcf = Abs( LC(prs[0])**( degree(prs[0], x) - degree(prs[1], x) ) )
# the signs of the first two polys in the sequence stay the same
sturm_seq = [prs[0], prs[1]]
# change the signs according to "-, -, +, +, -, -, +, +,..."
# and multiply times lcf if needed
flag = 0
m = len(prs)
i = 2
while i <= m-1:
if flag == 0:
sturm_seq.append( - prs[i] )
i = i + 1
if i == m:
break
sturm_seq.append( - prs[i] )
i = i + 1
flag = 1
elif flag == 1:
sturm_seq.append( prs[i] )
i = i + 1
if i == m:
break
sturm_seq.append( prs[i] )
i = i + 1
flag = 0
# subresultants or modified subresultants?
if method == 0 and lcf > 1:
aux_seq = [sturm_seq[0], sturm_seq[1]]
for i in range(2, m):
aux_seq.append(simplify(sturm_seq[i] * lcf ))
sturm_seq = aux_seq
return sturm_seq
def euclid_pg(p, q, x):
"""
p, q are polynomials in Z[x] or Q[x]. It is assumed
that degree(p, x) >= degree(q, x).
Computes the Euclidean sequence of p and q in Z[x] or Q[x].
If the Euclidean sequence is complete the coefficients of the polynomials
in the sequence are subresultants. That is, they are determinants of
appropriately selected submatrices of sylvester1, Sylvester's matrix of 1840.
In this case the Euclidean sequence coincides with the subresultant prs
of the polynomials p, q.
If the Euclidean sequence is incomplete the signs of the coefficients of the
polynomials in the sequence may differ from the signs of the coefficients of
the corresponding polynomials in the subresultant prs; however, the absolute
values are the same.
To compute the Euclidean sequence, no determinant evaluation takes place.
We first compute the (generalized) Sturm sequence of p and q using
sturm_pg(p, q, x, 1), in which case the coefficients are (in absolute value)
equal to subresultants. Then we change the signs of the remainders in the
Sturm sequence according to the pattern "-, -, +, +, -, -, +, +,..." ;
see Lemma 1 in the 1st reference or Theorem 3 in the 2nd reference as well as
the function sturm_pg(p, q, x).
References
==========
1. Akritas, A. G., G.I. Malaschonok and P.S. Vigklas: ``A Basic Result
on the Theory of Subresultants.'' Serdica Journal of Computing 10 (2016), No.1, 31-48.
2. Akritas, A. G., G.I. Malaschonok and P.S. Vigklas: ``On the Remainders
Obtained in Finding the Greatest Common Divisor of Two Polynomials.'' Serdica
Journal of Computing 9(2) (2015), 123-138.
3. Akritas, A. G., G.I. Malaschonok and P.S. Vigklas: ``Subresultant Polynomial
Remainder Sequences Obtained by Polynomial Divisions in Q[x] or in Z[x].''
Serdica Journal of Computing 10 (2016), No.3-4, 197-217.
"""
# compute the sturmian sequence using the Pell-Gordon (or AMV) theorem
# with the coefficients in the prs being (in absolute value) subresultants
prs = sturm_pg(p, q, x, 1) ## any other method would do
# defensive
if prs == [] or len(prs) == 2:
return prs
# the signs of the first two polys in the sequence stay the same
euclid_seq = [prs[0], prs[1]]
# change the signs according to "-, -, +, +, -, -, +, +,..."
flag = 0
m = len(prs)
i = 2
while i <= m-1:
if flag == 0:
euclid_seq.append(- prs[i] )
i = i + 1
if i == m:
break
euclid_seq.append(- prs[i] )
i = i + 1
flag = 1
elif flag == 1:
euclid_seq.append(prs[i] )
i = i + 1
if i == m:
break
euclid_seq.append(prs[i] )
i = i + 1
flag = 0
return euclid_seq
def euclid_q(p, q, x):
"""
p, q are polynomials in Z[x] or Q[x]. It is assumed
that degree(p, x) >= degree(q, x).
Computes the Euclidean sequence of p and q in Q[x].
Polynomial divisions in Q[x] are performed, using the function rem(p, q, x).
The coefficients of the polynomials in the Euclidean sequence can be uniquely
determined from the corresponding coefficients of the polynomials found
either in:
(a) the ``modified'' subresultant polynomial remainder sequence,
(references 1, 2)
or in
(b) the subresultant polynomial remainder sequence (references 3).
References
==========
1. Pell A. J., R. L. Gordon. The Modified Remainders Obtained in Finding
the Highest Common Factor of Two Polynomials. Annals of MatheMatics,
Second Series, 18 (1917), No. 4, 188-193.
2. Akritas, A. G., G.I. Malaschonok and P.S. Vigklas: ``Sturm Sequences
and Modified Subresultant Polynomial Remainder Sequences.''
Serdica Journal of Computing, Vol. 8, No 1, 29-46, 2014.
3. Akritas, A. G., G.I. Malaschonok and P.S. Vigklas: ``A Basic Result
on the Theory of Subresultants.'' Serdica Journal of Computing 10 (2016), No.1, 31-48.
"""
# make sure neither p nor q is 0
if p == 0 or q == 0:
return [p, q]
# make sure proper degrees
d0 = degree(p, x)
d1 = degree(q, x)
if d0 == 0 and d1 == 0:
return [p, q]
if d1 > d0:
d0, d1 = d1, d0
p, q = q, p
if d0 > 0 and d1 == 0:
return [p,q]
# make sure LC(p) > 0
flag = 0
if LC(p,x) < 0:
flag = 1
p = -p
q = -q
# initialize
a0, a1 = p, q # the input polys
euclid_seq = [a0, a1] # the output list
a2 = rem(a0, a1, domain=QQ) # first remainder
d2 = degree(a2, x) # degree of a2
euclid_seq.append( a2 )
# main loop
while d2 > 0:
a0, a1, d0, d1 = a1, a2, d1, d2 # update polys and degrees
a2 = rem(a0, a1, domain=QQ) # new remainder
d2 = degree(a2, x) # actual degree of a2
euclid_seq.append( a2 )
if flag: # change the sign of the sequence
euclid_seq = [-i for i in euclid_seq]
# gcd is of degree > 0 ?
m = len(euclid_seq)
if euclid_seq[m - 1] == nan or euclid_seq[m - 1] == 0:
euclid_seq.pop(m - 1)
return euclid_seq
def euclid_amv(f, g, x):
"""
f, g are polynomials in Z[x] or Q[x]. It is assumed
that degree(f, x) >= degree(g, x).
Computes the Euclidean sequence of p and q in Z[x] or Q[x].
If the Euclidean sequence is complete the coefficients of the polynomials
in the sequence are subresultants. That is, they are determinants of
appropriately selected submatrices of sylvester1, Sylvester's matrix of 1840.
In this case the Euclidean sequence coincides with the subresultant prs,
of the polynomials p, q.
If the Euclidean sequence is incomplete the signs of the coefficients of the
polynomials in the sequence may differ from the signs of the coefficients of
the corresponding polynomials in the subresultant prs; however, the absolute
values are the same.
To compute the coefficients, no determinant evaluation takes place.
Instead, polynomial divisions in Z[x] or Q[x] are performed, using
the function rem_z(f, g, x); the coefficients of the remainders
computed this way become subresultants with the help of the
Collins-Brown-Traub formula for coefficient reduction.
References
==========
1. Akritas, A. G., G.I. Malaschonok and P.S. Vigklas: ``A Basic Result
on the Theory of Subresultants.'' Serdica Journal of Computing 10 (2016), No.1, 31-48.
2. Akritas, A. G., G.I. Malaschonok and P.S. Vigklas: ``Subresultant Polynomial
remainder Sequences Obtained by Polynomial Divisions in Q[x] or in Z[x].''
Serdica Journal of Computing 10 (2016), No.3-4, 197-217.
"""
# make sure neither f nor g is 0
if f == 0 or g == 0:
return [f, g]
# make sure proper degrees
d0 = degree(f, x)
d1 = degree(g, x)
if d0 == 0 and d1 == 0:
return [f, g]
if d1 > d0:
d0, d1 = d1, d0
f, g = g, f
if d0 > 0 and d1 == 0:
return [f, g]
# initialize
a0 = f
a1 = g
euclid_seq = [a0, a1]
deg_dif_p1, c = degree(a0, x) - degree(a1, x) + 1, -1
# compute the first polynomial of the prs
i = 1
a2 = rem_z(a0, a1, x) / Abs( (-1)**deg_dif_p1 ) # first remainder
euclid_seq.append( a2 )
d2 = degree(a2, x) # actual degree of a2
# main loop
while d2 >= 1:
a0, a1, d0, d1 = a1, a2, d1, d2 # update polys and degrees
i += 1
sigma0 = -LC(a0)
c = (sigma0**(deg_dif_p1 - 1)) / (c**(deg_dif_p1 - 2))
deg_dif_p1 = degree(a0, x) - d2 + 1
a2 = rem_z(a0, a1, x) / Abs( (c**(deg_dif_p1 - 1)) * sigma0 )
euclid_seq.append( a2 )
d2 = degree(a2, x) # actual degree of a2
# gcd is of degree > 0 ?
m = len(euclid_seq)
if euclid_seq[m - 1] == nan or euclid_seq[m - 1] == 0:
euclid_seq.pop(m - 1)
return euclid_seq
def modified_subresultants_pg(p, q, x):
"""
p, q are polynomials in Z[x] or Q[x]. It is assumed
that degree(p, x) >= degree(q, x).
Computes the ``modified'' subresultant prs of p and q in Z[x] or Q[x];
the coefficients of the polynomials in the sequence are
``modified'' subresultants. That is, they are determinants of appropriately
selected submatrices of sylvester2, Sylvester's matrix of 1853.
To compute the coefficients, no determinant evaluation takes place. Instead,
polynomial divisions in Q[x] are performed, using the function rem(p, q, x);
the coefficients of the remainders computed this way become ``modified''
subresultants with the help of the Pell-Gordon Theorem of 1917.
If the ``modified'' subresultant prs is complete, and LC( p ) > 0, it coincides
with the (generalized) Sturm sequence of the polynomials p, q.
References
==========
1. Pell A. J., R. L. Gordon. The Modified Remainders Obtained in Finding
the Highest Common Factor of Two Polynomials. Annals of MatheMatics,
Second Series, 18 (1917), No. 4, 188-193.
2. Akritas, A. G., G.I. Malaschonok and P.S. Vigklas: ``Sturm Sequences
and Modified Subresultant Polynomial Remainder Sequences.''
Serdica Journal of Computing, Vol. 8, No 1, 29-46, 2014.
"""
# make sure neither p nor q is 0
if p == 0 or q == 0:
return [p, q]
# make sure proper degrees
d0 = degree(p,x)
d1 = degree(q,x)
if d0 == 0 and d1 == 0:
return [p, q]
if d1 > d0:
d0, d1 = d1, d0
p, q = q, p
if d0 > 0 and d1 == 0:
return [p,q]
# initialize
k = var('k') # index in summation formula
u_list = [] # of elements (-1)**u_i
subres_l = [p, q] # mod. subr. prs output list
a0, a1 = p, q # the input polys
del0 = d0 - d1 # degree difference
degdif = del0 # save it
rho_1 = LC(a0) # lead. coeff (a0)
# Initialize Pell-Gordon variables
rho_list_minus_1 = sign( LC(a0, x)) # sign of LC(a0)
rho1 = LC(a1, x) # leading coeff of a1
rho_list = [ sign(rho1)] # of signs
p_list = [del0] # of degree differences
u = summation(k, (k, 1, p_list[0])) # value of u
u_list.append(u) # of u values
v = sum(p_list) # v value
# first remainder
exp_deg = d1 - 1 # expected degree of a2
a2 = - rem(a0, a1, domain=QQ) # first remainder
rho2 = LC(a2, x) # leading coeff of a2
d2 = degree(a2, x) # actual degree of a2
deg_diff_new = exp_deg - d2 # expected - actual degree
del1 = d1 - d2 # degree difference
# mul_fac is the factor by which a2 is multiplied to
# get integer coefficients
mul_fac_old = rho1**(del0 + del1 - deg_diff_new)
# update Pell-Gordon variables
p_list.append(1 + deg_diff_new) # deg_diff_new is 0 for complete seq
# apply Pell-Gordon formula (7) in second reference
num = 1 # numerator of fraction
for u in u_list:
num *= (-1)**u
num = num * (-1)**v
# denominator depends on complete / incomplete seq
if deg_diff_new == 0: # complete seq
den = 1
for k in range(len(rho_list)):
den *= rho_list[k]**(p_list[k] + p_list[k + 1])
den = den * rho_list_minus_1
else: # incomplete seq
den = 1
for k in range(len(rho_list)-1):
den *= rho_list[k]**(p_list[k] + p_list[k + 1])
den = den * rho_list_minus_1
expo = (p_list[len(rho_list) - 1] + p_list[len(rho_list)] - deg_diff_new)
den = den * rho_list[len(rho_list) - 1]**expo
# the sign of the determinant depends on sg(num / den)
if sign(num / den) > 0:
subres_l.append( simplify(rho_1**degdif*a2* Abs(mul_fac_old) ) )
else:
subres_l.append(- simplify(rho_1**degdif*a2* Abs(mul_fac_old) ) )
# update Pell-Gordon variables
k = var('k')
rho_list.append( sign(rho2))
u = summation(k, (k, 1, p_list[len(p_list) - 1]))
u_list.append(u)
v = sum(p_list)
deg_diff_old=deg_diff_new
# main loop
while d2 > 0:
a0, a1, d0, d1 = a1, a2, d1, d2 # update polys and degrees
del0 = del1 # update degree difference
exp_deg = d1 - 1 # new expected degree
a2 = - rem(a0, a1, domain=QQ) # new remainder
rho3 = LC(a2, x) # leading coeff of a2
d2 = degree(a2, x) # actual degree of a2
deg_diff_new = exp_deg - d2 # expected - actual degree
del1 = d1 - d2 # degree difference
# take into consideration the power
# rho1**deg_diff_old that was "left out"
expo_old = deg_diff_old # rho1 raised to this power
expo_new = del0 + del1 - deg_diff_new # rho2 raised to this power
mul_fac_new = rho2**(expo_new) * rho1**(expo_old) * mul_fac_old
# update variables
deg_diff_old, mul_fac_old = deg_diff_new, mul_fac_new
rho1, rho2 = rho2, rho3
# update Pell-Gordon variables
p_list.append(1 + deg_diff_new) # deg_diff_new is 0 for complete seq
# apply Pell-Gordon formula (7) in second reference
num = 1 # numerator
for u in u_list:
num *= (-1)**u
num = num * (-1)**v
# denominator depends on complete / incomplete seq
if deg_diff_new == 0: # complete seq
den = 1
for k in range(len(rho_list)):
den *= rho_list[k]**(p_list[k] + p_list[k + 1])
den = den * rho_list_minus_1
else: # incomplete seq
den = 1
for k in range(len(rho_list)-1):
den *= rho_list[k]**(p_list[k] + p_list[k + 1])
den = den * rho_list_minus_1
expo = (p_list[len(rho_list) - 1] + p_list[len(rho_list)] - deg_diff_new)
den = den * rho_list[len(rho_list) - 1]**expo
# the sign of the determinant depends on sg(num / den)
if sign(num / den) > 0:
subres_l.append( simplify(rho_1**degdif*a2* Abs(mul_fac_old) ) )
else:
subres_l.append(- simplify(rho_1**degdif*a2* Abs(mul_fac_old) ) )
# update Pell-Gordon variables
k = var('k')
rho_list.append( sign(rho2))
u = summation(k, (k, 1, p_list[len(p_list) - 1]))
u_list.append(u)
v = sum(p_list)
# gcd is of degree > 0 ?
m = len(subres_l)
if subres_l[m - 1] == nan or subres_l[m - 1] == 0:
subres_l.pop(m - 1)
# LC( p ) < 0
m = len(subres_l) # list may be shorter now due to deg(gcd ) > 0
if LC( p ) < 0:
aux_seq = [subres_l[0], subres_l[1]]
for i in range(2, m):
aux_seq.append(simplify(subres_l[i] * (-1) ))
subres_l = aux_seq
return subres_l
def subresultants_pg(p, q, x):
"""
p, q are polynomials in Z[x] or Q[x]. It is assumed
that degree(p, x) >= degree(q, x).
Computes the subresultant prs of p and q in Z[x] or Q[x], from
the modified subresultant prs of p and q.
The coefficients of the polynomials in these two sequences differ only
in sign and the factor LC(p)**( deg(p)- deg(q)) as stated in
Theorem 2 of the reference.
The coefficients of the polynomials in the output sequence are
subresultants. That is, they are determinants of appropriately
selected submatrices of sylvester1, Sylvester's matrix of 1840.
If the subresultant prs is complete, then it coincides with the
Euclidean sequence of the polynomials p, q.
References
==========
1. Akritas, A. G., G.I. Malaschonok and P.S. Vigklas: "On the Remainders
Obtained in Finding the Greatest Common Divisor of Two Polynomials."
Serdica Journal of Computing 9(2) (2015), 123-138.
"""
# compute the modified subresultant prs
lst = modified_subresultants_pg(p,q,x) ## any other method would do
# defensive
if lst == [] or len(lst) == 2:
return lst
# the coefficients in lst are modified subresultants and, hence, are
# greater than those of the corresponding subresultants by the factor
# LC(lst[0])**( deg(lst[0]) - deg(lst[1])); see Theorem 2 in reference.
lcf = LC(lst[0])**( degree(lst[0], x) - degree(lst[1], x) )
# Initialize the subresultant prs list
subr_seq = [lst[0], lst[1]]
# compute the degree sequences m_i and j_i of Theorem 2 in reference.
deg_seq = [degree(Poly(poly, x), x) for poly in lst]
deg = deg_seq[0]
deg_seq_s = deg_seq[1:-1]
m_seq = [m-1 for m in deg_seq_s]
j_seq = [deg - m for m in m_seq]
# compute the AMV factors of Theorem 2 in reference.
fact = [(-1)**( j*(j-1)/S(2) ) for j in j_seq]
# shortened list without the first two polys
lst_s = lst[2:]
# poly lst_s[k] is multiplied times fact[k], divided by lcf
# and appended to the subresultant prs list
m = len(fact)
for k in range(m):
if sign(fact[k]) == -1:
subr_seq.append(-lst_s[k] / lcf)
else:
subr_seq.append(lst_s[k] / lcf)
return subr_seq
def subresultants_amv_q(p, q, x):
"""
p, q are polynomials in Z[x] or Q[x]. It is assumed
that degree(p, x) >= degree(q, x).
Computes the subresultant prs of p and q in Q[x];
the coefficients of the polynomials in the sequence are
subresultants. That is, they are determinants of appropriately
selected submatrices of sylvester1, Sylvester's matrix of 1840.
To compute the coefficients, no determinant evaluation takes place.
Instead, polynomial divisions in Q[x] are performed, using the
function rem(p, q, x); the coefficients of the remainders
computed this way become subresultants with the help of the
Akritas-Malaschonok-Vigklas Theorem of 2015.
If the subresultant prs is complete, then it coincides with the
Euclidean sequence of the polynomials p, q.
References
==========
1. Akritas, A. G., G.I. Malaschonok and P.S. Vigklas: ``A Basic Result
on the Theory of Subresultants.'' Serdica Journal of Computing 10 (2016), No.1, 31-48.
2. Akritas, A. G., G.I. Malaschonok and P.S. Vigklas: ``Subresultant Polynomial
remainder Sequences Obtained by Polynomial Divisions in Q[x] or in Z[x].''
Serdica Journal of Computing 10 (2016), No.3-4, 197-217.
"""
# make sure neither p nor q is 0
if p == 0 or q == 0:
return [p, q]
# make sure proper degrees
d0 = degree(p, x)
d1 = degree(q, x)
if d0 == 0 and d1 == 0:
return [p, q]
if d1 > d0:
d0, d1 = d1, d0
p, q = q, p
if d0 > 0 and d1 == 0:
return [p, q]
# initialize
i, s = 0, 0 # counters for remainders & odd elements
p_odd_index_sum = 0 # contains the sum of p_1, p_3, etc
subres_l = [p, q] # subresultant prs output list
a0, a1 = p, q # the input polys
sigma1 = LC(a1, x) # leading coeff of a1
p0 = d0 - d1 # degree difference
if p0 % 2 == 1:
s += 1
phi = floor( (s + 1) / 2 )
mul_fac = 1
d2 = d1
# main loop
while d2 > 0:
i += 1
a2 = rem(a0, a1, domain= QQ) # new remainder
if i == 1:
sigma2 = LC(a2, x)
else:
sigma3 = LC(a2, x)
sigma1, sigma2 = sigma2, sigma3
d2 = degree(a2, x)
p1 = d1 - d2
psi = i + phi + p_odd_index_sum
# new mul_fac
mul_fac = sigma1**(p0 + 1) * mul_fac
## compute the sign of the first fraction in formula (9) of the paper
# numerator
num = (-1)**psi
# denominator
den = sign(mul_fac)
# the sign of the determinant depends on sign( num / den ) != 0
if sign(num / den) > 0:
subres_l.append( simplify(expand(a2* Abs(mul_fac))))
else:
subres_l.append(- simplify(expand(a2* Abs(mul_fac))))
## bring into mul_fac the missing power of sigma if there was a degree gap
if p1 - 1 > 0:
mul_fac = mul_fac * sigma1**(p1 - 1)
# update AMV variables
a0, a1, d0, d1 = a1, a2, d1, d2
p0 = p1
if p0 % 2 ==1:
s += 1
phi = floor( (s + 1) / 2 )
if i%2 == 1:
p_odd_index_sum += p0 # p_i has odd index
# gcd is of degree > 0 ?
m = len(subres_l)
if subres_l[m - 1] == nan or subres_l[m - 1] == 0:
subres_l.pop(m - 1)
return subres_l
def compute_sign(base, expo):
'''
base != 0 and expo >= 0 are integers;
returns the sign of base**expo without
evaluating the power itself!
'''
sb = sign(base)
if sb == 1:
return 1
pe = expo % 2
if pe == 0:
return -sb
else:
return sb
def rem_z(p, q, x):
'''
Intended mainly for p, q polynomials in Z[x] so that,
on dividing p by q, the remainder will also be in Z[x]. (However,
it also works fine for polynomials in Q[x].) It is assumed
that degree(p, x) >= degree(q, x).
It premultiplies p by the _absolute_ value of the leading coefficient
of q, raised to the power deg(p) - deg(q) + 1 and then performs
polynomial division in Q[x], using the function rem(p, q, x).
By contrast the function prem(p, q, x) does _not_ use the absolute
value of the leading coefficient of q.
This results not only in ``messing up the signs'' of the Euclidean and
Sturmian prs's as mentioned in the second reference,
but also in violation of the main results of the first and third
references --- Theorem 4 and Theorem 1 respectively. Theorems 4 and 1
establish a one-to-one correspondence between the Euclidean and the
Sturmian prs of p, q, on one hand, and the subresultant prs of p, q,
on the other.
References
==========
1. Akritas, A. G., G.I. Malaschonok and P.S. Vigklas: ``On the Remainders
Obtained in Finding the Greatest Common Divisor of Two Polynomials.''
Serdica Journal of Computing, 9(2) (2015), 123-138.
2. https://planetMath.org/sturmstheorem
3. Akritas, A. G., G.I. Malaschonok and P.S. Vigklas: ``A Basic Result on
the Theory of Subresultants.'' Serdica Journal of Computing 10 (2016), No.1, 31-48.
'''
if (p.as_poly().is_univariate and q.as_poly().is_univariate and
p.as_poly().gens == q.as_poly().gens):
delta = (degree(p, x) - degree(q, x) + 1)
return rem(Abs(LC(q, x))**delta * p, q, x)
else:
return prem(p, q, x)
def quo_z(p, q, x):
"""
Intended mainly for p, q polynomials in Z[x] so that,
on dividing p by q, the quotient will also be in Z[x]. (However,
it also works fine for polynomials in Q[x].) It is assumed
that degree(p, x) >= degree(q, x).
It premultiplies p by the _absolute_ value of the leading coefficient
of q, raised to the power deg(p) - deg(q) + 1 and then performs
polynomial division in Q[x], using the function quo(p, q, x).
By contrast the function pquo(p, q, x) does _not_ use the absolute
value of the leading coefficient of q.
See also function rem_z(p, q, x) for additional comments and references.
"""
if (p.as_poly().is_univariate and q.as_poly().is_univariate and
p.as_poly().gens == q.as_poly().gens):
delta = (degree(p, x) - degree(q, x) + 1)
return quo(Abs(LC(q, x))**delta * p, q, x)
else:
return pquo(p, q, x)
def subresultants_amv(f, g, x):
"""
p, q are polynomials in Z[x] or Q[x]. It is assumed
that degree(f, x) >= degree(g, x).
Computes the subresultant prs of p and q in Z[x] or Q[x];
the coefficients of the polynomials in the sequence are
subresultants. That is, they are determinants of appropriately
selected submatrices of sylvester1, Sylvester's matrix of 1840.
To compute the coefficients, no determinant evaluation takes place.
Instead, polynomial divisions in Z[x] or Q[x] are performed, using
the function rem_z(p, q, x); the coefficients of the remainders
computed this way become subresultants with the help of the
Akritas-Malaschonok-Vigklas Theorem of 2015 and the Collins-Brown-
Traub formula for coefficient reduction.
If the subresultant prs is complete, then it coincides with the
Euclidean sequence of the polynomials p, q.
References
==========
1. Akritas, A. G., G.I. Malaschonok and P.S. Vigklas: ``A Basic Result
on the Theory of Subresultants.'' Serdica Journal of Computing 10 (2016), No.1, 31-48.
2. Akritas, A. G., G.I. Malaschonok and P.S. Vigklas: ``Subresultant Polynomial
remainder Sequences Obtained by Polynomial Divisions in Q[x] or in Z[x].''
Serdica Journal of Computing 10 (2016), No.3-4, 197-217.
"""
# make sure neither f nor g is 0
if f == 0 or g == 0:
return [f, g]
# make sure proper degrees
d0 = degree(f, x)
d1 = degree(g, x)
if d0 == 0 and d1 == 0:
return [f, g]
if d1 > d0:
d0, d1 = d1, d0
f, g = g, f
if d0 > 0 and d1 == 0:
return [f, g]
# initialize
a0 = f
a1 = g
subres_l = [a0, a1]
deg_dif_p1, c = degree(a0, x) - degree(a1, x) + 1, -1
# initialize AMV variables
sigma1 = LC(a1, x) # leading coeff of a1
i, s = 0, 0 # counters for remainders & odd elements
p_odd_index_sum = 0 # contains the sum of p_1, p_3, etc
p0 = deg_dif_p1 - 1
if p0 % 2 == 1:
s += 1
phi = floor( (s + 1) / 2 )
# compute the first polynomial of the prs
i += 1
a2 = rem_z(a0, a1, x) / Abs( (-1)**deg_dif_p1 ) # first remainder
sigma2 = LC(a2, x) # leading coeff of a2
d2 = degree(a2, x) # actual degree of a2
p1 = d1 - d2 # degree difference
# sgn_den is the factor, the denominator 1st fraction of (9),
# by which a2 is multiplied to get integer coefficients
sgn_den = compute_sign( sigma1, p0 + 1 )
## compute sign of the 1st fraction in formula (9) of the paper
# numerator
psi = i + phi + p_odd_index_sum
num = (-1)**psi
# denominator
den = sgn_den
# the sign of the determinant depends on sign(num / den) != 0
if sign(num / den) > 0:
subres_l.append( a2 )
else:
subres_l.append( -a2 )
# update AMV variable
if p1 % 2 == 1:
s += 1
# bring in the missing power of sigma if there was gap
if p1 - 1 > 0:
sgn_den = sgn_den * compute_sign( sigma1, p1 - 1 )
# main loop
while d2 >= 1:
phi = floor( (s + 1) / 2 )
if i%2 == 1:
p_odd_index_sum += p1 # p_i has odd index
a0, a1, d0, d1 = a1, a2, d1, d2 # update polys and degrees
p0 = p1 # update degree difference
i += 1
sigma0 = -LC(a0)
c = (sigma0**(deg_dif_p1 - 1)) / (c**(deg_dif_p1 - 2))
deg_dif_p1 = degree(a0, x) - d2 + 1
a2 = rem_z(a0, a1, x) / Abs( (c**(deg_dif_p1 - 1)) * sigma0 )
sigma3 = LC(a2, x) # leading coeff of a2
d2 = degree(a2, x) # actual degree of a2
p1 = d1 - d2 # degree difference
psi = i + phi + p_odd_index_sum
# update variables
sigma1, sigma2 = sigma2, sigma3
# new sgn_den
sgn_den = compute_sign( sigma1, p0 + 1 ) * sgn_den
# compute the sign of the first fraction in formula (9) of the paper
# numerator
num = (-1)**psi
# denominator
den = sgn_den
# the sign of the determinant depends on sign( num / den ) != 0
if sign(num / den) > 0:
subres_l.append( a2 )
else:
subres_l.append( -a2 )
# update AMV variable
if p1 % 2 ==1:
s += 1
# bring in the missing power of sigma if there was gap
if p1 - 1 > 0:
sgn_den = sgn_den * compute_sign( sigma1, p1 - 1 )
# gcd is of degree > 0 ?
m = len(subres_l)
if subres_l[m - 1] == nan or subres_l[m - 1] == 0:
subres_l.pop(m - 1)
return subres_l
def modified_subresultants_amv(p, q, x):
"""
p, q are polynomials in Z[x] or Q[x]. It is assumed
that degree(p, x) >= degree(q, x).
Computes the modified subresultant prs of p and q in Z[x] or Q[x],
from the subresultant prs of p and q.
The coefficients of the polynomials in the two sequences differ only
in sign and the factor LC(p)**( deg(p)- deg(q)) as stated in
Theorem 2 of the reference.
The coefficients of the polynomials in the output sequence are
modified subresultants. That is, they are determinants of appropriately
selected submatrices of sylvester2, Sylvester's matrix of 1853.
If the modified subresultant prs is complete, and LC( p ) > 0, it coincides
with the (generalized) Sturm's sequence of the polynomials p, q.
References
==========
1. Akritas, A. G., G.I. Malaschonok and P.S. Vigklas: "On the Remainders
Obtained in Finding the Greatest Common Divisor of Two Polynomials."
Serdica Journal of Computing, Serdica Journal of Computing, 9(2) (2015), 123-138.
"""
# compute the subresultant prs
lst = subresultants_amv(p,q,x) ## any other method would do
# defensive
if lst == [] or len(lst) == 2:
return lst
# the coefficients in lst are subresultants and, hence, smaller than those
# of the corresponding modified subresultants by the factor
# LC(lst[0])**( deg(lst[0]) - deg(lst[1])); see Theorem 2.
lcf = LC(lst[0])**( degree(lst[0], x) - degree(lst[1], x) )
# Initialize the modified subresultant prs list
subr_seq = [lst[0], lst[1]]
# compute the degree sequences m_i and j_i of Theorem 2
deg_seq = [degree(Poly(poly, x), x) for poly in lst]
deg = deg_seq[0]
deg_seq_s = deg_seq[1:-1]
m_seq = [m-1 for m in deg_seq_s]
j_seq = [deg - m for m in m_seq]
# compute the AMV factors of Theorem 2
fact = [(-1)**( j*(j-1)/S(2) ) for j in j_seq]
# shortened list without the first two polys
lst_s = lst[2:]
# poly lst_s[k] is multiplied times fact[k] and times lcf
# and appended to the subresultant prs list
m = len(fact)
for k in range(m):
if sign(fact[k]) == -1:
subr_seq.append( simplify(-lst_s[k] * lcf) )
else:
subr_seq.append( simplify(lst_s[k] * lcf) )
return subr_seq
def correct_sign(deg_f, deg_g, s1, rdel, cdel):
"""
Used in various subresultant prs algorithms.
Evaluates the determinant, (a.k.a. subresultant) of a properly selected
submatrix of s1, Sylvester's matrix of 1840, to get the correct sign
and value of the leading coefficient of a given polynomial remainder.
deg_f, deg_g are the degrees of the original polynomials p, q for which the
matrix s1 = sylvester(p, q, x, 1) was constructed.
rdel denotes the expected degree of the remainder; it is the number of
rows to be deleted from each group of rows in s1 as described in the
reference below.
cdel denotes the expected degree minus the actual degree of the remainder;
it is the number of columns to be deleted --- starting with the last column
forming the square matrix --- from the matrix resulting after the row deletions.
References
==========
Akritas, A. G., G.I. Malaschonok and P.S. Vigklas: ``Sturm Sequences
and Modified Subresultant Polynomial Remainder Sequences.''
Serdica Journal of Computing, Vol. 8, No 1, 29-46, 2014.
"""
M = s1[:, :] # copy of matrix s1
# eliminate rdel rows from the first deg_g rows
for i in range(M.rows - deg_f - 1, M.rows - deg_f - rdel - 1, -1):
M.row_del(i)
# eliminate rdel rows from the last deg_f rows
for i in range(M.rows - 1, M.rows - rdel - 1, -1):
M.row_del(i)
# eliminate cdel columns
for i in range(cdel):
M.col_del(M.rows - 1)
# define submatrix
Md = M[:, 0: M.rows]
return Md.det()
def subresultants_rem(p, q, x):
"""
p, q are polynomials in Z[x] or Q[x]. It is assumed
that degree(p, x) >= degree(q, x).
Computes the subresultant prs of p and q in Z[x] or Q[x];
the coefficients of the polynomials in the sequence are
subresultants. That is, they are determinants of appropriately
selected submatrices of sylvester1, Sylvester's matrix of 1840.
To compute the coefficients polynomial divisions in Q[x] are
performed, using the function rem(p, q, x). The coefficients
of the remainders computed this way become subresultants by evaluating
one subresultant per remainder --- that of the leading coefficient.
This way we obtain the correct sign and value of the leading coefficient
of the remainder and we easily ``force'' the rest of the coefficients
to become subresultants.
If the subresultant prs is complete, then it coincides with the
Euclidean sequence of the polynomials p, q.
References
==========
1. Akritas, A. G.:``Three New Methods for Computing Subresultant
Polynomial Remainder Sequences (PRS's).'' Serdica Journal of Computing 9(1) (2015), 1-26.
"""
# make sure neither p nor q is 0
if p == 0 or q == 0:
return [p, q]
# make sure proper degrees
f, g = p, q
n = deg_f = degree(f, x)
m = deg_g = degree(g, x)
if n == 0 and m == 0:
return [f, g]
if n < m:
n, m, deg_f, deg_g, f, g = m, n, deg_g, deg_f, g, f
if n > 0 and m == 0:
return [f, g]
# initialize
s1 = sylvester(f, g, x, 1)
sr_list = [f, g] # subresultant list
# main loop
while deg_g > 0:
r = rem(p, q, x)
d = degree(r, x)
if d < 0:
return sr_list
# make coefficients subresultants evaluating ONE determinant
exp_deg = deg_g - 1 # expected degree
sign_value = correct_sign(n, m, s1, exp_deg, exp_deg - d)
r = simplify((r / LC(r, x)) * sign_value)
# append poly with subresultant coeffs
sr_list.append(r)
# update degrees and polys
deg_f, deg_g = deg_g, d
p, q = q, r
# gcd is of degree > 0 ?
m = len(sr_list)
if sr_list[m - 1] == nan or sr_list[m - 1] == 0:
sr_list.pop(m - 1)
return sr_list
def pivot(M, i, j):
'''
M is a matrix, and M[i, j] specifies the pivot element.
All elements below M[i, j], in the j-th column, will
be zeroed, if they are not already 0, according to
Dodgson-Bareiss' integer preserving transformations.
References
==========
1. Akritas, A. G.: ``A new method for computing polynomial greatest
common divisors and polynomial remainder sequences.''
Numerische MatheMatik 52, 119-127, 1988.
2. Akritas, A. G., G.I. Malaschonok and P.S. Vigklas: ``On a Theorem
by Van Vleck Regarding Sturm Sequences.''
Serdica Journal of Computing, 7, No 4, 101-134, 2013.
'''
ma = M[:, :] # copy of matrix M
rs = ma.rows # No. of rows
cs = ma.cols # No. of cols
for r in range(i+1, rs):
if ma[r, j] != 0:
for c in range(j + 1, cs):
ma[r, c] = ma[i, j] * ma[r, c] - ma[i, c] * ma[r, j]
ma[r, j] = 0
return ma
def rotate_r(L, k):
'''
Rotates right by k. L is a row of a matrix or a list.
'''
ll = list(L)
if ll == []:
return []
for i in range(k):
el = ll.pop(len(ll) - 1)
ll.insert(0, el)
return ll if isinstance(L, list) else Matrix([ll])
def rotate_l(L, k):
'''
Rotates left by k. L is a row of a matrix or a list.
'''
ll = list(L)
if ll == []:
return []
for i in range(k):
el = ll.pop(0)
ll.insert(len(ll) - 1, el)
return ll if isinstance(L, list) else Matrix([ll])
def row2poly(row, deg, x):
'''
Converts the row of a matrix to a poly of degree deg and variable x.
Some entries at the beginning and/or at the end of the row may be zero.
'''
k = 0
poly = []
leng = len(row)
# find the beginning of the poly ; i.e. the first
# non-zero element of the row
while row[k] == 0:
k = k + 1
# append the next deg + 1 elements to poly
for j in range( deg + 1):
if k + j <= leng:
poly.append(row[k + j])
return Poly(poly, x)
def create_ma(deg_f, deg_g, row1, row2, col_num):
'''
Creates a ``small'' matrix M to be triangularized.
deg_f, deg_g are the degrees of the divident and of the
divisor polynomials respectively, deg_g > deg_f.
The coefficients of the divident poly are the elements
in row2 and those of the divisor poly are the elements
in row1.
col_num defines the number of columns of the matrix M.
'''
if deg_g - deg_f >= 1:
print('Reverse degrees')
return
m = zeros(deg_f - deg_g + 2, col_num)
for i in range(deg_f - deg_g + 1):
m[i, :] = rotate_r(row1, i)
m[deg_f - deg_g + 1, :] = row2
return m
def find_degree(M, deg_f):
'''
Finds the degree of the poly corresponding (after triangularization)
to the _last_ row of the ``small'' matrix M, created by create_ma().
deg_f is the degree of the divident poly.
If _last_ row is all 0's returns None.
'''
j = deg_f
for i in range(0, M.cols):
if M[M.rows - 1, i] == 0:
j = j - 1
else:
return j if j >= 0 else 0
def final_touches(s2, r, deg_g):
"""
s2 is sylvester2, r is the row pointer in s2,
deg_g is the degree of the poly last inserted in s2.
After a gcd of degree > 0 has been found with Van Vleck's
method, and was inserted into s2, if its last term is not
in the last column of s2, then it is inserted as many
times as needed, rotated right by one each time, until
the condition is met.
"""
R = s2.row(r-1)
# find the first non zero term
for i in range(s2.cols):
if R[0,i] == 0:
continue
else:
break
# missing rows until last term is in last column
mr = s2.cols - (i + deg_g + 1)
# insert them by replacing the existing entries in the row
i = 0
while mr != 0 and r + i < s2.rows :
s2[r + i, : ] = rotate_r(R, i + 1)
i += 1
mr -= 1
return s2
def subresultants_vv(p, q, x, method = 0):
"""
p, q are polynomials in Z[x] (intended) or Q[x]. It is assumed
that degree(p, x) >= degree(q, x).
Computes the subresultant prs of p, q by triangularizing,
in Z[x] or in Q[x], all the smaller matrices encountered in the
process of triangularizing sylvester2, Sylvester's matrix of 1853;
see references 1 and 2 for Van Vleck's method. With each remainder,
sylvester2 gets updated and is prepared to be printed if requested.
If sylvester2 has small dimensions and you want to see the final,
triangularized matrix use this version with method=1; otherwise,
use either this version with method=0 (default) or the faster version,
subresultants_vv_2(p, q, x), where sylvester2 is used implicitly.
Sylvester's matrix sylvester1 is also used to compute one
subresultant per remainder; namely, that of the leading
coefficient, in order to obtain the correct sign and to
force the remainder coefficients to become subresultants.
If the subresultant prs is complete, then it coincides with the
Euclidean sequence of the polynomials p, q.
If the final, triangularized matrix s2 is printed, then:
(a) if deg(p) - deg(q) > 1 or deg( gcd(p, q) ) > 0, several
of the last rows in s2 will remain unprocessed;
(b) if deg(p) - deg(q) == 0, p will not appear in the final matrix.
References
==========
1. Akritas, A. G.: ``A new method for computing polynomial greatest
common divisors and polynomial remainder sequences.''
Numerische MatheMatik 52, 119-127, 1988.
2. Akritas, A. G., G.I. Malaschonok and P.S. Vigklas: ``On a Theorem
by Van Vleck Regarding Sturm Sequences.''
Serdica Journal of Computing, 7, No 4, 101-134, 2013.
3. Akritas, A. G.:``Three New Methods for Computing Subresultant
Polynomial Remainder Sequences (PRS's).'' Serdica Journal of Computing 9(1) (2015), 1-26.
"""
# make sure neither p nor q is 0
if p == 0 or q == 0:
return [p, q]
# make sure proper degrees
f, g = p, q
n = deg_f = degree(f, x)
m = deg_g = degree(g, x)
if n == 0 and m == 0:
return [f, g]
if n < m:
n, m, deg_f, deg_g, f, g = m, n, deg_g, deg_f, g, f
if n > 0 and m == 0:
return [f, g]
# initialize
s1 = sylvester(f, g, x, 1)
s2 = sylvester(f, g, x, 2)
sr_list = [f, g]
col_num = 2 * n # columns in s2
# make two rows (row0, row1) of poly coefficients
row0 = Poly(f, x, domain = QQ).all_coeffs()
leng0 = len(row0)
for i in range(col_num - leng0):
row0.append(0)
row0 = Matrix([row0])
row1 = Poly(g,x, domain = QQ).all_coeffs()
leng1 = len(row1)
for i in range(col_num - leng1):
row1.append(0)
row1 = Matrix([row1])
# row pointer for deg_f - deg_g == 1; may be reset below
r = 2
# modify first rows of s2 matrix depending on poly degrees
if deg_f - deg_g > 1:
r = 1
# replacing the existing entries in the rows of s2,
# insert row0 (deg_f - deg_g - 1) times, rotated each time
for i in range(deg_f - deg_g - 1):
s2[r + i, : ] = rotate_r(row0, i + 1)
r = r + deg_f - deg_g - 1
# insert row1 (deg_f - deg_g) times, rotated each time
for i in range(deg_f - deg_g):
s2[r + i, : ] = rotate_r(row1, r + i)
r = r + deg_f - deg_g
if deg_f - deg_g == 0:
r = 0
# main loop
while deg_g > 0:
# create a small matrix M, and triangularize it;
M = create_ma(deg_f, deg_g, row1, row0, col_num)
# will need only the first and last rows of M
for i in range(deg_f - deg_g + 1):
M1 = pivot(M, i, i)
M = M1[:, :]
# treat last row of M as poly; find its degree
d = find_degree(M, deg_f)
if d is None:
break
exp_deg = deg_g - 1
# evaluate one determinant & make coefficients subresultants
sign_value = correct_sign(n, m, s1, exp_deg, exp_deg - d)
poly = row2poly(M[M.rows - 1, :], d, x)
temp2 = LC(poly, x)
poly = simplify((poly / temp2) * sign_value)
# update s2 by inserting first row of M as needed
row0 = M[0, :]
for i in range(deg_g - d):
s2[r + i, :] = rotate_r(row0, r + i)
r = r + deg_g - d
# update s2 by inserting last row of M as needed
row1 = rotate_l(M[M.rows - 1, :], deg_f - d)
row1 = (row1 / temp2) * sign_value
for i in range(deg_g - d):
s2[r + i, :] = rotate_r(row1, r + i)
r = r + deg_g - d
# update degrees
deg_f, deg_g = deg_g, d
# append poly with subresultant coeffs
sr_list.append(poly)
# final touches to print the s2 matrix
if method != 0 and s2.rows > 2:
s2 = final_touches(s2, r, deg_g)
pprint(s2)
elif method != 0 and s2.rows == 2:
s2[1, :] = rotate_r(s2.row(1), 1)
pprint(s2)
return sr_list
def subresultants_vv_2(p, q, x):
"""
p, q are polynomials in Z[x] (intended) or Q[x]. It is assumed
that degree(p, x) >= degree(q, x).
Computes the subresultant prs of p, q by triangularizing,
in Z[x] or in Q[x], all the smaller matrices encountered in the
process of triangularizing sylvester2, Sylvester's matrix of 1853;
see references 1 and 2 for Van Vleck's method.
If the sylvester2 matrix has big dimensions use this version,
where sylvester2 is used implicitly. If you want to see the final,
triangularized matrix sylvester2, then use the first version,
subresultants_vv(p, q, x, 1).
sylvester1, Sylvester's matrix of 1840, is also used to compute
one subresultant per remainder; namely, that of the leading
coefficient, in order to obtain the correct sign and to
``force'' the remainder coefficients to become subresultants.
If the subresultant prs is complete, then it coincides with the
Euclidean sequence of the polynomials p, q.
References
==========
1. Akritas, A. G.: ``A new method for computing polynomial greatest
common divisors and polynomial remainder sequences.''
Numerische MatheMatik 52, 119-127, 1988.
2. Akritas, A. G., G.I. Malaschonok and P.S. Vigklas: ``On a Theorem
by Van Vleck Regarding Sturm Sequences.''
Serdica Journal of Computing, 7, No 4, 101-134, 2013.
3. Akritas, A. G.:``Three New Methods for Computing Subresultant
Polynomial Remainder Sequences (PRS's).'' Serdica Journal of Computing 9(1) (2015), 1-26.
"""
# make sure neither p nor q is 0
if p == 0 or q == 0:
return [p, q]
# make sure proper degrees
f, g = p, q
n = deg_f = degree(f, x)
m = deg_g = degree(g, x)
if n == 0 and m == 0:
return [f, g]
if n < m:
n, m, deg_f, deg_g, f, g = m, n, deg_g, deg_f, g, f
if n > 0 and m == 0:
return [f, g]
# initialize
s1 = sylvester(f, g, x, 1)
sr_list = [f, g] # subresultant list
col_num = 2 * n # columns in sylvester2
# make two rows (row0, row1) of poly coefficients
row0 = Poly(f, x, domain = QQ).all_coeffs()
leng0 = len(row0)
for i in range(col_num - leng0):
row0.append(0)
row0 = Matrix([row0])
row1 = Poly(g,x, domain = QQ).all_coeffs()
leng1 = len(row1)
for i in range(col_num - leng1):
row1.append(0)
row1 = Matrix([row1])
# main loop
while deg_g > 0:
# create a small matrix M, and triangularize it
M = create_ma(deg_f, deg_g, row1, row0, col_num)
for i in range(deg_f - deg_g + 1):
M1 = pivot(M, i, i)
M = M1[:, :]
# treat last row of M as poly; find its degree
d = find_degree(M, deg_f)
if d is None:
return sr_list
exp_deg = deg_g - 1
# evaluate one determinant & make coefficients subresultants
sign_value = correct_sign(n, m, s1, exp_deg, exp_deg - d)
poly = row2poly(M[M.rows - 1, :], d, x)
poly = simplify((poly / LC(poly, x)) * sign_value)
# append poly with subresultant coeffs
sr_list.append(poly)
# update degrees and rows
deg_f, deg_g = deg_g, d
row0 = row1
row1 = Poly(poly, x, domain = QQ).all_coeffs()
leng1 = len(row1)
for i in range(col_num - leng1):
row1.append(0)
row1 = Matrix([row1])
return sr_list