Sylvester Matrix

We compute the Sylvester Matrix as a LIST OF THE LISTS OF ITS ROWS :
Example: (({{1, 2}, {3, 4}, {5, 6}})) is represented as { {1,2} , {3,4} , {5,6} }

Syl[A_,B_,x_]:=(dA=Exponent[A,x];(*compute the Sylvester Matrix of the polynomials A and B in x*)
dB=Exponent[B,x];(* dB=deg(B) *)
ca=Reverse[CoefficientList[A,x]];(*List of coefficients of A taken from increasing degree, them reversed =taken in decreaing degree *)
cb=Reverse[CoefficientList[B,x]];
ca=Join[ca,Table[0,{dB-1}]]; (* add dB-1 zeros to the list ca  *)
cb=Join[cb,Table[0,{dA-1}]]; 
Ta=List[ca];Tb=List[cb];(* Construct the first row with 'List[ca]' and the (dB+1)-th row with 'List[cb]' *)
Do[(ca=RotateRight[ca]; (* add successively the row 'RotateRioght[ca]' to the previously inserted lines *)
Ta=Join[Ta,List[ca]]),{dB-1}];  (* at this point we have Ta=(({{(Subscript[a, dA]), (Subscript[a, dA-1]), ⋯, (Subscript[a, 0]), 0, 0, 0}, {0, (Subscript[a, dA]), \[Placeholder], \[Placeholder], (Subscript[a, 0]), 0, 0}, {0, 0, ⋱, \[Placeholder], \[Placeholder], ⋱, 0}, {0, 0, 0, (Subscript[a, dA]), ⋯, ⋯, (Subscript[a, dA])}})) *)
Do[(cb=RotateRight[cb];Tb=Join[Tb,List[cb]]),{dA-1}];
Return[Join[Ta,Tb]];)

Example of computation of Row Echelon Form:

A=x^4 - 3 x^2+ 2 x

B=x^3 -1

Sab=Syl[A,B,x]

MatrixForm[Sab]

MatrixForm[RowReduce[%]]

PolynomialGCD[A,B]

Intersection of 2 plane curves

EXMAPLE 1: Non-degenerate case (no points at the infinity)

A=(X^2+Y^2)^3 - 4 X^2 Y^2

B=X^2 (1+Y)-(1-Y)^3

ContourPlot[{A==0,B==0},{X,-2,2},{Y,-1.5,1.5}]

Resultant[A,B,Y]

Resultant[A,B,X]

(* Resultant has degree 14, but we see only 4 solutions on the picture. Why ?*)

(* Why the resultant in X has only only even power of X ?*)

(* Why the resultant in Y has is a square ?*)

EXMAPLE 2: Degenrate case (points at the ∞)

A=X^2 Y + X + 1

B=X Y -1

ContourPlot[{A==0,B==0},{X,-2,2},{Y,-4,1}]

Resultant[A,B,Y]

Factor[%]

(*Is the root of the resultant x=0 a solution of the system {A=0, B=0 ?*)

Change of coordinates to go back to a non-degenerate case

(* We perform the change: X=X'+1  *)

AA=A/.X->X+Y

BB=B/.X->X+Y

ContourPlot[{AA==0,BB==0},{X,-10,10},{Y,-10,10}]

Resultant[AA,BB,Y]

(* The resultant has degree 1 only, corresponding to the
coordinte on the x-axis if the unique common affine solution
of AA and BB  *)

(* The linear change of cordinate do not change the number of affine solutions, so there is also only one common solution for
for A and B *)

aa=X Y^2 - 2 X^2 Y + X^3

aa/.X->X+Y

Simplify[%]

Another example: no complex solutions (we can see all of them on the picture )

aa=1+(x/2)-((x^2)/2)+7 y-8 x y+2 (x^2) y-3 (y^3)+4 x (y^3)-(x^2) (y^3)

bb=Expand[aa+2 b]

bb=-11+((45 x)/2)-((25 (x^2))/2)+2 (x^3)+7 y-8 x y+2 (x^2) y-3 (y^3)+4 x (y^3)-(x^2) (y^3)

Sab=Syl[aa,bb,y];MatrixForm[Sab]

r=Det[Sab]

Factor[r]

ContourPlot[{aa==0,bb==0},{x,0,6},{y,-5,5}]

(* ABOVE x=2 : 3 intersection points above x=2 ; this is the multiplicity of the factor (x-2) -> no complex solutions *)

(* ABOVE x=1 : 1 intersection point and another at the infinity *)

(* ABOVE x=3 : at the infinity only *)

Change of coordinates to avoid points at the infinity

aaa=aa/.x->x+y

bbb=bb/.x->x+y

ContourPlot[{aaa==0,bbb==0},{x,0,5},{y,-5,5}]

Resultant[aaa,bbb,y]

Factor[%]

CHECK THE FOLLOWING ON THE PICTURE :

(*  1 intersection point above x=3 *)
(* 2 intersectin points above x=2 *)
(* 2 intersection points above x=1 *)

Growth of coefficients in remainders sequence of the Euclidean algorithm

Standard Euclidean algorithm

Eucl[A_, B_, x_] := (A1 = A; A2 = B; Atmp = 0;
  While[A2 =!= 0,
   Atmp = PolynomialRemainder[A1, A2, x];
   Print[Atmp];
   A1 = A2;
   A2 = Simplify[Atmp];
   ];
  Return[A1];)

Euclidean algorithm for computing the resultant

ResEucl[A_,B_,x_]:=(A1=A;A2=B;A3=0;R1=1;R2=0;
d2=Exponent[A2,x];d1=Exponent[A1,x];d3=0;
While[d2>0,
A3=PolynomialRemainder[A1,A2,x];Print["r=",A3];
A3=Simplify[A3];
d3=Exponent[A3,x];
R2=(-1)^(d1 *d2) * R1*Coefficient[A2,x,d2]^(d1-d3);Print["R2=",R2];
(*Exchange   until here, i=2  *)
A1=A2;
A2=A3;
d1=d2;
d2=d3;
R1=R2;	
];
If[A2=!=0,
Return[R1*Coefficient[A2,x,Exponent[A2,x]]^d1],
Return[0];
];)

t1=(x+1)(x+2)(x+3)(x+15)(x-14)(x+10)^2

t2=x*(x-19)*14

ResEucl[t1,t2,x]

1) For Polynomials A,B in Q[X]

p=Eucl[(y+1)*(y+2)+3,y+1,y]

A=824 x^5 - 65 x^4 - 814 x^3 -741 x^2 -979x-764

B=216 x^4 + 663 x^3 +880 x^2 + 916 x+ 617

Eucl[A,B,x]

2) For Polynomials A,B in Q[X,Y] (  computed in Q(Y)[X]  )

(* A lot of repetitions in the successive remainders. Useless computations are done...much better method exit: 1) the subresultant algorithm (but tricky, so not discussed in this class)
2) modular method *)

p=(1-y+y^2)+(3+y)x + (3+2 y^2) x^2 - (3+y) x^3+(1+129 y)x^4

q=y-(y+1)x +(1+y+y^2)x^2

r=Eucl[p,q,x]

pp=(1-2*y+3*y^2)+ y* x + (2-3 y ) x^2 - (3-y+y^2)x^3 + (y-1)x^4

qq=y+(y+1) x -(1+y^2 +y^3)x^2

Eucl[pp,qq,x]

ResEucl[pp,qq,x]

Resultant[pp,qq,x]