Check-ins Using Download Hex

Artifact 9addf8ef5cf0e8cab979b7e5a7847f7b8ec1477ff1fceac89f5767547d57751d:

module bfauxil; % Support for the roots package and ROUNDED domain.

% Author: Stanley L. Kameny <valley!stan@rand.org>.
% Definitions of ilog2, irootn, icbrt, isqrt and support supplied by
% John Abbott.

% Copyright (c) 1988,1989,1990. Stanley L. Kameny. All Rights Reserved.

Comment   support for modules allroot and isoroot, and for ROUNDED
 domain logic;

exports !!shbinflp, bf2flr, bfdiffer, bfdivide, bfinverse, bflessp,
        bfminus, bfsqrt, cflot, conv!:bf2i, difbf, exptbf, fl2int,
        gf2bf, gf2fl, gfdiffer, gfdot, gfminus, gfplus, gfquotient,
        gfsqrt, gftimes, grpbf, icbrt, ilog2, invbf, irootn, isqrt,
        normbf, plubf, r2bf, r2fl, realrat;

imports abs!:, ashift, bflerrmsg, bfloat, bfminusp, bfnzp, bfp!:,
        bfzerop!:, bfzp, conv!:mt, cut!:ep, decprec!:, difference!:,
        divbf, divide!:, ep!:, eqcar, error1, errorp, errorset!*,
        evenp, fl2bf, gcdn, geq, gfim, gfrl, gfzerop, greaterp!:,
        hypot, i2bf!:, leq, lshift, make!:ibf, minus!:, minusp!:,
        msd!:, mt!:, order!:, plus!:, preci!:, read!:num, rndpwr,
        round!:mt, sgn, sqrt, terrlst, timbf, times!:, typerr;

fluid '(!:prec!: !:bprec!:);

global '(bfone!* bfhalf!* bfz!*);

global '(!!nfpd !!nbfpd !!shbinfl vv!! !!flbint);


global '(!!minflbf !!maxflbf);

symbolic procedure normbf x;
   begin scalar mt,s,r;integer ep,ep1;
      if (mt := mt!: x)=0 then go to ret;
      if mt<0 then <<mt := -mt; s := t>>;
      ep := ep!: x;
      while remainder(mt,1073741824)=0 do << % 2**30
        mt := lshift(mt,-30);
        ep := ep+30 >>;
      while remainder(mt,256)=0 do <<
        mt := lshift(mt,-8);
        ep := ep+8 >>;
      while not oddintp mt do <<
        mt := lshift(mt,-1);
        ep := ep+1>>;
      if s then mt := -mt;
ret:    return make!:ibf(mt,ep) end;

%symbolic procedure divbf(u,v); normbf divide!:(u,v,!:bprec!:);

symbolic procedure bfdivide(u,v);
   if atom u then u/v else divbf(u,v);

%symbolic procedure timbf(u,v); rndpwr times!:(u,v);

symbolic procedure bftimes(u,v); if atom u then u*v else timbf(u,v);

symbolic procedure plubf(a,b);
% this function calculates the normalized rounded sum of a and b,
% but avoids generating large numbers if magnitude difference is large.
  rndpwr
    begin scalar ma,mb,ea,eb,d,ld,p;
      if (ma:=mt!: a)=0 then return b;
      if (mb:=mt!: b)=0 then return a;
      if (d := (ea := ep!: a)-(eb := ep!: b))=0
        then return make!:ibf(ma+mb,ea);
      ld := d+msd!: abs ma - msd!: abs mb;
      p := !:bprec!:+1;
      if ld>p then return a;
      if ld<-p then return b;
      if d>0 then return make!:ibf(ashift(ma,d)+mb,eb)
        else return make!:ibf(ma+ashift(mb,-d),ea) end;

symbolic procedure bfplus(u,v); if atom u then u+v else plubf(u,v);

