Check-ins Using Download Hex

Artifact 00ce57f275aef13aa92841a680e10c595186db634af5d831583040f9723b7690:

module rdelem;  % Elementary functions in rounded domain.

exports deg2rad!*, quotient!:, rad2deg!*, rdacos!*, rdacosd!*,
        rdacosh!*, rdacot!*, rdacotd!*, rdacoth!*, rdacsc!*, rdacscd!*,
        rdacsch!*, rdarg!*, rdasec!*, rdasecd!*, rdasech!*, rdasin!*,
        rdasind!*, rdasinh!*, rdatan!*, rdatan2!*, rdatan2d!*,
        rdatand!*, rdatanh!*, rdcbrt!*, rdcos!*, rdcosd!*, rdcosh!*,
        rdcot!*, rdcotd!*, rdcoth!*, rdcsc!*, rdcscd!*, rdcsch!*,
        rde!*, rdexp!*, rdexpt!*, rdhalf!*, rdhypot!*, rdlog!*,
        rdlog10!*, rdlogb!*, rdnorm!*, rdone!*, rdpi!*, rdsec!*,
	rdsecd!*, rdsech!*, rdsin!*, rdsind!*, rdsinh!*,
        rdsqrt!*, rdtan!*, rdtand!*, rdtanh!*, rdtwo!*, rdzero!*,
        texpt!:, texpt!:any;

imports !*f2q, abs!:, acos, acos!*, acosd, acosh, acot, acotd, acoth,
        acsc, acscd, acsch, asec, asecd, asech, asin, asin!*, asind,
        asinh, atan, atan!*, atan2, atan2d, atand, atanh, bflerrmsg,
        bfloat, bfp!:, bfsqrt, cbrt, conv!:bf2i, conv!:bf2i, conv!:mt,
        convprec, cos, cos!*, cosd, cosh, cot, cotd, coth, csc, cscd,
        csch, difbf, divbf, e!*, ep!:, eqcar, equal!:, exp, exp!*,
        exp!:, exptbf, geq, greaterp!:, hypot, i2rd!*, incprec!:,
        invbf, leq, leq!:, lessp!:, log, log!*, log10, log!:,
        logb, logfp, lshift, make!:ibf, minus!:, minusp!:, mk!*sq,
        mkround, mt!:, neq, pi!*, plubf, preci!:, rd!:minus,
        rd!:minusp, read!:num, rndbfon, round!*, round!:last,
        round!:mt, sec, secd, sech, sgn, simprd, sin, sin!*, sind,
        sinh, sqrt, sqrt!:, tan, tan!*, tand, tanh, terrlst, timbf,
        times!:;

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

global '(bfz!* bfone!* bften!* bfhalf!* !:180!* !:bf1!.5!* bfthree!*
 !:bf60!* epsqrt!* bftwo!* !!pii !!flprec !!rdprec !!shbinfl
 pi!/180 !180!/pi !!ee !!maxarg);

pi!/180 := !!pii/180;  !180!/pi := 180/!!pii;

fluid '(!*numval);

deflist('((exp rdexp!*) (expt rdexpt!*) (log rdlog!*) (sin rdsin!*)
   (cos rdcos!*) (tan rdtan!*) (asin rdasin!*) (acos rdacos!*)
   (atan rdatan!*) (sqrt rdsqrt!*) (sinh rdsinh!*) (cosh rdcosh!*)
   (sec rdsec!*) (csc rdcsc!*) (cot rdcot!*) (tanh rdtanh!*)
   (coth rdcoth!*) (sech rdsech!*) (csch rdcsch!*) (asinh rdasinh!*)
   (acosh rdacosh!*) (acot rdacot!*) (asec rdasec!*) (acsc rdacsc!*)
   (atanh rdatanh!*) (acoth rdacoth!*) (asech rdasech!*)
   (acsch rdacsch!*) (logb rdlogb!*) (log10 rdlog10!*) (ln rdlog!*)
   (atan2 rdatan2!*) (hypot rdhypot!*) % (cbrt rdcbrt!*)
   (deg2rad deg2rad!*) (rad2deg rad2deg!*) (deg2dms deg2dms!*)
   (rad2dms rad2dms!*) (dms2deg dms2deg!*) (dms2rad dms2rad!*)
   (norm rdnorm!*) (arg rdarg!*) (e rde!*) (pi rdpi!*)),
   '!:rd!:);

