Artifact 5f5139231213ac406ace3804c8c9daa9549341006a5d66420dbda994663eef3c:
- Executable file
r37/packages/support/psl.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: 12605) [annotate] [blame] [check-ins using] [more...]
module psl; imports big2sys, bigp, floatloworder, floathighorder, gtneg, gtpos, i2bf!:, idifference, igetv, ilessp, iminus, inf, iplus, isub1, itimes, land, lshift, make!:ibf, neq, sys2int, trimbignum, vecinf, veclen, wand, wdifference, wminus, wor, wplus2, wputv, wquotient, wshift; exports ashift, msd!:, fl2bf, integerp!:, normbf, oddintp, preci!:; fluid '(bbits!*); global '(bfz!*); compiletime global '(!!fleps1exp !!plumaxexp !!pluminexp !!timmaxexp !!timminexp); remflag ('(ashift msd!: fl2bf ff0 ff1 bf!-bits bf!-bits!-mask integerp!: normbf oddintp preci!:), 'lose); flag('(cond),'eval); % Enable conditional compilation. %------------------------------------------------------------------- % The following routines support fast float operations by exploiting % the IEEE number format explicitly. compiletime if 'ieee member lispsystem!* then remflag('(safe!-fp!-plus safe!-fp!-times safe!-fp!-quot),'lose) else flag('(safe!-fp!-plus safe!-fp!-times safe!-fp!-quot),'lose); % Currently 32 and 64 bit IEEE machines are supported. % % The following macros assume that on 64 bit machines floathighorder % and floatloworder both load the full 64 bit floating point number. compiletime << define!-constant(ieeeshift,12 - bitsperword); % 32 bits:-20 define!-constant(signshift,1 - bitsperword); % 32 bits:-31 define!-constant(ieeebias,1023); define!-constant(ieeemask,2047); ds(floathiword,x(),floathighorder inf x); ds(floatloword,x(),floatloworder inf x); if bitsperword=32 then << ds(ieeezerop,u(), weq(0,floathiword u) and weq(0,floatloword u)); ds(ieeeequal,u(v), weq(floathiword u,floathiword v) and weq(floatloword u,floatloword v)); ds(ieeemant,f(), (lor(lshift( wor(wshift(wand (floathiword f, 1048575), % 16#FFFFF 6), wshift(lf,-26)), 26), wand(lshift(-1,-6), lf)) where lf := floatloword f)); >> else if bitsperword=64 then << ds(ieeezerop,u(), weq(0,floathiword u)); ds(ieeeequal,u(v), weq(floathiword u,floathiword v)); ds(ieeemant,f(), wand (floathiword f, 4503599627370495)); % 16#FFFFFFFFFFFFF >> else error(99,"#### unknown bit size"); ds(ieeeexpt,u(), wdifference(wand(ieeemask, wshift(floathiword u,ieeeshift)), ieeebias)); ds(ieeesign,u(),wshift(floathiword u,signshift)); % ieeemant is the mantissa part of the upper 32 bit group. define!-constant(!!plumaxexp,1018); define!-constant(!!pluminexp,-979); define!-constant(!!timmaxexp,509); define!-constant(!!timminexp,-510); define!-constant(!!fleps1exp,-40) >>; symbolic procedure safe!-fp!-plus(x,y); if ieeezerop x then y else if ieeezerop y then x else begin scalar u,ex,ey,sx,sy; ex := ieeeexpt x; ey := ieeeexpt y; if (sx := ieeesign x) eq (sy := ieeesign y) then if ilessp(ex,!!plumaxexp) and ilessp(ey,!!plumaxexp) then go to ret else return nil; if ilessp(ex,!!pluminexp) and ilessp(ey,!!pluminexp) then return nil; ret: u := floatplus2(x,y); return if sx eq sy or ieeezerop u then u else if ilessp(ieeeexpt u,iplus2(ex,!!fleps1exp)) then 0.0 else u end; symbolic procedure safe!-fp!-times(x,y); if ieeezerop x or ieeezerop y then 0.0 else if ieeeequal(x,1.0) then y else if ieeeequal(y,1.0) then x else begin scalar u,v; u := ieeeexpt x; v := ieeeexpt y; if igreaterp(u,!!timmaxexp) then if ilessp(v,0) then go to ret else return nil; if igreaterp(u,0) then if ilessp(v,!!timmaxexp) then go to ret else return nil; if igreaterp(u,!!timminexp) then if igreaterp(v,!!timminexp) then go to ret else return nil; if ilessp(v,0) then return nil; ret: return floattimes2(x,y) end; symbolic procedure safe!-fp!-quot(x,y); if ieeezerop y then rdqoterr() else if ieeezerop x then 0.0 else if ieeeequal(y,1.0) then x else begin scalar u,v; u := ieeeexpt x; v := ieeeexpt y; if igreaterp(u,!!timmaxexp) then if igreaterp(v,0) then go to ret else return nil; if igreaterp(u,0) then if igreaterp(v,!!timminexp) then go to ret else return nil; if igreaterp(u,!!timminexp) then if ilessp(v,!!timmaxexp) then go to ret else return nil; if igreaterp(v,0) then return nil; ret: return floatquotient(x,y) end; compiletime if 'ieee member lispsystem!* then flag('(safe!-fp!-plus safe!-fp!-times safe!-fp!-quot),'lose) else remflag('(safe!-fp!-plus safe!-fp!-times safe!-fp!-quot),'lose); %--------------------------------------------------------------- deflist('((iminus iminus)),'unary); symbolic smacro procedure ashift (m,d); if negintp m then -lshift(-m,d) else lshift(m,d); symbolic smacro procedure oddintp x; wand(if bigp x then wgetv(inf x,2) else if fixnp x then fixval inf x else x,1) eq 1; symbolic macro procedure bf!-bits (x); {'quote, bbits!*}; %symbolic macro procedure bf!-bits!-mask (x); % {'quote, lshift(1, bf!-bits()) - 1}; %symbolic procedure ff1 (w,n); % if n eq 0 then w else % if wshift (w, wminus n) eq 0 then % ff1 (w,wquotient(n,2)) % else iplus2(ff1 (wshift (w, wminus n),wquotient(n,2)),n) ; symbolic smacro procedure ff1 (ww,nn); <<while not (n eq 0) do << x := wshift(w,wminus n); if not (x eq 0) then % left half << w := x; y := iplus2(y,n) >>; % Iplus2 etc. used for n := wquotient (n,2) % bootstrapping. >>; iplus2(y,w) >> where w=ww,n=nn,x=nil,y=0; %symbolic procedure ff0 (w,n); %% returns the number of 0 bits at the least significant end % if n eq 0 then w else % begin scalar lo; % lo := wand(w,isub1 wshift(1,n)); % return if lo eq 0 % then iplus2(n,ff0 (wshift(w,wminus n),wquotient(n,2))) % else ff0 (lo,wquotient(n,2)) ; % end; comment ff0 determines the number of 0 bits at the least significant end of an integer, ie. the largest power of two by which the integer is divisible; compiletime put('hu_hu_hu,'opencode,'((!*move (reg 1) (reg 1)))); symbolic smacro procedure ff0 (ww,nn); <<while not (n eq 0) do << lo := wand(w,isub1 wshift(1,n)); if lo eq 0 then % left half << w := wshift(w,wminus n); y := iplus2(y,n) >>; % Iplus2 etc. used for n := wquotient (n,2) % bootstrapping. >>; if w neq 0 then << w := 17; hu_hu_hu (w); y >> else iadd1 y >> % we have to destroy w for gc !! where w=ww,n=nn,lo=nil,y=0; % use wshift(bitsperword,-1) rather than bitsperword/2 as the former % is open compiled comment we split msd!: into two parts: one for bignums, one for machine words. That will greatly reduce the size of preci!: below; symbolic smacro procedure word!-msd!: u; ff1(u,wshift(bitsperword,-1)); symbolic smacro procedure big!-msd!: u; iplus2(itimes2(bf!-bits(),isub1 s),word!-msd!: igetv(u,s)) where s := veclen vecinf u; symbolic smacro procedure msd!: u; if bigp u then big!-msd!: u else if fixnp u then word!-msd!: fixval inf u else word!-msd!: u; %symbolic smacro procedure msd!: u; % % returns the most significant (binary) digit of a positive integer u % if bigp u % then iplus2(itimes2(bf!-bits(),isub1 s), % ff1(igetv(u,s),wshift(bitsperword,-1))) % where s := veclen vecinf u % else if fixnp u then ff1 (fixval inf u,wshift(bitsperword,-1)) % else ff1 (u,wshift(bitsperword,-1)); symbolic smacro procedure mt!: u; cadr u; symbolic smacro procedure ep!: u; cddr u; symbolic smacro procedure preci!: nmbr; % This function counts the precision of a number "n". NMBR is a % binary bigfloat representation of "n". % msd!: abs mt!: nmbr (if bigp m then big!-msd!: m else if fixnp m then (word!-msd!:(if iminusp n then iminus n else n) where n = fixval inf m) else if iminusp m then word!-msd!:(iminus m) else word!-msd!: m) where m = mt!: nmbr; %symbolic smacro procedure preci!: nmbr; % % This function counts the precision of a number "n". NMBR is a % % binary bigfloat representation of "n". % % msd!: abs mt!: nmbr % (if bigp m then msd!: m % else if fixnp m % then (ff1(if iminusp n then iminus n else n, % wshift(bitsperword,-1)) % where n = fixval inf m) % else if iminusp m then ff1(iminus m,wshift(bitsperword,-1)) % else ff1(m,wshift(bitsperword,-1))) % where m = mt!: nmbr; symbolic smacro procedure make!:ibf (mt, ep); '!:rd!: . (mt . ep); if not('ieee memq lispsystem!*) then flag('(fl2bf),'lose); symbolic procedure fl2bf f; % u is a floating point number % result is a binary bigfloat if fixp f then i2bf!: f else begin scalar m,e; m := ieeemant f; e := ieeeexpt f; % if exponent <> -1023 add 16#10000000000000, implicit highest bit if e neq -1023 then m := lor (m, lshift(1,52)); return if izerop m then bfz!* else normbf make!:ibf (if ieeesign f eq 1 then -m else m, idifference(e,52)) end; symbolic procedure normbf x; begin scalar mt,s;integer ep,ep1; if (mt := mt!: x)=0 then go to ret; if mt<0 then <<mt := -mt; s := t>>; ep := ep!: x; % ep1 := remainder(ep,bf!-bits()); % if ep1 < 0 then ep1 := ep1 + bf!-bits(); % if ep1 neq 0 then <<ep := ep - ep1; mt := lshift(mt,ep1)>>; while bigp mt and wgetv(inf mt,2) eq 0 do << mt := lshift(mt,-bf!-bits()); ep := ep+bf!-bits() >>; ep1 := ff0(if bigp mt then wgetv(inf mt,2) else if fixnp mt then fixval inf mt else mt,wshift(bitsperword,-1)); if not (ep1 eq 0) then <<mt := lshift(mt,wminus ep1); ep := ep + ep1>>; if s then mt := -mt; ret: return make!:ibf(mt,ep) end; %symbolic procedure normbf x; % begin scalar mt,s;integer ep,ep1; % if (mt := mt!: x)=0 then go to ret; % if mt<0 then <<mt := -mt; s := t>>; % ep := ep!: x; % while bigp mt and land(mt,bf!-bits!-mask())=0 do << % mt := lshift(mt,-bf!-bits()); % ep := ep+bf!-bits() >>; % while land(mt,255)=0 do << % mt := lshift(mt,-8); % ep := ep+8 >>; % while land(mt,1)=0 do << % mt := lshift(mt,-1); % ep := ep+1>>; %% ep1 := remainder(ep,bf!-bits()); %% if ep1 < 0 then ep1 := ep1 + bf!-bits(); %% if ep1 neq 0 then <<ep := ep - ep1; mt := lshift(mt,ep1)>>; % if s then mt := -mt; %ret: return make!:ibf(mt,ep) end; symbolic procedure integerp!: x; % This function returns T if X is a BINARY BIG-FLOAT % representing an integer, else it returns NIL. % X is any LISP entity. bfp!: x and (ep!: x >= 0 or preci!: x > - ep!: x and land(abs mt!: x,lshift(2,-ep!: x) - 1) = 0); flag ('(ashift lshift msd!: fl2bf ff0 ff1 bf!-bits bf!-bits!-mask integerp!: normbf oddintp preci!:), 'lose); if not('ieee memq lispsystem!*) then remflag('(fl2bf),'lose); % This belong in $pxu/nbig30a. symbolic(bigfloathi!* := (2 ** 53 - 1) * 2 ** 971); symbolic(bigfloatlow!* := - bigfloathi!*); remflag('(cond),'eval); % HP-Risc and IBM RS architectures need special handling of fltinf in % fastmath.red if 'HP!-Risc member lispsystem!* then <<remflag('(fltinf),'lose); ds(fltinf,x(),mkitem(vector!-tag,x)); flag('(fltinf),'lose)>>; if 'IBMRS member lispsystem!* then <<remflag('(fltinf),'lose); ds(fltinf,x(),mkstr x); flag('(fltinf),'lose)>>; endmodule; end;