Following is some code I wrote to provide a factorization of a quartic
into quadratics, derived from A&S. Comments?
Dan Stanger
define_variable(qfactor_signs,[[-1,-1],[-1,1],[1,-1],[1,1]],list);
qfactor(quartic,var,[sfun]):=block([a0,a1,a2,a3,u,alpha,beta,gamma,delta,f1,f2],
if(length(sfun) > 0) then sfun:sfun[1] else sfun:ratsimp,
l:coeff(quartic,var^4),
quartic:expand(quartic/l),
a0:first(partition(quartic,var)),
a3:coeff(quartic,var,3),
a2:coeff(quartic,var,2),
a1:coeff(quartic,var,1),
u:qfactor_solve(var^3-a2*var^2+(a1*a3-4*a0)*var-(a1^2+a0*a3^2-4*a0*a2),var),
block([flag:false,res,u1],
for u1 in u do (
print("solving with root",u1),
alpha:a3/2,beta:sqrt(a3^2/4+u1-a2),
gamma:u1/2,delta:sqrt((u1/2)^2-a0),
for i in qfactor_signs do (
for j in qfactor_signs do (
res:qfactor_prod(var,i,j,alpha,beta,gamma,delta,a0,a1,a2,a3,sfun),
flag:first(res),
if flag = true then (f1:res[2],f2:res[3],return(0))
),
if flag = true then return(0)
),
if flag = true then return(0)
)
),
l*f1*f2
);
qfactor_prod(var,sign1,sign2,alpha,beta,gamma,delta,a0,a1,a2,a3,sfun):=block(
[p1,q1,p2,q2],
p1:alpha+sign1[1]*beta,
q1:gamma+sign1[2]*delta,
p2:alpha+sign2[1]*beta,
q2:gamma+sign2[2]*delta,
if apply(sfun,[(p1+p2)=a3]) and
apply(sfun,[(p1*p2+q1+q2)=a2]) and
apply(sfun,[(p1*q2+p2*q1)=a1]) and
apply(sfun,[(q1*q2)=a0])
then [true,var^2+p1*var+q1,var^2 +p2*var+q2] else [false]
);
qfactor_solve(cubic,var):=block([r],
r:[],
for i in map(rhs,solve(cubic,var)) do
if freeof(%i,i) then push(i,r),
r);