symbolic procedure difbf(a,b);
% this function calculates the normalized rounded difference of a and b,
% but avoids generating large numbers if magnitude difference is large.
  rndpwr
    begin scalar ma,mb,ea,eb,d,ld,p;
      if (ma:=mt!: a)=0 then return minus!: b;
      if (mb:=mt!: b)=0 then return a;
      if (d := (ea := ep!: a)-(eb := ep!: b))=0
        then return make!:ibf(ma - mb,ea);
      ld := d+msd!: abs ma - msd!: abs mb;
      p := !:bprec!:+1;
      if ld>p then return a;
      if ld<-p then return minus!: b;
      if d>0 then return make!:ibf(ashift(ma,d) - mb,eb)
        else return make!:ibf(ma - ashift(mb,-d),ea) end;

symbolic procedure bfdiffer(u,v); if atom u then u - v else difbf(u,v);

symbolic procedure invbf u; divbf(bfone!*,u);

symbolic procedure bfinverse u; if atom u then 1.0/u else invbf u;

symbolic procedure bfminus u; if atom u then -u else minus!: u;

symbolic procedure bflessp(a,b); if atom a then a<b else grpbf(b,a);

symbolic procedure grpbf(a,b);
  % this function calculates a > b, but avoids generating large numbers
  % if magnitude difference is large.
   <<if ma=0 then mb<0
     else if mb=0 then ma>0
     else if ma>0 and mb<0 then t
     else if ma<0 and mb>0 then nil
     else (if do>0 then ma>0           % the case |a| > |b|
            else if do<0 then ma<0     % the case |a| < |b|
            else if de=0 then ma>mb    % exponents are the same
            else if de>0 then ashift(ma,de)>mb
            else ma>ashift(mb,-de))
          where do=order!: a - order!: b,
                de=ep!: a - ep!: b>>
    where ma=mt!: a,mb=mt!: b;

%symbolic procedure bfminusp u; if atom u then minusp u else minusp!: u;

%symbolic procedure bfzp u; if atom u then zerop u else mt!: u=0;

%symbolic procedure bf!:zerop u; if atom u then zerop u else mt!: u=0;

%symbolic procedure bfnzp u; not bfzp u;

%symbolic procedure bfloat x; if floatp x then fl2bf x else
%normbf(if atom x then if fixp x then i2bf!: x else read!:num x else x);

