Artifact a59144d60834c1dfa6518f51649fb1c6de4005b0ec9ba8e77aca6d6c3118a272:
- Executable file
r37/packages/int/jpatches.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: 10987) [annotate] [blame] [check-ins using] [more...]
- Executable file
r38/packages/int/jpatches.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: 10987) [annotate] [blame] [check-ins using]
module jpatches; % Routines for manipulating sf's with power folding. % Author: James H. Davenport. fluid '(!*noncomp sqrtflag); exports !*addsq,!*multsq,!*invsq,!*multf,!*exptsq,!*exptf; %squashsqrtsq % !*MULTF(A,B) multiplies the polynomials (standard forms) U and V % in much the same way as MULTF(U,V) would, EXCEPT... % (1) !*MULTF inhibits the action of OFF EXP and of non-commutative % multiplications % (2) Within !*MULTF powers of square roots, and powers of % exponential kernels are reduced as if substitution rules % such as FOR ALL X LET SQRT(X)**2=X were being applied; % Note that !*MULTF comes between MULTF and !*Q2F SUBS2F MULTF in its % behaviour, and that it is the responsibility of the user to call it % in sensible places where its services are needed. % Similarly for the other functions defined here. %symbolic procedure !*addsq(u,v); %U and V are standard quotients. % %Value is canonical sum of U and V; % if null numr u then v % else if null numr v then u % else if denr u=1 and denr v=1 then addf(numr u,numr v) ./ 1 % else begin scalar nu,du,nv,dv,x; % x := gcdf(du:=denr u,dv:=denr v); % du:=quotf(du,x); dv:=quotf(dv,x); % nu:=numr u; nv:=numr v; % u:=addf(!*multf(nu,dv),!*multf(nv,du)); % if u=nil then return nil ./ 1; % v:=!*multf(du,denr v); % return !*ff2sq(u,v) % end; %symbolic procedure !*multsq(a,b); % begin % scalar n,d; % n:=!*multf(numr a,numr b); % d:=!*multf(denr a,denr b); % return !*ff2sq(n,d) % end; %symbolic procedure !*ff2sq(a,b); % begin % scalar gg; % if null a then return nil ./ 1; % gg:=gcdf(a,b); % if not (gg=1) then << % a:=quotf(a,gg); % b:=quotf(b,gg) >>; % if minusf b then << % a:=negf a; % b:=negf b >>; % return a ./ b % end; symbolic procedure !*addsq(u,v); %U and V are standard quotients. %Value is canonical sum of U and V; if null numr u then v else if null numr v then u else if denr u=1 and denr v=1 then addf(numr u,numr v) ./ 1 else begin scalar du,dv,x,y,z; x := gcdf(du:=denr u,dv:=denr v); du:=quotf(du,x); dv:=quotf(dv,x); y:=addf(!*multf(dv,numr u),!*multf(du,numr v)); if null y then return nil ./ 1; z:=!*multf(denr u,dv); if minusf z then <<y := negf y; z := negf z>>; % In this case (as opposed to ADDSQ), Y and Z may have % developed common factors from SQRT expansion, so a % gcd of Y and Z is needed. x := gcdf(y,z); return if x=1 then y ./ z else quotf(y,x) ./ quotf(z,x) end; symbolic procedure !*multsq(u,v); %U and V are standard quotients. Result is the canonical product of %U and V with surd powers suitably reduced. if null numr u or null numr v then nil ./ 1 else if denr u=1 and denr v=1 then !*multf(numr u,numr v) ./ 1 else begin scalar w,x,y; x := gcdf(numr u,denr v); y := gcdf(numr v,denr u); w := !*multf(quotf(numr u,x),quotf(numr v,y)); x := !*multf(quotf(denr u,y),quotf(denr v,x)); if minusf x then <<w := negf w; x := negf x>>; y := gcdf(w,x); % another factor may have been generated. return if y=1 then w ./ x else quotf(w,y) ./ quotf(x,y) end; symbolic procedure !*invsq a; % Note that several examples (e.g., int(1/(x**8+1),x)) give a more % compact result when SQRTFLAG is true if SQRT2TOP is not called. if sqrtflag then sqrt2top invsq a else invsq a; symbolic procedure !*multf(u,v); % U and V are standard forms % Value is SF for U*V; begin scalar x,y; if null u or null v then return nil else if u=1 then return squashsqrt v else if v=1 then return squashsqrt u else if domainp u then return multd(u,squashsqrt v) else if domainp v then return multd(v,squashsqrt u) else if !*noncomp then return multf(u,v); x:=mvar u; y:=mvar v; if x eq y then go to c else if ordop(x,y) then go to b; x:=!*multf(u,lc v); y:=!*multf(u,red v); return if null x then y else if not domainp lc v and mvar u eq mvar lc v and not atom mvar u and car mvar u memq '(expt sqrt) then addf(!*multf(x,!*p2f lpow v),y) else makeupsf(lpow v,x,y); b: x:=!*multf(lc u,v); y:=!*multf(red u,v); return if null x then y else if not domainp lc u and mvar lc u eq mvar v and not atom mvar v and car mvar v memq '(expt sqrt) then addf(!*multf(!*p2f lpow u,x),y) else makeupsf(lpow u,x,y); c: y:=addf(!*multf(list lt u,red v),!*multf(red u,v)); if eqcar(x,'sqrt) then return addf(squashsqrt y,!*multfsqrt(x, !*multf(lc u,lc v),ldeg u + ldeg v)) else if eqcar(x,'expt) and prefix!-rational!-numberp caddr x then return addf(squashsqrt y,!*multfexpt(x, !*multf(lc u,lc v),ldeg u + ldeg v)); x:=mkspm(x,ldeg u + ldeg v); return if null x or null (u:=!*multf(lc u,lc v)) then y else addf(multpf(x,u),y) end; symbolic procedure makeupsf(u,x,y); % Makes u .* x .+ y except when u is not a valid lpow (because of % sqrts). if atom car u or cdr u = 1 then addf(multpf(u,x),y) else if caar u eq 'sqrt then addf(!*multfsqrt(car u,x,cdr u),y) else if <<begin scalar v; v:=car u; if car v neq 'expt then return nil; v:=caddr v; if atom v then return nil; return (car v eq 'quotient and numberp cadr v and numberp caddr v) end >> then addf(!*multfexpt(car u,x,cdr u),y) else addf(multpf(u,x),y); symbolic procedure !*multfsqrt(x,u,w); % This code (Due to Norman a& Davenport) squashes SQRT(...)**2. begin scalar v; w:=divide(w,2); v:=!*q2f simp cadr x; u:=!*multf(u,exptf(v,car w)); if cdr w neq 0 then u:=!*multf(u,!*p2f mksp(x,1)); return u end; symbolic procedure !*multfexpt(x,u,w); begin scalar expon,v; expon:=caddr x; x:=cadr x; w:=w * cadr expon; expon:=caddr expon; v:=gcdn(w,expon); w:=w/v; v:=expon/v; if not (w > 0) then rerror(int,8,"Invalid exponent") else if v = 1 then return !*multf(u,exptf(if numberp x then x else if atom x then !*k2f x else !*q2f if car x eq '!*sq then argof x else simp x, w)); expon:=0; while not (w < v) do <<expon:=expon + 1; w:=w-v>>; if expon>0 then u:=!*multf(u,exptf(!*q2f simp x,expon)); if w = 0 then return u; x:=list('expt,x,list('quotient,1,v)); return multf(squashsqrt u,!*p2f mksp(x,w)) % Cannot be *MULTF. end; symbolic procedure prefix!-rational!-numberp u; % Tests for m/n in prefix representation. eqcar(u,'quotient) and numberp cadr u and numberp caddr u; % symbolic procedure squashsqrtsq sq; % !*multsq(squashsqrt numr sq ./ 1,1 ./ squashsqrt denr sq); symbolic procedure squashsqrt sf; if (not sqrtflag) or atom sf or atom mvar sf then sf else if car mvar sf eq 'sqrt and ldeg sf > 1 then addf(squashsqrt red sf,!*multfsqrt(mvar sf,lc sf,ldeg sf)) else if car mvar sf eq 'expt and prefix!-rational!-numberp caddr mvar sf and ldeg sf >= caddr caddr mvar sf then addf(squashsqrt red sf,!*multfexpt(mvar sf,lc sf,ldeg sf)) else (if null x then squashsqrt red sf else lpow sf .* x .+ squashsqrt red sf) where x = squashsqrt lc sf; %remd 'simpx1; % The following definition requires frlis!* declared global. %symbolic procedure simpx1(u,m,n); % %u,m and n are prefix expressions; % %value is the standard quotient expression for u**(m/n); % begin scalar flg,z; % if null frlis!* or null intersection(frlis!*,flatten (m . n)) % then go to a; % exptp!* := t; % return !*k2q list('expt,u,if n=1 then m % else list('quotient,m,n)); % a: if numberp m and fixp m then go to e % else if atom m then go to b % else if car m eq 'minus then go to mns % else if car m eq 'plus then go to pls % else if car m eq 'times and numberp cadr m and fixp cadr m % and numberp n % then go to tms; % b: z := 1; % c: if atom u and not numberp u then flag(list u,'used!*); % u := list('expt,u,if n=1 then m else list('quotient,m,n)); % if not(u member exptl!*) then exptl!* := u . exptl!*; % d: return mksq(u,if flg then -z else z); %u is already in lowest %% %terms; % e: if numberp n and fixp n then go to int; % z := m; % m := 1; % go to c; % mns: m := cadr m; % if !*mcd then return invsq simpx1(u,m,n); % flg := not flg; % go to a; % pls: z := 1 ./ 1; % pl1: m := cdr m; % if null m then return z; % z := multsq(simpexpt list(u, % list('quotient,if flg then list('minus,car m) % else car m,n)), % z); % go to pl1; % tms: z := gcdn(n,cadr m); % n := n/z; % z := cadr m/z; % m := retimes cddr m; % go to c; % int:z := divide(m,n); % if cdr z<0 then z:= (car z - 1) . (cdr z+n); % if 0 = cdr z % then return simpexpt list(u,car z); % if n = 2 % then return multsq(simpexpt list(u,car z), % simpsqrti u); % return multsq(simpexpt list(u,car z), % mksq(list('expt,u,list('quotient,1,n)),cdr z)) % end; symbolic procedure !*exptsq(a,n); % Raises A to the power N using !*MULTSQ. if n=0 then 1 ./ 1 else if n=1 then a else if n<0 then !*exptsq(invsq a,-n) else begin scalar q,r; q:=divide(n,2); r:=cdr q; q:=car q; q:=!*exptsq(!*multsq(a,a),q); if r=0 then return q else return !*multsq(a,q) end; symbolic procedure !*exptf(a,n); % Raises A to the power N using !*MULTF. if n=0 then 1 else if n=1 then a else begin scalar q,r; q:=divide(n,2); r:=cdr q; q:=car q; q:=!*exptf(!*multf(a,a),q); if r=0 then return q else return !*multf(a,q) end; endmodule; end;