/* The aim is to study the subresultant PRS algorithm due to Brown,
which is implemented in maxima. References for resultants and subresultants is
https://en.wikipedia.org/wiki/Polynomial_greatest_common_divisor
and for a detailed study of the algorithm:
"On Euclid's Algorithm and the Theory of Subresultants", Brown & Traub,
JACM, Vol. 18, No. 4, October 1971, pp. 505-514.
"The subresultant PRS algorithm", Brown,
ACM Trans. math.soft., Vol 4, No 3, September 1978, pp. 237-249 */
/* Resultants and subresultants */
F:a_4+a_3*x+a_2*x^2+a_1*x^3+a_0*x^4;
G:b_3+b_2*x+b_1*x^2+b_0*x^3;
coeffs(P,x):=block (local(l),l:[],
for i from hipow(P,x) step -1 thru 0 do l:cons(coeff(P,x,i),l),
reverse(l))$
/* Beware all matrices are transposed with respect to those in Wikipedia
for printing convenience */
result(P,Q,x):=block(local(mat,len1,len2,ll1,ll2),
len1:hipow(P,x)+1,
len2:hipow(Q,x)+1, /* assume len1 >= len2 */
mat:zeromatrix(len1+len2-2,len1+len2-2),
ll1:coeffs(P,x),
ll2:coeffs(Q,x),
for i from 1 thru len2-1 do (
for j from i thru i+len1-1 do (
mat[i,j]:ll1[j-i+1])),
for i from len2 thru len2+len1-2 do (
for j from i-len2+1 thru i do (
mat[i,j]:ll2[j-i+len2])),
mat)$ /* This is the usual matrix resultant */
/* To do one step of reduction, remove 2 lines and last column */
subres1(mat,d1,d2):= submatrix(d2,d1+d2,mat,length(transpose(mat)))$
subres(P,Q,j,x):=block (local(mat,d1,d2,V),
mat:result(P,Q,x),
d1:hipow(P,x),
d2:hipow(Q,x),
V:zeromatrix(d1+d2-2*j,d1+d2-j), /* See Wikipedia */
for i from 1 thru d1+d2-2*j-1 do V[i,i]:1,
for i from 0 thru j do V[d1+d2-2*j,d1+d2-j-i]:x^i,
if (j >= min(d1,d2)) then print(j," too big")
else for i from 1 thru j do (
mat:subres1(mat,d1,d2),
d1:d1-1, d2:d2-1),
mat.transpose(V))$ /* Matrix product */
for j from 0 thru hipow(G,x)-1 do print(subres(F,G,j,x));
/* It is the determinants of these objects that we need */
/* This is the example in Brown */
r[0]: x^8+x^6-3*x^4-3*x^3+8*x^2+2*x-5$
r[1]: 3*x^6+5*x^4-4*x^2-9*x+21$
/* newdet is much faster than determinant, in this computation. */
for j from 0 thru 5 do print (j," ---> ",expand(newdet(subres(r[0],r[1],j,x))));
/* Result:
0 ---> 260708
1 ---> 9326 x - 12300
2
2 ---> 169 x + 325 x - 637
2
3 ---> 65 x + 125 x - 245
4 2
4 ---> 25 x - 5 x + 15
4 2
5 ---> 15 x - 3 x + 9
*/
/* If the "fundamental theorem in Brown and Traub, or Brown was correct one
should get two vanishing polynomials here, while polynomials of degree 2 and 4
are in fact repeated */
/* Brown's polynomial remainder sequence */
/* Builds a sequence of polynomials of decreasing degrees G_1,G_2,G_3 such that
G_(i+1) is proportional to the remainder of G_(i-1) by G_i. the coefficient is such that
all the G_i have integral coeffs if G_O and G_1 have, and each G_i is up to sign a subresultant
of G_0 and G_1 thus the last non vanishing one is the GCD of G_0 and G_1.
Moreover the growth of the coeffs in G_i is limited thanks to a similar fact for
subresultants.
In the algorithm g_i is the leading coeff of G_i and the numerical factor
g[i]^(d[i]+1)/b[i] achieves the desired result. Finally b[i] is computed recursively
using an intermediate h[i] described in Brown's paper.
Comment in maxima source code (from Rupert Swarbrick)
;; Write G[1] = P, G[2] = Q and write g[i] for the leading coefficient of
;; G[i]. On the k'th iteration of the loop, we are computing G[k+2], which is
;; given by the formula
;;
;; G[k+2] = (-1)^(delta[k]+1) prem(G[k], G[k+1]) / (g[k] h[k]^delta[k])
;;
;; except we set the denominator to 1 when k=1. Here, h[2] = g[2]^delta[1] and,
;; for i >= 3, satisfies the recurrence:
;;
;; h[i] = g[i]^delta[i-1] h[i-1]^(1 - delta[i-1])
;;
;; Here, delta[i] = deg(G[i]) - deg(G[i+1]), which is non-negative.
;;
*/
r[0]: x^8+x^6-3*x^4-3*x^3+8*x^2+2*x-5$
r[1]: 3*x^6+5*x^4-4*x^2-9*x+21$
/* an example without holes in the sequence of degrees */
r[2]: x^6+x^5-3*x^4-3*x^3+8*x^2+2*x-5$
r[3]: 3*x^4+5*x^3-4*x^2-9*x+21$
lc(P,x):=coeff(P,x,hipow(P,x))$
prs_result(P,Q,x):= block(local(delta,g,G,b,h),
delta[k]:=hipow(G[k],x)-hipow(G[k+1],x),
G[1]:P,
G[2]:Q,
print(1, "---> ", G[1]),
for k from 1 unless ( G[k+1] = 0) do (
g[k]: lc(G[k],x),
g[k+1]: lc(G[k+1],x),
if k=1 then (
b[1]: 1, h[1]: 1) else (
h[k]: g[k]^delta[k-1]*h[k-1]^(1 - delta[k-1]),
b[k]: g[k]*h[k]^delta[k]),
G[k+2]: (-1)^(delta[k]+1)*
expand(remainder(g[k+1]^(delta[k]+1)*G[k],G[k+1],x)/b[k]),
print(k+1, "---> ", G[k+1])))$
prs_result(r[0],r[1],x)$
prs_result(r[2],r[3],x)$
/* The result for Brown's example
(%i15) 6 4 2
1 ---> 3 x + 5 x - 4 x - 9 x + 21
4 2
2 ---> 15 x - 3 x + 9
2
3 ---> 65 x + 125 x - 245
4 ---> 9326 x - 12300
5 ---> 260708
The second example
(%i16) 4 3 2
1 ---> 3 x + 5 x - 4 x - 9 x + 21
3 2
2 ---> (- x ) + 47 x - 135 x + 30
2
3 ---> 239 x - 727 x + 163
4 ---> 1152 - 34516 x
5 ---> 2899120
*/
/* Note one gets the same results as by direct computation of determinants,
a remark first noted by Habicht (1948), then Brown & Traub and Collins. The
explanation comes from studying how one step of polynomial division translates
on subresultants, namely when F=QG+R one can basically replace F by R in
result(F,G) up to trivial signs and factors, because GQ leads to repeated
lines in the determinant, so one gets to result(G,R) and iterate downto a
single polynomial in the remainder sequence. Here is an example:
subres(r[0],r[1],2,x);
content(expand(newdet(%)));
yields
[13, 13*x^2 + 25*x - 49]
But the third polynomial in the prs chain and remainder of the division
of r[0] by r[1] is:
r[2]: 15*x^4-3*x^2+9;
subres(r[1],r[2],2,x);
content(expand(newdet(%)));
yields the thourth polynomial of the prs chain
[85293, 13*x^2 + 25*x - 49]
r[3]: 13*x^2+25*x-49;
subres(r[2],r[3],2,x); yields
2 too big
But one can do:
(%i32) subres(r[2],r[3],1,x);
[ 15 0 - 3 9 ]
[ ]
[ 13 25 - 49 0 ]
(%o32) [ ]
[ 0 13 25 - 49 x ]
[ ]
[ 0 0 13 25 x - 49 ]
content(expand(newdet(%))) yields
[150, 4663*x - 6150]
and twice %[2] is just the fifth polynomial above.
So a chain of subresultants can also end up on a polynomial
of the prs chain without ending on zero. The main point in Brown
is to explicit the factors and signs appearing in this computation. */