reduce-historical/r37/packages/solve/quartic.red

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;