module quartic; % Procedures for solving cubic, quadratic and quartic
% eqns.
% Author: Anthony C. Hearn.
% Modifications by: Stanley L. Kameny.
fluid '(!*sub2 !*rounded !*trigform dmode!*);
!*trigform := t; % Default value.
switch trigform;
symbolic procedure multfq(u,v);
% Multiplies standard form U by standard quotient V.
begin scalar x;
x := gcdf(u,denr v);
return multf(quotf(u,x),numr v) ./ quotf(denr v,x)
end;
symbolic procedure quotsqf(u,v);
% Forms quotient of standard quotient U and standard form V.
begin scalar x;
x := gcdf(numr u,v);
return quotf(numr u,x) ./ multf(quotf(v,x),denr u)
end;
symbolic procedure cubertq u;
% Rationalizing the value in this and the following function leads
% usually to neater results.
% rationalizesq
simpexpt list(mk!*sq subs2!* u,'(quotient 1 3));
% simprad(u,3);
symbolic procedure sqrtq u;
% rationalizesq
simpexpt list(mk!*sq subs2!* u,'(quotient 1 2));
% simprad(u,2);
% symbolic procedure subs2!* u; <<!*sub2 := t; subs2 u>>;
symbolic procedure solvequadratic(a2,a1,a0);
% A2, a1 and a0 are standard quotients.
% Solves a2*x**2+a1*x+a0=0 for x.
% Returns a list of standard quotient solutions.
% Modified to use root_val to compute numeric roots. SLK.
if !*rounded and numcoef a0 and numcoef a1 and numcoef a2
then for each z in cdr root_val list mkpolyexp2(a2,a1,a0)
collect simp!* z else
begin scalar d;
d := sqrtq subtrsq(quotsqf(exptsq(a1,2),4),multsq(a2,a0));
a1 := quotsqf(negsq a1,2);
return list(subs2!* quotsq(addsq(a1,d),a2),
subs2!* quotsq(subtrsq(a1,d),a2))
end;
symbolic procedure numcoef a; denr a = 1 and domainp numr a;
symbolic procedure mkpolyexp2(a2,a1,a0);
% The use of 'x is arbitrary here, since it is not used by root_val.
<<a0 := numr a0;
if numr a1 then a0 := (('x . 1) . numr a1) . a0;
mk!*sq(((('x . 2) . numr a2) . a0) . 1)>>;
symbolic procedure solvecubic(a3,a2,a1,a0);
% A3, a2, a1 and a0 are standard quotients.
% Solves a3*x**3+a2*x**2+a1*x+a0=0 for x.
% Returns a list of standard quotient solutions.
% See Abramowitz and Stegun, Sect. 3.8.2, for details.
begin scalar q,r,sm,sp,s1,s2,x;
a2 := quotsq(a2,a3);
a1 := quotsq(a1,a3);
a0 := quotsq(a0,a3);
q := subtrsq(quotsqf(a1,3),quotsqf(exptsq(a2,2),9));
r := subtrsq(quotsqf(subtrsq(multsq(a1,a2),multfq(3,a0)),6),
quotsqf(exptsq(a2,3),27));
if null numr q or not !*trigform or not all_real(a0,a1,a2)
then go to cbr;
% this section uses trig functions, but only when a0,a1,a2 are real.
x := sqrtq negsq addsq(exptsq(q,3),exptsq(r,2));
if one_real simp list('times,'i,mk!*sq x) and not pos_num q
then x := negsq x;
s1 := quotsqf(atan2q(x,r),3);
s2 := negsq sqrtq negsq q;
sp := negsq multfq(2,multsq(s2,cossq s1));
sm := multsq(simp '(sqrt 3),multsq(s2,sinsq s1));
% sp := -2*sqrt(-q)*cos(atan2(sqrt( - q**3 - r**2),r)/3)$
% sm := - sqrt(-q)*sqrt(3)*sin(atan2(sqrt( - q**3 - r**2),r)/3)$
go to com;
cbr: x := sqrtq addsq(exptsq(q,3),exptsq(r,2));
s1 := cubertq addsq(r,x);
s2 := if numr s1 then negsq quotsq(q,s1)
else cubertq subtrsq(r,x);
% This optimization only works if s1 is non zero.
sp := addsq(s1,s2);
sm := quotsqf(multsq(simp '(times i (sqrt 3)),subtrsq(s1,s2)),2);
com: x := subtrsq(sp,quotsqf(a2,3));
sp := negsq addsq(quotsqf(sp,2),quotsqf(a2,3));
return list(subs2!* x,subs2!* addsq(sp,sm),
subs2!* subtrsq(sp,sm))
end;
symbolic procedure pos_num a;
begin scalar dmode,!*msg,!*numval;
dmode := dmode!*;
!*numval := t;
on rounded,complex;
a := resimp a;
a := real_1 a and (numr simp list('sign,mk!*sq a)=1);
off rounded,complex;
if dmode then onoff(get(dmode,'dname),t);
return a
end;
symbolic procedure sinsq a;
simpiden list('sin,mk!*sq subs2!* a);
symbolic procedure cossq a;
simpiden list('cos,mk!*sq subs2!* a);
symbolic procedure all_real(a,b,c);
begin scalar dmode,!*ezgcd,!*msg,!*numval;
% If ezgcd is on, modular arithmetic with rounded numbers can be
% attempted.
dmode := dmode!*;
!*numval := t;
on complex,rounded;
% We should probably put an errorset here, so mode is correctly
% reset with an error.
a := real_1(a := resimp a) and real_1(b := resimp b)
and real_1(c := resimp c);
off rounded,complex;
if dmode then onoff(get(dmode,'dname),t);
return a
end;
symbolic procedure real_1 x;
numberp denr x and domainp numr x and null numr impartsq x;
symbolic procedure one_real a;
begin scalar dmode,!*msg,!*numval;
dmode := dmode!*;
!*numval := t;
on complex,rounded;
a := real_1 resimp a;
off rounded,complex;
if dmode then onoff(get(dmode,'dname),t);
return a end;
symbolic procedure atan2q(b,a);
% Used by solvecubic to set up trig form expressions for atan2 in
% terms of atan and, where necessary, a bias of +/- pi; or to call
% atan2 directly if numerical solution is called for.
((begin scalar !*msg,x,y,dmode,q,fg,s1,s2,s3,s4,s5;
y := b := simp!*(b := mk!*sq subs2!* b);
x := a := simp!*(a := mk!*sq subs2!* a);
if domainp numr y and domainp numr x
and numberp denr y and numberp denr x then go to aret;
dmode := dmode!*;
on complex,rounded;
y := resimp b; x := resimp a;
if not(domainp numr y and domainp numr x
and numberp denr y and numberp denr x) then go to ret;
q := sqrtq addsq(exptsq(x,2),exptsq(y,2));
x := quotsq(x,q); y := quotsq(y,q);
if null numr x then
<<s1 := t;
if numr simp list('sqn,list('repart,mk!*sq y))=-1
then s2 := t;
go to ret>>;
s3 := t;
x := numr simp list('sign,list('repart,mk!*sq x));
if x=-1 then
<<y := numr simp list('sign,list('repart,mk!*sq y));
if y=-1 then s4 := t else s5 := t>>;
ret: off rounded,complex;
if dmode then onoff(get(dmode,'dname),t);
if s1 then
fg := quotsqf(simp 'pi,2);
if s2 then fg := negsq fg;
if s3 then fg := simpiden list('atan,mk!*sq quotsq(b,a));
if s4 then fg := subtrsq(fg,simp 'pi);
if s5 then fg := addsq(fg,simp 'pi);
aret: return if fg then fg else
simpiden list('atan2,mk!*sq b,mk!*sq a) end)
where !*numval=t);
symbolic procedure solvequartic(a4,a3,a2,a1,a0);
% Solve the quartic equation a4*x**4+a3*x**3+a2*x**2+a1*x+a0 = 0,
% where the ai are standard quotients, using technique described in
% Section 3.8.3 of Abramowitz and Stegun;
begin scalar x,y,yy,cx,z,s,l,zz1,zz2,dmode,neg,!*msg,!*numval;
% Convert equation to monomial form.
dmode := dmode!*;
a3 := quotsq(a3,a4);
a2 := quotsq(a2,a4);
a1 := quotsq(a1,a4);
a0 := quotsq(a0,a4);
% Build and solve the resultant cubic equation. We select the
% real root if there is only one; or if there are three, we choose
% one that yields real coefficients for the quadratics. If no
% roots are known to be real, we use an arbitrary one.
yy := subtrsq(exptsq(a3,2),multfq(4,a2));
x := solvecubic(!*f2q 1,
negsq a2,
subs2!* subtrsq(multsq(a1,a3),multfq(4,a0)),
subs2!* negsq addsq(exptsq(a1,2),
multsq(a0,yy)));
cx := car x;
% Now check for real roots of the cubic.
for each rr in x do if one_real rr then s := append(s,list rr);
x := if (l := length s)=1 then car s else cx;
% Now solve the two equivalent quadratic equations.
a3 := quotsqf(a3,2); yy := quotsqf(yy,4);
% select real coefficient for quadratic if possible.
y := addsq(yy,x);
if l<2 then go to zz;
loop: if not pos_num negsq y then go to zz else if l=1 then
<<x := cx; y := addsq(yy,x); go to zz>>;
l := l-1; s := cdr s; x := car s;
y := addsq(yy,x); go to loop;
zz: y := sqrtq y;
x := quotsqf(x,2);
z := sqrtq subtrsq(exptsq(x,2),a0);
% the following test is needed, according to some editions of
% Abramowitz and Stegun, to select the correct signs
% (for the terms z) in the quadratics to produce correct roots.
% Unfortunately, this test may fail for coefficients which are not
% numeric because of the inability to recognize zero.
!*numval := t;
on rounded,complex;
if null numr
(zz1 :=
resimp subtrsq(a1,addsq(multsq(subtrsq(a3,y),addsq(x,z)),
multsq(addsq(a3,y),subtrsq(x,z))))) then go to rst;
if null numr
(zz2 :=
resimp subtrsq(a1,addsq(multsq(subtrsq(a3,y),subtrsq(x,z)),
multsq(addsq(a3,y),addsq(x,z)))))
then <<neg := t; go to rst>>;
if domainp numr zz1 and domainp numr zz2
and numberp denr zz1 and numberp denr zz2 and
numr simp list('sign,list('difference,list('norm,mk!*sq zz1),
list('norm,mk!*sq zz2)))=1 then neg := t;
rst: off rounded,complex;
if dmode then onoff(get(dmode,'dname),t);
if neg then z := negsq z;
return append(solvequadratic(!*f2q 1,subtrsq(a3,y),subtrsq(x,z)),
solvequadratic(!*f2q 1,addsq(a3,y),addsq(x,z)))
end;
endmodule;
end;