Artifact 28e6a0b1d7e5ef09a10cf3bfa61bbf3e4ebd4bcee0bb9c7218e43e67a7020cfa:
- Executable file
r37/packages/poly/conj.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: 6425) [annotate] [blame] [check-ins using] [more...]
module conj; % Rationalize denoms of standard quotients by conjugate % computation. % Author: Anthony C. Hearn. % Modifications by: Eberhard Schruefer. % Copyright (c) 1992 RAND. All rights reserved. fluid '(!*algint !*rationalize !*structure dmode!* kord!* powlis!*); put('rationalize,'simpfg,'((t (rmsubs)) (nil (rmsubs)))); symbolic smacro procedure subtrf(u,v); % Returns u - v for standard forms u and v. addf(u,negf v); symbolic procedure rationalizesq u; % Rationalize the standard quotient u. begin scalar !*structure,!*sub2,v,x; % Modified by R. Liska. % We need structure off to form rationalized denominator properly % in subs2f1. % ACH had hoped that the cost of having GCD on here was small, % since the consequences can be large (e.g., df(log((sqrt(a^2+x^2) % +2*sqrt(sqrt(a^2+x^2)*a+a*x)+a+x)/(sqrt(a^2+x^2) - a + x)),x)). % However, limit((sqrt(x^(2/5) +1) - x^(1/3)-1)/x^(1/3),x,0) takes % too long. if x := get(dmode!*,'rationalizefn) then u := apply1(x,u); % We need the following in case we are in sparse_bareiss. powlis!* := '(i 2 (nil . t) -1 nil) . powlis!*; v := subs2q u; powlis!* := cdr powlis!*; % We need the subs2 to get rid of surd powers. % We also need to check if u has changed from the example % df((1/x)**(2/3),x). return if domainp denr v then v else if (x := rationalizef denr v) neq 1 then <<v := multf(numr v,x) ./ multf(denr v,x); % We need the gcdchk so that df((1/x)**(2/3),x) % is not in a loop. However, algint needs all % factors for some reason. if null !*algint and null !*rationalize then v := gcdchk v; % There used to be an exptchk here, but that led % to loops (e.g., in df(int(1/(3*x**4+7),x),x)). v := subs2q v; if not domainp numr quotsq(u,v) then rationalizesq v else v>> else u end; symbolic procedure lowertowerp(u,v); % True if v is potentially an algebraic component of a member of v. if null u then nil else if atom car u or cdar u = v then lowertowerp(cdr u,v) else if caar u eq 'expt and eqcar(caddar u,'quotient) and cadr caddar u = cadr cadr v % numerator of quotient. and fixp caddr caddar u and fixp caddr cadr v and cdr divide(caddr caddar u,caddr cadr v) = 0 % denominator. and lowertowerp1(cadar u,car v) then car u else lowertowerp(cdr u,v); symbolic procedure lowertowerp1(u,v); % This procedure decides if u can be an algebraic extension of v. % The = case is decidedly heuristic at the moment. % We could think of this as a membership test (including =). % However, different SQRT representations complicate things. (if x>y then t else if numberp u and numberp v then not(gcdn(u,v)=1) else x=y) where x=exprsize u,y=exprsize v; symbolic procedure exprsize u; % Get size of u. Iterative to avoid excessive recursion. begin integer n; a: if null u then return n else if atom u then return n+1; n := exprsize car u + n; u := cdr u; go to a end; symbolic procedure rationalizef u; % Look for I and sqrts, cbrts, quartics at present. % I'm not sure I in the presence of (-1)^(1/4) say is handled % properly. % It is assumed that any surd powers have been reduced before % entering this procedure. begin scalar x,y,z; x := z := kernels u; a: if null x then return 1; y := car x; if eqcar(y,'expt) and eqcar(caddr y,'quotient) and lowertowerp(z,cdr y) then nil else if y eq 'i or eqcar(y,'expt) and caddr y = '(quotient 1 2) or eqcar(y,'sqrt) then return conjquadratic(mkmain(u,y),y) else if eqcar(y,'expt) and caddr y = '(quotient 1 3) then return conjcubic(mkmain(u,y),y) else if eqcar(y,'expt) and caddr y = '(quotient 1 4) then return conjquartic(mkmain(u,y),y); x := cdr x; go to a end; symbolic procedure conjquadratic(u,v); if ldeg u = 1 then subtrf(multf(!*k2f v,reorder lc u),reorder red u) else errach list(ldeg u,"invalid power in rationalizef"); symbolic procedure conjcubic(u,v); begin scalar c1,c2,c3,w; if ldeg u = 2 then <<c1 := reorder lc u; if degr(red u,v) = 1 then <<c2 := reorder lc red u; c3 := reorder red red u>> else c3 := reorder red u>> else <<c2 := reorder lc u; c3 := reorder red u>>; w := conj2 v; if w eq 'failed then return u; v := !*k2f v; return addf(multf(exptf(v,2),subtrf(exptf(c2,2),multf(c1,c3))), addf(multf(v,subtrf(multf(w,exptf(c1,2)), multf(c2,c3))), subtrf(exptf(c3,2),multf(w,multf(c1,c2))))) end; symbolic procedure conj2 u; % (if not domainp denr v then errach list("conj2",u) (if not domainp denr v then 'failed else if denr v neq 1 then multd(!:recip denr v,numr v) else numr v) where v = simp cadr u; symbolic procedure conjquartic(u,v); begin scalar c1,c3,c4,q1,q2,q3,q4,w; if ldeg u = 3 then <<c1 := reorder lc u; if degr(red u,v) = 1 then <<c3 := reorder lc red u; c4 := reorder red red u>> else c4 := reorder red u>> else if ldeg u = 1 then <<c3 := reorder lc u; c4 := reorder red u>>; w := conj2 v; if w eq 'failed then return u; v := !*k2f v; q1 := subtrf(addf(exptf(c3,3),multf(c1,exptf(c4,2))), multf(w,multf(c3,exptf(c1,2)))); q2 := negf addf(multf(w,multf(c4,exptf(c1,2))), multf(exptf(c3,2),c4)); q3 := addf(multf(c3,exptf(c4,2)), subtrf(multf(exptf(w,2),exptf(c1,3)), multf(w,multf(c1,exptf(c3,2))))); q4 := subtrf(multf(w,multf(multd(2,c1),multf(c3,c4))),exptf(c4,3)); return addf(multf(exptf(v,3),q1), addf(multf(exptf(v,2),q2),addf(multf(v,q3),q4))) end; symbolic procedure mkmain(u,var); % Make kernel var the main variable of u. begin scalar kord!*; kord!* := list var; return reorder u end; endmodule; end;