Artifact 7e3eee161db61273fda42975d2aa3909d5fbd21795cc9571cc1de1d60cc9eb90:
- Executable file
r38/packages/solve/quartic.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 8457) [annotate] [blame] [check-ins using] [more...]
module quartic; % Procedures for solving cubic, quadratic and quartic % eqns. % Author: Anthony C. Hearn. % Modifications by: Stanley L. Kameny, Eberhard Schruefer. 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!* (if eqcar(z,'equal) then caddr z else errach {"Quadratic confusion",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 solutions in % terms of surds and Bronstein for solutions in terms of % trig functions. 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. s2 := sqrtq simp {'abs,prepsq q}; if pos_num r then s2 := negsq s2; if pos_num q then <<s1 := quotsqf(trigsq(quotsq(negsq r,exptsq(s2,3)),'asinh),3); sp := trigsq(s1,'sinh); sm := multsq(simp '(times (sqrt 3) i),trigsq(s1,'cosh))>> else if pos_num addsq(exptsq(q,3),exptsq(r,2)) then <<s1 := quotsqf(trigsq(quotsq(negsq r,exptsq(s2,3)),'acosh),3); sp := trigsq(s1,'cosh); sm := multsq(simp '(times (sqrt 3) i),trigsq(s1,'sinh))>> else <<s1 := quotsqf(trigsq(quotsq(negsq r,exptsq(s2,3)),'acos),3); sp := trigsq(s1,'cos); sm := multsq(simp '(sqrt 3),trigsq(s1,'sin))>>; return {subs2!* subtrsq(multsq(s2,multsq(-2 ./ 1,sp)),quotsqf(a2,3)), subs2!* subtrsq(multsq(s2,addsq(sp,sm)),quotsqf(a2,3)), subs2!* subtrsq(multsq(s2,subtrsq(sp,sm)),quotsqf(a2,3))}; 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 trigsq(a,fn); simpiden list(fn,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 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;