Deflating by a quadratic of complex roots...

Suppose you’ve found a root z=a+bi, complex. Then, provided the coefficients of the polynomial p(x) are real, you’ve found a second root: the complex conjugate, (OverBar[z])=a-bi. Thus, we might deflate by the quadratic
(x-z)(x-(OverBar[z]))=(x^2)-2Re[z]x+z(OverBar[z]).
Let’s call the coefficients of the quadratic [1,α=-2Re[z],β=z(OverBar[z])].

Now we want to factor p(x), using the quadratic:

p(x)=r(x) ((x^2)-2Re[z]x+z(OverBar[z])) + f(x,z,(OverBar[z]))

It used to be that we just had a “p(t)” there, where we now have an f(x,z,(OverBar[z])). 

In general that function is 

f(x,z,(OverBar[z])) = (p(z)-p((OverBar[z]))/(z-(OverBar[z]))x - (p(z)(OverBar[z])-p((OverBar[z]))z/(z-(OverBar[z]))

which is obtained using Lagrange polynomials, which we will study later on: 

f(x,z,(OverBar[z])) = p(z)(x-(OverBar[z]))/(z-(OverBar[z]))+p((OverBar[z]))(x-z)/((OverBar[z])-z)

But because z and (OverBar[z]) are roots, f(x,z,(OverBar[z]))=0 (or so we say!).

Complex pair deflation:

Having found a complex root, we remove both it and its conjugate by synthetically dividing out a quadratic.

NOTE: the coefficients go in in order, from leading coefficient on down. So the polynomial

p(x)=-4x^3+2x^2+3x-6

is represented as {-4,2,3,-6}.

deflateComplexPair[c_,z_,verbose_]:=
Module[{n=Length[c],newcoefs,
next,nexter,nextest,previous,
alpha=-2Re[z],beta=Re[z*Conjugate[z]]},
next=c[[1]];
nexter=c[[2]]-alpha*next;
newcoefs={next};
AppendTo[newcoefs,nexter];
For[i=3,i<=n-2,i++,
If[verbose,Print[c[[i]]," ",next]];
nextest=c[[i]]-alpha*nexter-beta*next;
If[verbose,Print[i," ",c[[i]]," ",nexter," ",next]];
AppendTo[newcoefs,nextest];
previous=next;
next=nexter;
nexter=nextest;
If[verbose,Print[{previous,next,nexter}]];
];
If[verbose,Print[{c[[n-1]],next,nexter}]];
If[verbose,Print["Should be zero: ",c[[n-1]]-alpha*nexter-beta*next]];
If[verbose,Print[{c[[n]],next,nexter}]];
If[verbose,Print["Should be zero: ",c[[n]]-beta*nexter]];
newcoefs
]
deflateComplexPair[{1,0,0,0,1},-(((-1))^(1/4)),True]
deflateComplexPair[{1,0,0,0,1},((1-I)/(Sqrt[2])),True]
deflateComplexPair[{1,0,0,0,-1},I,True]
deflateComplexPair[{1,1,1,1},I,True]

deflateComplexPair[{1,0,0,0,-1},I,True]

(* Tea time polynomial for section 2.6, p. 88: *)
remainder=deflateComplexPair[{-4,2,3,-6},.84493+.73994I,True]
-remainder[[2]]/remainder[[1]]

quadraticRoots[a_,b_,c_]:=
Module[{tmp=b+If[b==0,1,Sign[b]]Sqrt[b^2-4a*c]},
-{2*c/tmp,tmp/(2*a)}
]
quadraticRoots[1,0,-1]

quadraticRootsList[list_]:=quadraticRoots[list[[1]],list[[2]],list[[3]]]
quadraticRootsList[{1,-(Sqrt[2]),1}]
quadraticRootsList[{1,(Sqrt[2]),1}]

coefs={1,0,0,0,1};
p[x_,list_]:=Total[Table[list[[i]]x^(i-1),{i,Length[list],1,-1}]]
Plot[p[x,coefs],{x,-2,2}]
Solve[p[x,coefs]==0,x]

coefs={1,1,1,1};
p[x_,list_]:=Total[Table[list[[i]]x^(i-1),{i,Length[list],1,-1}]]
Plot[p[x,coefs],{x,-2,2}]
Solve[p[x,coefs]==0,x]