/* 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. */