Artifact d1df22435478c9b509cd53d6ae6d0451029a8a9c705378cecb74a179fd09b07b:
- Executable file
r37/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: 9778) [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. 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;