Artifact 74106d0dc5698bb764395ca306b0bedacfc771c2a023b491abe63d0d9fc94f45:
- Executable file
r37/packages/arith/comprd.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: 9744) [annotate] [blame] [check-ins using] [more...]
- Executable file
r38/packages/arith/comprd.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: 9744) [annotate] [blame] [check-ins using]
module comprd; % *** Support for Complex Rounded Arithmetic. % Authors: Anthony C. Hearn and Stanley L. Kameny. % Last updated: 23 June 1993. % Copyright (c) 1989, 1993 The RAND Corporation. All rights reserved. Comment this module defines a complex rounded as: (<tag>. (<structure> . <structure>>) or ('!:cr!: . (rl . im)) The <structure> depends on the precision. It is either a floating point number or the stripped bfloat (mt . ep); exports !*cr2crn, !*cr2rd, !*cr2rn, !*crn2cr, !*gi2cr, !*rd2cr, !*rn2cr, cr!:differ, cr!:minus, cr!:minusp, cr!:onep, cr!:plus, cr!:prep, cr!:prin, cr!:quotient, cr!:times, cr!:zerop, cr2i!*, crhalf!*, cri!*, cri!/2, crone!*, crprcd, gf2cr!:, i2cr!*, mkcr; imports bfloat, bfnzp, bftrim!:, bfzp, chkint!*, chkrn!*, convprec, convprec!*, crim, crrl, ep!:, errorp, errorset!*, gf2bf, gfdiffer, gfminus, gfplus, gfquotient, gftimes, gfzerop, initdmode, leq, lessp!:, make!:cr, make!:rd, maprin, mkcrn, mkquote, mkrn, mkround, normbf, over, plubf, preci!:, prin2!*, r2bf, rd!:minus, rd!:minusp, rd!:onep, rd!:prep, rd!:zerop, rdprep1, rdhalf!*, rdone!*, rdqoterr, rdtwo!*, rdzchk, rdzero!*, realrat, retag, rndbfon, round!:mt, safe!-fp!-plus, safe!-fp!-times, timbf, union; fluid '(!:prec!: !:bprec!:); global '(bfone!* epsqrt!*); fluid '(dmode!* !*bfspace !*numval !*roundbf !*!*roundbf); global '(cr!-tolerance!* domainlist!* !!nfpd !!flprec !!rdprec bfz!* yy!!); domainlist!* := union('(!:cr!:),domainlist!*); fluid '(!*complex!-rounded); put('complex!-rounded,'tag,'!:cr!:); put('!:cr!:,'dname,'complex!-rounded); flag('(!:cr!:),'field); put('!:cr!:,'i2d,'i2cr!*); put('!:cr!:,'plus,'cr!:plus); put('!:cr!:,'times,'cr!:times); put('!:cr!:,'difference,'cr!:differ); put('!:cr!:,'quotient,'cr!:quotient); put('!:cr!:,'zerop,'cr!:zerop); put('!:cr!:,'onep,'cr!:onep); put('!:cr!:,'prepfn,'cr!:prep); put('!:cr!:,'prifn,'cr!:prin); put('!:cr!:,'minus,'cr!:minus); put('!:cr!:,'minusp,'cr!:minusp); % put('!:cr!:,'rationalizefn,'girationalize!:); % Needs something % different. put('!:cr!:,'!:rn!:,'!*cr2rn); put('!:rn!:,'!:cr!:,'!*rn2cr); put('!:rd!:,'!:cr!:,'!*rd2cr); put('!:cr!:,'!:rd!:,'!*cr2rd); put('!:cr!:,'!:crn!:,'!*cr2crn); put('!:crn!:,'!:cr!:,'!*crn2cr); put('!:gi!:,'!:cr!:,'!*gi2cr); put('!:cr!:,'cmpxfn,'mkcr); put('!:cr!:,'ivalue,'mkdcrn); put('!:cr!:,'realtype,'!:rd!:); put('!:rd!:,'cmpxtype,'!:cr!:); symbolic procedure cr!:minusp u; (if atom x then zerop y and x<0 else zerop car y and car x<0) where x=cadr u,y=cddr u; symbolic procedure striptag u; if atom u then u else cdr u; symbolic procedure mkcr(u,v); make!:cr (striptag u, striptag v); symbolic procedure gf2cr!: x; make!:cr (striptag car x, striptag cdr x); symbolic procedure crprcd u; (rl . im) where rl=convprec!* crrl u,im=convprec!* crim u; symbolic procedure crprcd2(x,y); <<x := crprcd x; yy!! := crprcd y; x>>; % simp must call convprec!*, since precision may have changed. symbolic procedure cr!:simp u; (gf2cr!: crprcd u) ./ 1; put('!:cr!:,'simpfn,'cr!:simp); %symbolic procedure mkdcr u; cri!*() ./ 1; symbolic procedure i2cr!* u; %converts integer U to tagged cr form. <<u := chkint!* u; mkcr(u,if atom u then 0.0 else bfz!*)>>; symbolic procedure trimcrrl n; trimcr crrl n; symbolic procedure trimcr n; bftrim!: if atom n then bfloat n else retag n; symbolic procedure cr2rderr; error(0, "complex to real type conversion requires zero imaginary part"); symbolic procedure !*cr2rn n; % Converts a cr number n into a rational if possible. if bfnzp retag crim n then cr2rderr() else <<n := realrat trimcrrl n; mkrn(car n,cdr n)>>; symbolic procedure !*rn2cr u; % Converts the (tagged) rational u/v into a (tagged) rounded complex % number to the system precision. <<u := chkrn!* r2bf cdr u; mkcr(u,if atom u then 0.0 else bfz!*)>>; symbolic procedure !*cr2crn u; % Converts a (tagged) cr number u into a (tagged) crn. mkcrn(realrat trimcrrl u,realrat trimcr crim u); symbolic procedure !*crn2cr u; % Converts a (tagged) crn number u into a (tagged) cr. mkcr(rl,if !*roundbf then bfloat im else im) where rl=chkrn!* r2bf cadr u where im=chkrn!* r2bf cddr u; symbolic procedure !*cr2rd n; if bfnzp retag crim n then cr2rderr() else make!:rd crrl n; symbolic procedure !*rd2cr u; mkcr(x,if atom x then 0.0 else bfz!*) where x=convprec u; symbolic procedure !*gi2cr u; mkcr(rl,if !*roundbf then bfloat im else im) where rl=chkint!* cadr u where im=chkint!* cddr u; symbolic procedure bfrsq u; (if atom x then x*x+y*y else plubf(timbf(x,x),timbf(y,y))) where x=car u,y=cdr u; symbolic procedure crzchk(u,x,y); begin if atom car u then if bfrsq u<(bfrsq x)*!!fleps2 then return 0.0 . 0.0 else go to ck; if lessp!:(bfrsq u,timbf(bfrsq x,cr!-tolerance!*)) then return bfz!* . bfz!*; ck: return rdzchk(car u,car x,car y) . rdzchk(cdr u,cdr x,cdr y) end; symbolic procedure cr!:plus(u,v); begin scalar x,y; x := crprcd2(u,v); y := yy!!; u := if !*!*roundbf then gfplus(x,y) else if (v := safe!-crfp!-plus(x,y)) then v else ((if errorp r then <<rndbfon(); gfplus(x := gf2bf x,y := gf2bf y)>> else car r) where r=errorset(list('gfplus,mkquote x,mkquote y),nil,nil)); return gf2cr!: crzchk(u,x,y) end; symbolic procedure cr!:differ(u,v); begin scalar x,y; x := crprcd2(u,v); y := yy!!; u := if !*!*roundbf then gfdiffer(x,y) else if (v := safe!-crfp!-diff(x,y)) then v else ((if errorp r then <<rndbfon(); gfplus(x := gf2bf x,y := gf2bf y)>> else car r) where r=errorset(list('gfdiffer,mkquote x,mkquote y),nil,nil)); return gf2cr!: crzchk(u,x,gfminus y) end; symbolic procedure cr!:times(u,v); gf2cr!: (if !*!*roundbf then gftimes(x,yy!!) else if (u := safe!-crfp!-times(x,yy!!)) then u else ((if errorp r then <<rndbfon(); gftimes(gf2bf x,gf2bf yy!!)>> else car r) where r=errorset!*(list('gftimes,mkquote x,mkquote yy!!),nil))) where x=crprcd2(u,v); symbolic procedure cr!:quotient(u,v); gf2cr!: (if gfzerop yy!! then rdqoterr() else if !*!*roundbf then gfquotient(x,yy!!) else if (u := safe!-crfp!-quot(x,yy!!)) then u else ((if errorp r then <<rndbfon(); gfquotient(gf2bf x,gf2bf yy!!)>> else car r) where r=errorset!*(list('gfquotient,mkquote x,mkquote yy!!),nil))) where x=crprcd2(u,v); symbolic procedure safe!-crfp!-plus(u,v); (if x and y then crzchk(x . y,u,v)) where x=safe!-fp!-plus(car u,car v),y=safe!-fp!-plus(cdr u,cdr v); symbolic procedure safe!-crfp!-diff(u,v); (if x and y then crzchk(x . y,u,gfminus v)) where x=safe!-fp!-plus(car u,-car v),y=safe!-fp!-plus(cdr u,-cdr v); symbolic procedure safe!-crfp!-times(u,v); begin scalar ru,iu,rv,iv,a,b; ru := car u; iu := cdr u; rv := car v; iv := cdr v; if not (a := safe!-fp!-times(ru,rv)) or not (b := safe!-fp!-times(iu,iv)) then return nil; if not(u := safe!-fp!-plus(a,-b)) then return nil; u := rdzchk(u,a,-b); if not (a := safe!-fp!-times(ru,iv)) or not (b := safe!-fp!-times(iu,rv)) then return nil; if not(v := safe!-fp!-plus(a,b)) then return nil; return u . rdzchk(v,a,b) end; symbolic procedure safe!-crfp!-quot(u,v); % compute u * inverse v. begin scalar ru,iu,rv,iv,a,b,dd; ru := car u; iu := cdr u; rv := car v; iv := cdr v; if not (a := safe!-fp!-times(rv,rv)) or not (b := safe!-fp!-times(iv,iv)) or not (dd := safe!-fp!-plus(a,b)) then return nil; rv := rv/dd; iv := iv/dd; if not (a := safe!-fp!-times(ru,rv)) or not (b := safe!-fp!-times(iu,iv)) or not (u := safe!-fp!-plus(a,b)) then return nil; u := rdzchk(u,a,b); if not (a := safe!-fp!-times(ru,-iv)) or not (b := safe!-fp!-times(iu,rv)) or not (v := safe!-fp!-plus(a,b)) then return nil; return u . rdzchk(v,a,b) end; symbolic procedure cr!:minus u; gf2cr!: gfminus crprcd u; symbolic procedure cr!:zerop u; bfzp retag crrl u and bfzp retag crim u; symbolic procedure cr!:onep u; bfzp retag crim u and rd!:onep mkround retag crrl u; % prep works entirely in bfloat, to avoid floating point conversion % errors. symbolic procedure cr!:prep u; crprep1((rd!:prep tagrl u) . rd!:prep tagim u); symbolic procedure crprep1 u; % a and d are 1,-1,or rounded. (if not numberp d and rd!:zerop d then a else <<d := if numberp d or not rd!:minusp d then crprimp d else {'minus,crprimp rd!:minus d}; if not numberp a and rd!:zerop a then d else <<if a = -1 then a := {'minus,1} else if not numberp a and rd!:minusp a then a := {'minus,rd!:minus a}; {'plus,a,d}>> >>) where a=car u,d=cdr u; symbolic procedure crprimp u; if u=1 then 'i else if u= -1 then {'minus,'i} else {'times,u,'i}; symbolic procedure cr!:prin v; if atom (v := cr!:prep v) or car v eq 'times or car v memq domainlist!* then maprin v else <<prin2!* "("; maprin v; prin2!* ")">>; initdmode 'complex!-rounded; symbolic procedure crone!*; mkcr(rdone!*(),rdzero!*()); symbolic procedure crhalf!*; mkcr(rdhalf!*(),rdzero!*()); symbolic procedure cri!*; mkcr(rdzero!*(),rdone!*()); symbolic procedure cri!/2; mkcr(rdzero!*(),rdhalf!*()); symbolic procedure cr2i!*; mkcr(rdzero!*(),rdtwo!*()); endmodule; end;