reduce-historical/r37/packages/support/psl.red

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;