module specfac; % Splitting of low degree polynomials.
% Author: Anthony C. Hearn.
% Copyright (c) 1991 RAND. All rights reserved.
fluid '(!*keepsqrts !*sub2 !*surds knowndiscrimsign kord!* zlist);
% switch surds;
exports cubicf,quadraticf,quarticf;
symbolic procedure coeffs pol;
% Extract coefficients of polynomial wrt its main variable and leading
% degree. Result is a list of coefficients.
begin integer degree,deg1; scalar cofs,mv;
mv := mvar pol;
degree := ldeg pol;
while not domainp pol and mvar pol eq mv do
<<deg1 := ldeg pol;
for i:= 1:(degree-deg1-1) do cofs := nil . cofs;
cofs := lc pol . cofs;
pol := red pol;
degree := deg1>>;
for i:=1:degree-1 do cofs := nil . cofs;
return reversip(pol . cofs)
end;
symbolic procedure shift!-pol pol;
% Shifts main variable, mv, of square free nth degree polynomial pol so
% that coefficient of mv**(n-1) is zero.
% Does not assume pol is univariate.
begin scalar lc1,ld,mv,pol1,redp,shift,x;
mv := mvar pol;
ld := ldeg pol;
redp := red pol;
if domainp redp or not(mvar redp eq mv) or ldeg redp<(ld-1)
then return list(pol,1,nil ./ 1);
lc1 := lc pol;
x := lc redp;
shift := quotsq(!*f2q x,!*f2q multd(ld,lc1));
pol1 := subf1(pol,list(mv . mk!*sq addsq(!*k2q mv,negsq shift)));
return list(numr pol1,denr pol1,shift)
end;
symbolic procedure quadraticf!*(pol,var);
if domainp pol then errach "invalid quadratic to factr"
else if mvar pol = var then quadraticf pol
else begin scalar kord,w;
kord := kord!*;
kord!* := list var;
w := coeffs !*q2f resimp(pol ./ 1);
kord!* := kord;
w := quadraticf1(car w,cadr w,caddr w);
if w eq 'failed then return list(1,pol);
var := !*k2f var;
return list(if car w neq 1 then mkrn(1,car w) else 1,
addf(multf(var,cadr w),caddr w),
addf(multf(var,cadddr w),car cddddr w))
end;
symbolic procedure quadraticf pol;
% Finds factors of square free quadratic polynomial pol (if they
% exist). Does not assume pol is univariate.
(if x eq 'failed then list(1,pol)
else if not domainp car x then list(1,pol)
% Answer would be rational.
else list(if car x neq 1 then mkrn(1,car x) else 1,
y .* cadr x .+ caddr x,y .* cadddr x .+ car cddddr x)
where y = (mvar pol .** 1))
where x = quadraticf1(car w,cadr w,caddr w) where w = coeffs pol;
symbolic procedure quadraticf1(a,b,c);
begin scalar a1,denom,discrim,w;
if null b and minusf c and not minusf a
then <<a := rootxf(a,2);
c := rootxf(negf c,2);
return if a eq 'failed or c eq 'failed then 'failed
else list(1,a,c,a,negf c)>>;
discrim := powsubsf addf(exptf(b,2),multd(-4,multf(a,c)));
% A null discriminator can arise from a polynomial such as
% 16x^2+(32i-8)*x-8i-15;
if null discrim then nil
else <<if knowndiscrimsign
then <<if knowndiscrimsign eq 'negative
then return 'failed>>
% else if not clogflag and minusf discrim
% then return 'failed;
else if minusf discrim then return 'failed;
discrim:=rootxf(discrim,2);
if discrim='failed then return discrim>>;
denom := multd(4,a);
a := a1 := multd(2,a);
w := addf(b,discrim);
c := addf(b,negf discrim);
b := w;
if (w := gcdf(a,b)) neq 1
then <<a1 := quotf(a,w); b := quotf(b,w);
denom := quotf(denom,w)>>;
if (w := gcdf(a,denom)) neq 1 and (w := gcdf(c,denom))
then <<a := quotf(a,w);
c := quotf(c,w);
denom := quotf(denom,w)>>;
return list(denom,a1,b,a,c)
end;
symbolic procedure rootxf(u,n);
% Return either polynomial nth root of u or "failed".
begin scalar x,y,z,w;
if domainp u
then return if minusf u then 'failed
else if atom u and (y := irootn(u,n))**n=u then y
else if not atom u and (x := get(car u,'rootfn))
then apply2(x,u,n)
else if !*surds and not(u member zlist)
then nrootn!*(u,n)
else 'failed;
x := comfac u;
u := quotf(u,comfac!-to!-poly x);
z := 1;
if car x then if cdr(y := divide(cdar x,n)) = 0
then z := multpf(caar x .** car y,z)
else if !*surds
then <<z := multf(mkrootf(caar x,n,cdr y),z);
if car y neq 0
then z := multpf(caar x .** car y,z)>>
else return 'failed;
x := cdr x;
if domainp x
then if minusf x then return 'failed
else if fixp x and (y := irootn(x,n))**n=x
then z := multd(y,z)
else if !*surds and fixp x
then z := multf(nrootn!*(x,n),z)
else if not atom x and (w := get(car x,'rootfn))
then apply2(w,x,n)
else return 'failed
else if (y := rootxf(x,n)) eq 'failed then return y
else z := multf(y,z);
if u=1 then return z;
x := sqfrf u;
c: if null x then return z
else if cdr(y := divide(cdar x,n)) = 0
then <<z := multf(exptf(caar x,car y),z); x := cdr x>>
else if !*surds
then <<z := multf(mkrootf(prepf caar x,n,cdr y),
multf(exptf(caar x,car y),z));
x := cdr x>>
else return 'failed;
go to c
end;
symbolic procedure mkrootf(u,m,n);
if m neq 2 or null !*keepsqrts
then !*p2f mksp(list('expt,u,list('quotient,1,m)),n)
else if n neq 1 then errach 'mkrootf
else !*q2f simpsqrt list u;
symbolic procedure nrootn!*(u,n);
% Returns a standard form representation of the nth root of u.
begin scalar x;
if null u then return nil;
u := nrootn(u,n);
x := cdr u; % surd part.
u := car u; % rational part.
if x=1 then return x;
x := mkrootf(prepf x,n,1);
return powsubsf multf(u,x)
end;
symbolic procedure cubicf pol;
% Split the cubic pol if a change of origin puts it in the form
% (x-a)**3-b=0.
begin scalar a,a0,a1,b,neg,p;
p := shift!-pol pol;
a := coeffs car p;
if cadr a then return list(1,pol)
% Cadr a non nil probably means there are some surds in the
% coefficients that don't reduce to 0.
else if caddr a then return list(1,pol);
% Factorization not possible by this method.
a0 := cadddr a;
a := car a;
if minusf a0 then <<neg := t; a0 := negf a0>>;
if (a := rootxf(a,3)) eq 'failed
or (a0 := rootxf(a0,3)) eq 'failed
then return list(1,pol);
if neg then a0 := negf a0;
a := !*f2q a;
a0 := !*f2q a0;
p := addsq(!*k2q mvar pol,caddr p);
% Now numr (a*(mv+shift)+a0) is a factor of pol.
a1 := numr addsq(multsq(a,p),a0);
% quotf(pol,a) is quadratic factor. However, the surd division may
% not work properly, so we calculate factor directly.
b := multsq(a0,a0);
b := addsq(b,multsq(negsq multsq(a,a0),p));
b := numr addsq(b,multsq(multsq(a,a),exptsq(p,2)));
return aconc!*(quadraticf b,a1)
end;
symbolic procedure powsubsf u;
% We believe that the result of this operation must be a polynomial.
% If subs2q returns a rational, it must be because there are
% unsimplified surds. Hopefully rationalizesq can fix those.
begin scalar !*sub2;
u := subs2q !*f2q u;
if denr u neq 1
then <<u := rationalizesq u;
if denr u neq 1 then errach list('powsubsf,u)>>;
return numr u
end;
symbolic procedure quarticf pol;
% Splits quartics that can be written in the form
% (x-a)**4+b*(x-a)**2+c.
% Note that any call of rootxf can lead to a result "failed."
begin scalar !*sub2,a,a2,a0,b,dsc,p,p1,p2,q,shift,var;
var := mvar pol;
p := shift!-pol pol;
a := coeffs car p;
shift := caddr p;
if cadr a % pol not correctly shifted, possibly due to sqrt.
% e.g., 729para^4*be^4 - 81para^3*sqrt(27*be^2*para^2 - 8cte1^3)*
% sqrt(3)*be^3 - 216para^2*be^2*cte1^3 + 12para*sqrt(27be^2*para^2
% - 8*cte1^3)*sqrt(3) *be*cte1^3 + 8*cte1^6.
or cadddr a then return list(1,pol);
% Factorization not possible by this method.
a2 := cddr a;
a0 := caddr a2;
a2 := car a2;
a := car a;
q := quadraticf1(a,a2,a0);
if not(q eq 'failed)
then <<a2 := car q; q := cdr q;
a := exptsq(addsq(!*k2q mvar pol,shift),2);
b := numr subs2q quotsq(addsq(multsq(!*f2q car q,a),
!*f2q cadr q),
!*f2q cadr p);
a := numr subs2q quotsq(addsq(multsq(!*f2q caddr q,a),
!*f2q cadddr q),
!*f2q cadr p);
a := quadraticf!*(a,var);
b := quadraticf!*(b,var);
return multf(a2,multf(car a,car b))
. nconc!*(cdr a,cdr b)>>
else if null !*surds or denr shift neq 1
then return list(1,pol);
% Factorization not possible by this method.
shift := numr shift;
if knowndiscrimsign eq 'negative then go to complex;
dsc := powsubsf addf(exptf(a2,2),multd(-4,multf(a,a0)));
p2 := minusf a0;
if not p2 and minusf dsc then go to complex;
p1 := not a2 or minusf a2;
if not p1 then if p2 then p1 := t else p2 := t;
p1 := if p1 then 'positive else 'negative;
p2 := if p2 then 'negative else 'positive;
a := rootxf(a,2);
if a eq 'failed then return list(1,pol);
dsc := rootxf(dsc,2);
if dsc eq 'failed then return list(1,pol);
p := invsq !*f2q addf(a,a);
q := multsq(!*f2q addf(a2,negf dsc),p);
p := multsq(!*f2q addf(a2,dsc),p);
b := multf(a,exptf(addf(!*k2f mvar pol,shift),2));
a := powsubsf addf(b,q);
b := powsubsf addf(b,p);
knowndiscrimsign := p1;
a := quadraticf!*(a,var);
knowndiscrimsign := p2;
b := quadraticf!*(b,var);
knowndiscrimsign := nil;
return multf(car a,car b) . nconc!*(cdr a,cdr b);
% Complex case.
complex:
a := rootxf(a,2);
if a eq 'failed then return list(1,pol);
a0 := rootxf(a0,2);
if a0 eq 'failed then return list(1,pol);
a2 := powsubsf addf(multf(2,multf(a,a0)),negf a2);
a2 := rootxf(a2,2);
if a2 eq 'failed then return list(1,pol);
% Now a*(x+shift)**2 (+/-) b*(x+shift) + c is a factor.
p := addf(!*k2f mvar pol,shift);
q := addf(multf(a,exptf(p,2)),a0);
p := multf(a2,p);
a := powsubsf addf(q,p);
b := powsubsf addf(q,negf p);
knowndiscrimsign := 'negative;
a := quadraticf!*(a,var);
b := quadraticf!*(b,var);
knowndiscrimsign := nil;
return multf(car a,car b) . nconc!*(cdr a,cdr b)
end;
endmodule;
end;