symbolic procedure !!shbinflp;
   begin integer n; vv!! := 9.0;
     while n<300 and not errorp errorset!*('(vv!!!*1e10),nil)
        do n := n+10;
     return n<300 end;

symbolic procedure vv!!!*1e10; vv!! := vv!!*1.0e10;

symbolic (!!shbinfl := !!shbinflp());

symbolic procedure bfsqrt x;
 % computes sqrt x by Newton's method.
  if minusp!: x then terrlst(x,'bfsqrt) else
  begin  scalar nx,dx,dc,k7,nf;
         if bfzerop!: x then return bfz!*;
         k7 := !:bprec!: + 7;
         dc := make!:ibf (1, (-k7+(order!: x + 10)/2));
         nx := if not oddintp ep!:(nx := conv!:mt(x,2))
            then make!:ibf((2+3*mt!: nx)/5, (ep!: nx/2))
            else make!:ibf((9+5*mt!: nx)/10, ((ep!: nx - 1)/2));
         nf := 1;
   loop: if (nf := 2*nf)>k7 then nf := k7;
         dx := times!:(bfhalf!*,plus!:(divide!:(x,nx,nf),nx));
         if nf>=k7 and not greaterp!:(abs!: difference!:(dx,nx),dc)
            then return rndpwr nx;
         nx := dx; go to loop end;

symbolic procedure realrat x;
   begin scalar d,g;
         if bfp!: x then go to bf;
         if eqcar(x,'quotient) then
            if fixp cadr x and fixp caddr x then
               <<x := if (d := caddr x)<0 then -cadr x else cadr x;
                 d := abs d; go to ret>>
            else x := cadr x/caddr x;
         if zerop x then return (0 . 1);
         if not floatp x then return (x . 1);
         x := bfloat x;
     bf: d := cddr(x := normbf x); x := cadr x;
         if x=0 then return (0 . 1);
         if d< 0 then d := lshift(1,-d)
          else <<x := ashift(x,d); d := 1>>;
    ret: g := gcdn(abs x,d); return (x/g) . (d/g) end;

remflag ('(fl2int),'lose);

symbolic procedure fl2int x;
   <<x := fl2bf x;
     (if d=0 then m else ashift(m,d))
      where m=mt!: x,d=ep!: x>>;

flag ('(fl2int),'lose);

symbolic procedure cflot x;
   if floatp x then x else if atom x then float x else bf2flr x;

symbolic procedure conv!:bf2i nmbr;
% This function converts a <BIG-FLOAT>, i.e., a BINARY BIG-FLOAT
%      representation of "n", to an integer.  The result
%      is the integer part of "n".
% **** For getting the nearest integer to "n", please use
% ****      the combination MT!:( CONV!:EP(NMBR,0)).
% NMBR is a BIG-FLOAT representation of the number "n".
%   if ep!:(nmbr := cut!:ep(nmbr, 0)) = 0 then mt!: nmbr
%    else
     ashift (mt!: nmbr, ep!: nmbr);

symbolic procedure bf2flr u;
   % u is always bigfloat.
   % Converts bfloat to float by rounding at !!nbfpd binary digits.
   % We use error1 rather than rerror, because we want to catch such an
   % error in an errorset.
   begin scalar ep,m,y;
      if bfzerop!: u then return 0.0;
      ep := ep!:(u := round!:mt(u,!!nbfpd));
      if grpbf(!!minflbf,y := abs!: u) or grpbf(y,!!maxflbf)
        then error1();
      if ep<0 then <<ep := ep+!!nbfpd; m := t>>;
      ep := 2.0**ep;
      if ep = 0.0 then error1(); % underflow
      return if not m then ep * mt!: u else ep * mt!: u / !!flbint
   end;

symbolic procedure gf2fl a; % force into float format.
   if atom a then a else if bfp!: a then bf2flr a
      else (gf2fl car a) . gf2fl cdr a;

symbolic procedure gf2bf a; if a then % force into bfloat format.
   if atom a then bfloat a else if bfp!: a then a
      else (gf2bf car a) . gf2bf cdr a;

symbolic procedure r2bf u;
  % translate any real number object to bigfloat.
    if atom u then bfloat u
       else if bfp!: u then u
       else if numberp car u then divbf(i2bf!: car u,i2bf!: cdr u)
       else if eqcar(u,'quotient) then
            divbf(i2bf!: cadr u,i2bf!: caddr u)
       else if eqcar(u,'!:rn!:) then r2bf cdr u
       else r2bf cadr u;

symbolic procedure r2fl u;
 % translate any real number object to float.
   if u=0 then 0.0
      else if atom u then float u
      else if numberp car u then (float car u)/cdr u
      else if eqcar(u,'quotient) then (float cadr u)/caddr u
      else if bfp!: u then bf2flr u
      else if eqcar(u,'!:rn!:) then r2fl cdr u
      else r2fl cadr u;

symbolic procedure gfplus(u,v);
   if atom car u then gffplus(u,v) else gbfplus(u,v);

symbolic procedure gffplus(u,v); (car u+car v) . (cdr u+cdr v);

symbolic procedure gbfplus(u,v);
   (plubf(car u,car v)) . plubf(cdr u,cdr v);

symbolic procedure gfdiffer(u,v);
   if atom car u then gffdiff(u,v) else gbfdiff(u,v);

symbolic procedure gffdiff(u,v); (car u - car v) . (cdr u - cdr v);

symbolic procedure gbfdiff(u,v);
   (difbf(car u,car v)) . difbf(cdr u,cdr v);

symbolic procedure gftimes(u,v);
   if atom car u then gfftimes(u,v) else gbftimes(u,v);

symbolic procedure gfftimes(u,v);
   begin scalar ru,iu,rv,iv;
         ru := car u; iu := cdr u; rv := car v; iv := cdr v;
         return (ru*rv - iu*iv) . (ru*iv+iu*rv) end;

symbolic procedure gbftimes(u,v);
   begin scalar ru,iu,rv,iv;
         ru := car u; iu := cdr u; rv := car v; iv := cdr v;
         return (difbf(timbf(ru,rv),timbf(iu,iv))) .
            plubf(timbf(ru,iv),timbf(iu,rv)) end;

symbolic procedure gfquotient(u,v);
   if atom car u then gffquot(u,v) else gbfquot(u,v);

symbolic procedure gffquot(u,v);
   begin scalar ru,iu,rv,iv,d;
         ru := car u; iu := cdr u; rv := car v; iv := cdr v;
         d := rv*rv+iv*iv;
         return ((ru*rv+iu*iv)/d) . ((iu*rv - ru*iv)/d) end;

symbolic procedure gbfquot(u,v);
   begin scalar ru,iu,rv,iv,d;
         ru := car u; iu := cdr u; rv := car v; iv := cdr v;
         d := plubf(timbf(rv,rv),timbf(iv,iv));
         return divbf(plubf(timbf(ru,rv),timbf(iu,iv)),d) .
            divbf(difbf(timbf(iu,rv),timbf(ru,iv)),d) end;

symbolic procedure gfminus u; (bfminus car u) . (bfminus cdr u);

symbolic procedure gfrotate u; (bfminus cdr u) . (car u);

%symbolic procedure gfrl u; car u;

%symbolic procedure gfim u; cdr u;

%symbolic procedure gfzerop u;
%   if not atom gfrl u then mt!: gfrl u = 0 and mt!: gfim u = 0
%      else equal(u,(0.0 . 0.0));

symbolic procedure gfdot(u,v);
   if atom car u then gffdot(u,v) else gbfdot(u,v);

symbolic procedure gffdot(u,v); car u*car v+cdr u*cdr v;

symbolic procedure gbfdot(u,v);
   plubf(timbf(car u,car v),timbf(cdr u,cdr v));

symbolic procedure gfrsq u; gfdot(u,u);

symbolic procedure gffrsq u; car u*car u+cdr u*cdr u;

symbolic procedure gbfrsq u;
   plubf(timbf(car u,car u),timbf(cdr u,cdr u));

symbolic procedure gffmult(r,u); (r*car u) . (r*cdr u);

symbolic procedure gffsqrt x;
   begin scalar x0,nx,xd,xd0,rl,im; rl := gfrl x; im := gfim x;
     rl := sqrt(hypot(rl,im)/2+rl/2); im := im/(2*rl); nx := rl . im;
     repeat
       <<x0 := nx;
         nx := gffmult(0.5,gffplus(x0,gffquot(x,x0)));
         xd0 := xd; xd := gffrsq gffdiff(x,gfftimes(nx,nx))>>
     until xd0 and xd0 - xd<=0.0; return x0 end;

symbolic procedure gbfmult(r,u);
   <<r := bfloat r; (timbf(r,car u)) . (timbf(r,cdr u))>>;

symbolic procedure gbfsqrt x;
   begin scalar x0,nx,xd,xd0,rl;
     nx :=
      <<rl := (bfsqrt timbf(bfhalf!*,plubf(bfsqrt gfrsq x,gfrl x)));
        rl . timbf(bfhalf!*,divbf(gfim x,rl))>>;
     repeat
       <<x0 := nx;
         nx := gbfmult(bfhalf!*,gbfplus(x0,gbfquot(x,x0)));
         xd0 := xd; xd := gbfrsq gbfdiff(x,gbftimes(nx,nx))>>
     until xd0 and mt!: difbf(xd0,xd)<=0; return x0 end;

symbolic smacro procedure rl2gfc x;
   x . if atom x then 0.0 else bfz!*;

symbolic procedure gfsqrt x;
 % computes gfsqrt x by Newton's method, for both gf and gbf.
   begin scalar xn,neg,negi;
         if gfzerop x then return x;
         if bfminusp gfrl x
           then <<x := gfminus x; neg := t;
                  if bfminusp gfim x then nil else negi := t >>;
         if bfzp gfim x then
            <<x := gfrl x;
              xn := rl2gfc(if atom x then sqrt x else bfsqrt x);
              go to ret>>;
         xn := if atom gfrl x then gffsqrt x else gbfsqrt x;
         if negi then xn := gfminus xn;
    ret: return if neg then gfrotate xn else xn end;

symbolic procedure sgn x; if x>0 then 1 else if x<0 then -1 else 0;

symbolic procedure exptbf(x,n,a);
  % Computes a*x**n in bfloat arithmetic for positive x
  % and positive integer n.
   begin
  lp: if oddintp n then a := timbf(a,x); % not evenp n
      n := lshift (n, -1);
      if n=0 then return a;
      x := timbf(x,x); go to lp end;

symbolic procedure icbrt(x);
   % x is a number : result is integer s approx cube root of x,
   % i.e. if x > 0 then s**3 <= n < (s+1)**3 o/w s**3 >= n > (s-1)**3.
   irootn(fix2(x),3);

symbolic procedure fix2 x; if fixp x then x else fl2int x;

symbolic procedure ilog2 n;
   % n is an integer.  Result is an integer r, s.t. 2**r <= abs(n)
   %  < 2**(r+1).
   begin scalar ans, powers!-of!-2, pwr;
      if n<=0 then terrlst(n,'ilog2);
      pwr := 2;
      powers!-of!-2 := pwr . nil;
      while n>pwr do
	 <<powers!-of!-2 := pwr . powers!-of!-2; pwr := pwr*pwr>>;
      ans := 0;
      while(pwr := car(powers!-of!-2)) neq 2 do
	<<powers!-of!-2 := cdr(powers!-of!-2);
	  if n >= pwr then << n := n/pwr; ans := ans + 1; >>;
	  ans := ans*2>>;
      if n >= 2 then ans := ans + 1;
      return ans
   end;

symbolic procedure isqrt(x);
   % x is a number : result is integer s approx square root of x,
   % i.e. if x > 0 then s**2 <= n < (s+1)**2 o/w s**2 >= n > (s-1)**2.
   if x<=0 then terrlst(x,'isqrt) else irootn(fix2(x), 2);

symbolic procedure qroundup(m,n);
   % m, n are integers, n>0 : result is least integer >= m/n.
   if m<0 then -((-m)/n) else (m+n-1)/n;

symbolic procedure irootn(n,r);
    % n, r are integers : result is integer s approx r'th root of n,
    % i.e. if n > 0 then s**r <= n < (s+1)**r o/w s**r >= n > (s-1)**r.
    if not fixp n then typerr(n,"integer")
     else if not fixp r or r<=0 then typerr(r,"positive integer")
     else if n<0
      then if evenp r then typerr(r,"odd integer") else -irootn(-n,r)
     else if r = 1 then n
     else if n = 0 then 0
     else begin scalar ans;
	     ans := irootn1(n,r,ilog2 n,(ilog2 r)/2);
	     if ans**r>n then return ans - 1 else return ans
       end;

symbolic procedure irootn!-power2(p, q);
   % p, q are positive integers.
   % Result is an integer (slightly) greater than 2**(p/q).
   % Uses the first few terms of the Taylor expansion (tweaked a bit).
   % Error (of numans/denans) is at most 0.03%.
   begin scalar whole!-part, p1, numans, denans;
      whole!-part := (p+q/2)/q;
      p1 := p-q*whole!-part;
      numans := q^3/100 + 1000*q^3 + 693*q^2*p1 + 243*q*p1^2 + 57*p1^3;
      denans := 1000*q^3;
      return 1+ (2^whole!-part * numans)/denans
      % add 1 to force rounding up
   end;

symbolic procedure irootn1(n,r,logn,xs);
   % n, r integers >0, logn is ilog2(n), xs is the excess # bits in the
   % top 1/2.
   % result is s, s.t. s**r <= n < (s+1)**r or (s-1)**r <= n < s**r.
   begin scalar x, upb, size, tmp;
      size := logn / r;
      if size < 17 then upb := irootn!-power2(1+logn, r)
       else <<x := size/2 - xs;
	      upb := irootn1(n/2**(x*r), r, logn-r*x, xs)*2**x;
	      tmp := upb**(r-1);
	      return (((r-1)*upb)*tmp+n)/(r*tmp)>>;
      repeat <<x := upb; upb := x - qroundup(x - n/x**(r-1),r)>>
	 until upb >= x;
      return x
   end;

put('irootn,'number!-of!-args,2);   %  For VALUECHK.

endmodule;

end;
REDUCE Historical
[ REDUCE Sourceforge Project | Historical SVN Repository | GitHub Mirror | SourceHut Mirror | NotABug Mirror | Chisel Mirror | Chisel RSS ]