% deflist('((sind rdsind!*) (cosd rdcosd!*) (asind rdasind!*) (acosd
%    rdacosd!*) (tand rdtand!*) (cotd rdcotd!*) (atand rdatand!*) (acotd
%    rdacotd!*) (secd rdsecd!*) (cscd rdcscd!*) (asecd rdasecd!*) (acscd
%    rdacscd!*) (atan2d rdatan2d!*)),'!:rd!:);



for each n in '(exp sin cos tan asin acos atan sinh cosh  % log
    sec csc cot tanh coth sech csch asinh acosh acot asec acsc atanh
    acoth asech acsch logb log10 ln atan2 hypot
%   sind cosd asind acosd tand cotd atand acotd secd cscd asecd acscd
%   atan2d cbrt
    deg2rad rad2deg deg2dms rad2dms dms2deg dms2rad norm arg argd)
       do put(n,'simpfn,'simpiden);

flag('(dms2deg dms2rad),'listargp);

deflist('((dms2deg!* simpdms) (dms2rad!* simpdms)), 'simparg);

deflist('((atan2 2) (hypot 2) (atan2d 2) (logb 2)),
  'number!-of!-args);

flag('(acsc sind asind tand atand cotd acotd cscd acscd csch
       acsch deg2rad rad2deg),'odd);   % sgn.

flag('(cosd secd),'even);

flag('(cotd sech),'nonzero);

symbolic procedure rdexp!* u; mkround
  (if not atom x then exp!* x
   else if x>!!maxarg then <<rndbfon(); exp!* bfloat x>>
   else if x<-!!maxarg then 0.0 else exp x)
   where x=convprec u;

symbolic procedure rdsqrt!* u;
   mkround(if atom x then sqrt x else bfsqrt x)
   where x=convprec u;

symbolic procedure rdexpt!*(u,v);
   mkround
     (if not atom x then texpt!:any(x,y) else
       if zerop x then if zerop y then rederr "0**0 formed" else u else
      ((if z>!!maxarg then <<rndbfon(); texpt!:any(bfloat x,bfloat y)>>
        else if z<-!!maxarg then 0.0 else rexpt(x,y))
        where z=y*logfp bfloat abs x))
      where x=convprec u,y=convprec v;

symbolic procedure rdlog!* u;
   mkround(if atom x then log x else log!* x)
   where x=convprec u;

% symbolic procedure rdsgn!* u;
%   (if atom x then sgn x else sgn mt!: x) where x=round!* u;

symbolic procedure rdatan2!*(u,v);
   if !:zerop u and !:zerop v
     then rerror(arith,8,"0/0 formed")
    else (mkround(if atom x then atan2(x,y) else atan2!*(x,y))
	  where x=convprec u,y=convprec v);

% symbolic procedure rdatan2d!*(u,v);
%    mkround(if atom x then atan2d(x,y) else rad2deg!: atan2!*(x,y))
%    where x=convprec u,y=convprec v;

symbolic procedure atan2!*(y,x);
   if mt!: x=0 then if (y := mt!: y)=0 then bfz!* else
      <<x := pi!/2!*(); if y<0 then minus!: x else x>>
   else <<(if mt!: x>0 then a
      else if mt!: y<0 then difbf(a,pi!*())
         else plubf(a,pi!*()))
     where a=atan!* divbf(y,x)>>;

% symbolic procedure atan2d!*(y,x);
%    if mt!: x=0 then if (y := mt!: y)=0 then bfz!* else
%       <<x := timbf(!:180!*,bfhalf!*); if y<0 then minus!: x else x>>
%    else <<(if mt!: x>0 then a
%       else if mt!: y<0 then difbf(a,!:180!*) else plubf(a,!:180!*))
%      where a=rad2deg!: atan!* divbf(y,x)>>;

symbolic procedure rde!*; mkround if !*!*roundbf then e!*() else !!ee;

symbolic procedure rdpi!*;
   mkround if !*!*roundbf then pi!*() else !!pii;

symbolic procedure pi!/2!*; timbf(bfhalf!*,pi!*());

symbolic procedure deg2rad!* u;
   mkround(if atom x then deg2rad x else deg2rad!: x)
   where x=convprec u;

symbolic procedure rad2deg!* u;
   mkround(if atom x then rad2deg x else rad2deg!: x)
   where x=convprec u;

symbolic procedure deg2rad x; x*pi!/180;

symbolic procedure rad2deg x; x*!180!/pi;

symbolic procedure deg2rad!: x; divbf(timbf(x,pi!*()),!:180!*);

symbolic procedure rad2deg!: x; divbf(timbf(x,!:180!*),pi!*());

symbolic procedure rdsin!* u;
   mkround (if atom x then sin x else sin!* x)
   where x=convprec u;

% symbolic procedure rdsind!* u;
%    mkround (if atom x then sind x else sin!* deg2rad!: x)
%    where x=convprec u;

symbolic procedure rdcos!* u;
   mkround(if atom x then cos x else cos!* x)
   where x=convprec u;

% symbolic procedure rdcosd!* u;
%    mkround(if atom x then cosd x else cos!* deg2rad!: x)
%   where x=convprec u;

symbolic procedure rdtan!* u;
   mkround(if atom x then tan x else tan!* x)
   where x=convprec u;

% symbolic procedure rdtand!* u;
%    mkround(if atom x then tand x else tan!* deg2rad!: x)
%   where x=convprec u;

symbolic procedure rdasin!* u;
   mkround(if atom x then asin x else asin!* x)
   where x=convprec u;

% symbolic procedure rdasind!* u;
%    mkround(if atom x then asind x else rad2deg!: asin!* x)
%    where x=convprec u;

symbolic procedure rdacos!* u;
   mkround(if atom x then acos x else acos!* x)
   where x=convprec u;

% symbolic procedure rdacosd!* u;
%    mkround(if atom x then acosd x else rad2deg!: acos!* x)
%    where x=convprec u;

symbolic procedure rdatan!* u;
   mkround(if atom x then atan x else atan!* x)
   where x=convprec u;

% symbolic procedure rdatand!* u;
%    mkround(if atom x then atand x else rad2deg!: atan!* x)
%   where x=convprec u;

symbolic procedure rdsinh!* u;
   mkround(if atom x then sinh x else sinh!* x)
   where x=convprec u;

symbolic procedure rdcosh!* u;
   mkround(if atom x then cosh x else cosh!* x)
   where x=convprec u;

% these redefine functions that are in bfelem, and are faster.

symbolic procedure sinh!* x;
   timbf(bfhalf!*,difbf(y,invbf y)) where y=exp!* x;

symbolic procedure cosh!* x;
   timbf(bfhalf!*,plubf(y,invbf y)) where y=exp!* x;


% no bfelem functions after this point.

symbolic procedure rdsec!* u;
   mkround(if atom x then sec x else invbf cos!* x)
   where x=convprec u;

% symbolic procedure rdsecd!* u;
%    mkround(if atom x then secd x else invbf cos!* deg2rad!: x)
%    where x=convprec u;

symbolic procedure rdcsc!* u;
   mkround(if atom x then csc x else invbf sin!* x)
   where x=convprec u;

% symbolic procedure rdcscd!* u;
%    mkround(if atom x then cscd x else invbf sin!* deg2rad!: x)
%   where x=convprec u;

symbolic procedure rdcot!* u;
   mkround(if atom x then cot x else tan!* difbf(pi!/2!*(),x))
   where x=convprec u;

% symbolic procedure rdcotd!* u;
%   mkround(if atom x then cotd x else tan!* difbf(pi!/2!*(),
%           deg2rad!: x))
%    where x=convprec u;

symbolic procedure rdtanh!* u;
   mkround(if atom x then tanh x else divbf(sinh!* x,cosh!* x))
   where x=convprec u;

symbolic procedure rdcoth!* u;
   mkround(if atom x then coth x else divbf(cosh!* x,sinh!* x))
   where x=convprec u;

symbolic procedure rdsech!* u;
   mkround(if atom x then sech x else invbf cosh!* x)
   where x=convprec u;

symbolic procedure rdcsch!* u;
   mkround(if atom x then csch x else invbf sinh!* x)
   where x=convprec u;

symbolic procedure rdasinh!* u;
   mkround(if atom x then asinh x else asinh!* x)
   where x=convprec u;

symbolic procedure rdacosh!* u;
   mkround(if atom x then acosh x else acosh!* x)
   where x=convprec u;

symbolic procedure asinh!* x; begin scalar s;
   if minusp!: x then x := minus!: x else s := t;
   x := if leq!:(x,epsqrt!*) then x
      else log!* plubf(x,
         if lessp!:(x,bftwo!*) then bfsqrt plubf(timbf(x,x),bfone!*)
         else if lessp!:(invbf x,epsqrt!*) then x
         else timbf(x,bfsqrt plubf(bfone!*,divbf(bfone!*,timbf(x,x)))));
   return if s then x else minus!: x end;

symbolic procedure acosh!* x;
   if lessp!:(x,bfone!*) then terrlst(x,'acosh)
   else log!* plubf(x,if leq!:(invbf x,epsqrt!*) then x
      else timbf(x,bfsqrt difbf(bfone!*,divbf(bfone!*,timbf(x,x)))));

symbolic procedure rdacot!* u;
   mkround(if atom x then acot x
      else difbf(pi!/2!*(),atan!* x))
   where x=convprec u;

% symbolic procedure rdacotd!* u;
%   mkround(if atom x then acotd x
%       else rad2deg!: difbf(pi!/2!*(),atan!* x))
%    where x=convprec u;

symbolic procedure rdasec!* u;  % not yet
   mkround(if atom x then asec x else
      difbf(pi!/2!*(),asin!* invbf x))
   where x=convprec u;

% symbolic procedure rdasecd!* u;  % not yet
%    mkround(if atom x then asecd x else
%       rad2deg!: difbf(pi!/2!*(),asin!* invbf x))
%    where x=convprec u;

symbolic procedure rdacsc!* u;
   mkround(if atom x then acsc x else asin!* invbf x)
   where x=convprec u;

% symbolic procedure rdacscd!* u;
%   mkround(if atom x then acscd x else rad2deg!: asin!* invbf x)
%   where x=convprec u;

symbolic procedure rdatanh!* u;
   mkround(if atom x then atanh x else atanh!* x)
   where x=convprec u;

symbolic procedure atanh!* x;
   if not greaterp!:(bfone!*,abs!: x) then terrlst(x,'atanh)
   else if leq!:(abs!: x,epsqrt!*) then x
   else timbf(bfhalf!*,
      log!* divbf(plubf(bfone!*,x),difbf(bfone!*,x)));

symbolic procedure rdacoth!* u;
   mkround(if atom x then acoth x else atanh!* invbf x)
   where x=convprec u;

symbolic procedure rdasech!* u;   % not from here down
   mkround(if atom x then asech x
      else if leq!:(x,bfz!*) or greaterp!:(x,bfone!*)
         then terrlst(x,'asech) else acosh!* invbf x)
   where x=convprec u;

symbolic procedure rdacsch!* u;
   mkround(if atom x then acsch x
      else if mt!: x=0 then terrlst(x,'acsh) else asinh!* invbf x)
   where x=convprec u;

symbolic procedure rdlogb!*(u,v);
   mkround(if atom x then logb(x,b) else logb!*(x,b))
   where x=convprec u,b=convprec v;

symbolic procedure rdlog10!* u;
   mkround(if atom x then log10 x else logb!*(x,bften!*))
   where x=convprec u;

symbolic procedure logb!* (x,b); %log x to base b;
   begin scalar a,s;
      a := greaterp!:(x,bfz!*);
      s := not(leq!:(b,bfz!*) or equal!:(b,bfone!*));
      if a and s then return divbf(log!* x,log!* b)
         else terrlst((if a then list ('base,b)
            else if s then list('arg,x) else list(x,b)),'logb) end;

% symbolic procedure rdcbrt!* u;
%    mkround(if atom x then cbrt x else cbrt!* x)
%    where x=convprec u;

% symbolic procedure cbrt!* x;
%    begin scalar s,l,g,u,nx,r; u := bfone!*;
%          if mt!: x=0 or equal!:(abs!: x,u) then return x
%          else if minusp!: x then x := minus!: x else s := t;
%          if lessp!:(x,u) then <<x := divbf(u,x); l := t>>
%             else if equal!:(x,u) then go to ret;
%          nx := '!:bf!: .
%             <<r := remainder(ep!:(nx := conv!:mt(x,3)),3);
%               if r=0 then (5+mt!: nx/179) . (ep!: nx/3)
%               else if r=1 or r=-2
%                then (10+mt!: nx/74) . ((ep!: nx-1)/3)
%               else (22+mt!: nx/35) . ((ep!: nx-2)/3)>>;
%    loop: r := nx;
%          nx := plubf(divbf(r,!:bf1!.5!*),
%             divbf(x,timbf(r,timbf(r,bfthree!*))));
%          if g and leq!:(r,nx) then go to ret;
%          g := t; go to loop;
%     ret: if l then nx := divbf(u,nx);
%          return if s then nx else minus!: nx end;

symbolic procedure rdhypot!*(u,v);
   mkround(if atom p then hypot(p,q) else hypot!*(p,q))
   where p=convprec u,q=convprec v;

symbolic procedure hypot!*(p,q);
 % Hypot(p,q)=sqrt(p*p+q*q) but avoids intermediate swell.
 begin scalar r;
   if minusp!: p then p := minus!: p; if minusp!: q then q := minus!: q;
   if mt!: p=0 then return q else if mt!: q=0 then return p
   else if lessp!:(p,q) then <<r := p; p := q; q := r>>;
   r := divbf(q,p);
   return if lessp!:(r,epsqrt!*) then p
      else timbf(p,bfsqrt plubf(bfone!*,timbf(r,r))) end;

symbolic procedure simpdms l;
   % Converts argument of form ({d,m,s}) to rd ((d m s)) if possible.
   if cdr l or atom (l := car l) or not eqcar(l,'list)
      or length l neq 4 then nil else
   begin scalar fl;
      l := for each a in cdr l collect
          if not (null(a := simprd list a) and (fl := t))
             then a := car a;
      if not fl then return list l end;

symbolic procedure round2a!* a; if atom a then a else round!* a;

symbolic procedure dms2rad!* u; deg2rad!* dms2deg!* u;

symbolic procedure dms2deg!* u;
   mkround(if atom caddr l then dms2deg l else dms2deg!: l)
   where l=list(round2a!* car u,round2a!* cadr u,round!* caddr u);

symbolic procedure dms2deg l; ((caddr l/60.0+cadr l)/60.0+car l);

symbolic procedure dms2deg!: l;
   plubf(bfloat car l,divbf(plubf(bfloat cadr l,
      divbf(bfloat caddr l,!:bf60!*)),!:bf60!*));

symbolic procedure rad2dms x; deg2dms rad2deg x;

symbolic procedure rad2dms!* u; deg2dms!* rad2deg!* u;

symbolic procedure deg2dms!* u;
   mklist3!*(if atom x then deg2dms x else deg2dms!: x)
   where x=round2a!* u;

symbolic procedure mklist3!* x; % floats seconds if not integer.
   'list . list(car x,cadr x,
      <<(if atom s and zerop(s-fix s) then fix s
         else if not atom s and integerp!: s then conv!:bf2i s
         else mk!*sq !*f2q mkround s) where s=caddr x>>);

symbolic procedure deg2dms x; % dms output in form list(d,m,s);
   begin integer d,m;
%     m := fix(x := 60.0*(x-(d := fix2 x)));
      m := fix(x := 60.0*(x-(d := fix x)));
      return list(d,m,60.0*(x-m)) end;

symbolic procedure deg2dms!: x; % dms output in form list(d,m,s).
   begin integer d,m;
      d := conv!:bf2i x;
      m := conv!:bf2i(x := timbf(!:bf60!*,difbf(x,bfloat d)));
      return list(d,m,timbf(!:bf60!*,difbf(x,bfloat m))) end;

symbolic procedure rdnorm!* u; if rd!:minusp u then rd!:minus u else u;

symbolic procedure rdarg!* u;
   if rd!:minusp u then rdpi!*() else rdzero!*();

% the following bfloat definitions are needed in addition to files
% smbflot and bfelem.red to support rdelem.

global '(!:bfone!* bftwo!* bfhalf!* bfz!* !:bf!-0!.25);

symbolic procedure rdone!*; if !*!*roundbf then bfone!* else 1.0;

symbolic procedure rdtwo!*; if !*!*roundbf then bftwo!* else 2.0;

symbolic procedure rdhalf!*; if !*!*roundbf then bfhalf!* else 0.5;

symbolic procedure rdzero!*; if !*!*roundbf then bfz!* else 0.0;

symbolic procedure texpt!:(nmbr, k);
% This function calculates the Kth power of "n" up to the
% binary precision specified by !:BPREC!:. %SK
% NMBR is a BINARY BIG-FLOAT representation of "n" and K an integer.
   if not fixp k then bflerrmsg 'texpt!:  % use texpt!:any in this case.
    else if k=0 then bfone!*
    else if k=1 then nmbr
    else if k<0 then invbf texpt!:(nmbr,-k) %SK
    else exptbf(nmbr,k,bfone!*); %SK

symbolic procedure quotient!:(n1, n2);
% This function calculates the integer quotient of "n1"
%      and "n2", just as the "QUOTIENT" for integers does.
% **** For calculating the quotient up to a necessary
% ****      precision, please use DIVIDE!:.
% N1 and N2 are BIG-FLOAT representations of "n1" and "n2".
begin integer e1, e2;
  if (e1 := ep!: n1) = (e2 := ep!: n2) then return
             make!:ibf(mt!: n1 / mt!: n2, 0)
   else if e1 > e2 then return
             quotient!:(incprec!:(n1, e1 - e2) , n2)
   else return
             quotient!:(n1, incprec!:(n2, e2 - e1));
end$

symbolic procedure texpt!:any(x, y);
  %modified by SK to use bfsqrt and exp!*, invbf and timbf.
% This function calculates the power x**y, where "x"
%      and "y" are any numbers.  The precision of
%      the result is specified by !:PREC!:. % SK
% **** For a negative "x", this function returns
% ****      -(-x)**y unless "y" is an integer.
% X is a BIG-FLOAT representation of "x".
% Y is either an integer, a floating-point number,
%      or a BIG-FLOAT number, i.e., a BIG-FLOAT
%      representation of "y".
    if equal!:(x,e!*()) then exp!* bfloat y
    else if fixp y then texpt!:(x, y)
    else if integerp!: y then texpt!:(x,conv!:bf2i y)
    else if not(bfp!: y or bfp!:(y := read!:num y))
     then bflerrmsg 'texpt!:any     % read!:num probably not necessary.
    else if minusp!: y then invbf texpt!:any(x,minus!: y) %SK
    else if equal!:(y,bfhalf!*) then bfsqrt x   %SK
    else if equal!:(y,!:bf!-0!.25) then bfsqrt bfsqrt x   %SK
    else begin integer n;  scalar xp, yp;
          n := (if !:bprec!: then !:bprec!:
                else max(preci!: x, preci!: y));
%          if minusp!: x then xp:=minus!: x else xp := x;
          if minusp!: x then bflerrmsg 'texpt!:any
           else xp := x;
          if integerp!: times!:(y,bftwo!*) then
             << xp := incprec!:(xp, 1);
                yp := texpt!:(xp, conv!:bf2i y);
                yp := times!:(yp, sqrt!:(xp, n + 1));
                yp := round!:mt(yp, n) >>
          else
             << yp := timbf(y, log!:(xp, n + 1)); %SK
                yp := exp!:(yp, n) >>;
          return (if minusp!: x then minus!: yp else yp);
     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 remainder(mt!: x,lshift(1,-ep!: x)) = 0);

endmodule;

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