OS/2 Shareware BBS: 10 Tools/10-Tools.zip/sa104os2.zip/SATHR104.ZIP/SATHER/LIBRARY/FLT.SA

Text File | 1995-02-05 | 19KB | 597 lines

-- Copyright (C) International Computer Science Institute, 1994.  COPYRIGHT  --
-- NOTICE: This code is provided "AS IS" WITHOUT ANY WARRANTY and is subject --
-- to the terms of the SATHER LIBRARY GENERAL PUBLIC LICENSE contained in    --
-- the file "Doc/License" of the Sather distribution.  The license is also   --
-- available from ICSI, 1947 Center St., Suite 600, Berkeley CA 94704, USA.  --
--------> Please email comments to "[e-mail]". <----------

value class FLT < $IS_EQ{FLT}, $IS_LT{FLT}, $NIL{FLT} is
   -- IEEE 754-1984 "single" format 32-bit floating point.
   -- Most of these functions presently just call the FLTD versions
   -- because of C double/float weirdness.

   create(f:FLT):SAME is return f; end;

   plus(f:SAME):SAME is -- The sum of self and `f'. Built-in.
      raise "FLT::plus(SAME):SAME undefined.";
   end;

   minus(f:SAME):SAME is -- The difference between self and `f'. Built-in.
      raise "FLT::minus(SAME):SAME undefined.";
   end;

   negate:SAME is -- The negation of self. Same as zero minus self,
      -- except that IEEE does funny things for the sign bit and
      -- different rounding modes.
      raise "FLT::negate:SAME undefined.";
   end;

   times(f:SAME):SAME is -- The signed product of self and `f'. Built-in.
      raise "FLT::times(SAME):SAME undefined.";
   end;

   div(f:SAME):SAME is -- The quotient of self and `f'. Built-in.
      raise "FLT::div(SAME):SAME undefined.";
   end;

   is_eq(f:SAME):BOOL is -- True if self and `f' represent the same value.  Built-in.
      raise "FLT::is_eq(f:SAME):BOOL undefined.";
   end;

   is_neq(f:SAME):BOOL is
      -- True if self and `f' represent different values.  Built-in.
      -- See FLTD::is_neq for more about IEEE awkwardness.
      raise "FLT::is_neq(f:SAME):BOOL undefined.";
   end;

   is_lt(f:SAME):BOOL is -- True if self is less than `f'. Built-in.
      raise "FLT::is_lt(SAME):BOOL undefined.";
   end;

   is_leq(f:SAME):BOOL is -- True if self is less than or equal to `f'.  Built-in.
      raise "FLT::is_leq(SAME):BOOL undefined.";
   end;

   is_gt(f:SAME):BOOL is -- True if self is greater than `f' as signed integers. Built-in.
      raise "FLT::is_gt(SAME):BOOL undefined.";
   end;

   is_geq(f:SAME):BOOL is -- True if self is greater than or equal to `f'.  Built-in.
      raise "FLT::is_geq(SAME):BOOL undefined.";
   end;

   is_within(tolerance,val:SAME):BOOL is
       return (self-val).abs<=tolerance;
   end;

   -- IEEE functions.

   is_finite:BOOL is -- returns true if zero, subnormal or normal.
       return C_FLT::ir_finite_(self);
   end;

   is_inf:BOOL is -- returns true if infinite
       return C_FLT::ir_isinf_(self);
   end;

   is_nan:BOOL is -- returns true if NaN
       return C_FLT::ir_isnan_(self);
   end;

   is_normal:BOOL is -- returns true if normal
       return C_FLT::ir_isnormal_(self);
   end;

   is_subnormal:BOOL is -- returns true if subnormal
       return C_FLT::ir_issubnormal_(self);
   end;

   is_zero:BOOL is -- returns true is zero
       return C_FLT::ir_iszero_(self);
   end;

   signbit_set:BOOL is -- returns true if sign bit of self is set
       return C_FLT::ir_signbit_(self);
   end;

   unbiased_exponent:INT is
       -- return unbiased exponent of self as an INT;
       -- for zero this is INT::maxint.negate, for an
       -- infinite it is INT::maxint.  If subnormal,
       -- normalization occurs first.
       return C_FLT::ir_ilogb_(self);
   end;

   copysign(y:SAME):SAME is
       -- return self with the sign bit set to the same as y's sign bit.
       return C_FLT::r_copysign_(self,y);
   end;

   nextup:SAME is -- return next representable number from self.
       return C_FLT::r_nextafter_(self,1.flt);
   end;

   nextdown:SAME is -- return previous representable number from self.
       return C_FLT::r_nextafter_(self,-1.flt);
   end;

   --  x.remainder(y) and x.mod(y) return a remainder of x with respect
   --  to y; that is, the result r is one of the numbers that differ from
   --  x by an integral multiple of y.  Thus (x-r)/y  is an integral
   --  value, even though it might exceed INT::maxint if it were
   --  explicitly computed as an INT.  Both functions return one  of the
   --  two such r smallest in magnitude.  remainder(x,y) is the operation
   --  specified in ANSI/IEEE Std 754-1985; the result of x.mod(y) may
   --  differ from remainder's result by +-y.  The magnitude of
   --  remainder's result can not exceed half that of y; its sign might
   --  not agree with either x or y.  The magnitude of mod's result is
   --  less than that of y; its sign agrees with that of x.  Neither
   --  function will raise an exception as long as both arguments are
   --  normal or subnormal.  x.remainder(0), x.mod(0), infinity.remainder(y),
   --  and infinity.mod(y) are invalid operations that produce a NaN.

   remainder(y:SAME):SAME is
      return C_FLT::r_remainder_(self,y);
   end;

   mod(y:SAME):SAME is
      return C_FLT::r_fmod_(self,y);
   end;

   scale_by(n:INT):SAME is
      -- return x*2.pow(n) computed by exponent manipulation rather
      -- than by actually performing an exponentiation or a multiplication.
      -- 1 <= x.abs.scale_by(-x.unbiased_exponent) < 2 for every x
      -- except 0, infinity, and NaN.
      return C_FLT::r_scalbn_(self,n);
   end;

   -- Bessel functions of the first and second kinds.  y0, y1 and yn have
   -- logarithmic singularities at the origin, so they treat zero and
   -- negative arguments the way log does.

   bessel_j0:SAME is
      return self.fltd.bessel_j0.flt;
   end;

   bessel_j1:SAME is
      return self.fltd.bessel_j1.flt;
   end;

   bessel_jn(n:INT):SAME is
      return self.fltd.bessel_jn(n).flt;
   end;

   bessel_y0:SAME is
      return self.fltd.bessel_y0.flt;
   end;

   bessel_y1:SAME is
      return self.fltd.bessel_y1.flt;
   end;

   bessel_yn(n:INT):SAME is
      return self.fltd.bessel_yn(n).flt;
   end;

   -- Error functions

   erf:SAME is
      -- error function x.erf = (1/sqrt(pi))*integrate(0,x,exp(-t^2)dt)
      return self.fltd.erf.flt;
   end;

   one_minus_erf:SAME is
      -- 1.0-self.erf, but computed in a way to avoid cancellation for large self.
      return self.fltd.one_minus_erf.flt;
   end;

   -- Exponential, logarithm, power functions.  All these functions handle
   -- exceptional arguments in the spirit of IEEE 754-1985.  So:
   -- 0.log is -infinity with a division by zero exception
   -- For x<0, including -infinity, x.log is a quiet NaN with an invalid op exception
   -- For x=+infinity or a quiet NaN, x.log is x without exception
   -- For a signaling NaN, x.log is a quiet NaN with an invalid op exception
   -- 1.log is zero without exception
   -- For any other positive x, x.log is a normalized number with an inexact exception

   exp:SAME is -- The exponential e^self.
      return self.fltd.exp.flt;
   end;

   exp_minus_one:SAME is -- e^self-1.0, accurate even for tiny self.
      return self.fltd.exp_minus_one.flt;
   end;

   exp2:SAME is -- 2^self
      return self.fltd.exp2.flt;
   end;

   exp10:SAME is -- 10^self
      return self.fltd.exp10.flt;
   end;

   log:SAME is -- The natural logarithm of self.
      return self.fltd.log.flt;
   end;

   plus_one_log:SAME is -- (self+1).log, accurate even for tiny self.
      return self.fltd.plus_one_log.flt;
   end;

   log2:SAME is -- The logarithm base two of self.
      return self.fltd.log2.flt;
   end;

   log10:SAME is -- The logarithm base ten of self.
      return self.fltd.log10.flt;
   end;

   pow(arg:SAME):SAME is
      -- self raised to the arg'th power.  x.pow(0.0)=1.0 for all x.
      return self.fltd.pow(arg.fltd).flt;
   end;

   -- Hyperbolic functions.  They handle exceptional arguments in the
   -- spirit of IEEE 754-1985.  So:
   -- sinh and cosh return +-infinity on overflow
   -- acosh returns a NaN if its argument is less than 1.0
   -- atanh returns a NaN if its argument has an absolute value >1.0

   sinh:SAME is -- The hyperbolic sine of self.
       return self.fltd.sinh.flt;
   end;

   cosh:SAME is -- The hyperbolic cosine of self.
       return self.fltd.cosh.flt;
   end;

   tanh:SAME is -- The hyperbolic tangent of self.
       return self.fltd.tanh.flt;
   end;

   asinh:SAME is -- The inverse hyperbolic sine of self.
       return self.fltd.asinh.flt;
   end;

   acosh:SAME is -- The inverse hyperbolic cosine of self.
       return self.fltd.acosh.flt;
   end;

   atanh:SAME is -- The inverse hyperbolic tangent of self.
       return self.fltd.asinh.flt;
   end;

   -- Trigonometric functions.  These functions handle exceptional arguments
   -- in the spirit of IEEE 754-1985.  So:
   -- +-infinity.sin, +-infinity.cos, +-infinity.tan return NaN
   -- x.asin and x.acos with x.abs>1 return NaN
   -- sinpi etc. are similar except they compute self.sinpi=(self*pi).sin avoiding
   -- range-reduction issues because their definition permits range reduction
   -- that is fast and exact for all self.  The corresponding inverse functions
   -- compute asinpi(x)

   hypot(arg:SAME):SAME is
       -- sqrt(self*self+arg*arg), taking precautions against unwarranted
       -- IEEE exceptions.  +-infinity.hypot(arg) is +infinity for any arg,
       -- even a NaN, and is exceptional only for a signaling NaN.
       return self.fltd.hypot(arg.fltd).flt;
   end;

   sin:SAME is
       return self.fltd.sin.flt;
   end;

   cos:SAME is
       return self.fltd.cos.flt;
   end;

   -- This is commented out because I haven't implemented compiler-wise
   -- TUP classes yet.
   --sincos:TUP{SAME,SAME} is
   --    --Simultaneous computation of self.sin and self.cos.  This is faster
   --    --than independently computing them.
   --    return(C_FLT::c_flt_sincos(self));
   --end;

   tan:SAME is
       return self.fltd.tan.flt;
   end;

   asin:SAME is -- The arc sine of self in the range [-pi/2, pi/2]
       return self.fltd.asin.flt;
   end;

   acos:SAME is -- The arc sine of self in the range [0.0, pi]
       return self.fltd.acos.flt;
   end;

   atan:SAME is -- The arc tangent of self in the range [-pi/2, pi/2].
       return self.fltd.atan.flt;
   end;

   atan2(f:SAME):SAME is
      -- The arc tangent of self divided by arg in the range [-pi, pi].
      -- It chooses the quadrant specified by (self, arg).
      return self.fltd.atan2(f.fltd).flt;
   end;

   sinpi:SAME is
       return self.fltd.sinpi.flt;
   end;

   cospi:SAME is
       return self.fltd.cospi.flt;
   end;

   -- This is commented out because I haven't implemented compiler-wise
   -- TUP classes yet.
   --sincospi:TUP{SAME,SAME} is
   --  -- Simultaneous computation of self.sinpi and self.cospi.  This is faster
   --  -- than independently computing them.
   --    return(C_FLT::c_flt_sincospi(self));
   --end;

   tanpi:SAME is
       return self.fltd.tanpi.flt;
   end;

   asinpi:SAME is
       -- (1/pi) times the arc sine of self.  Result is in the range [-1/2,1/2]
       return self.fltd.asinpi.flt;
   end;

   acospi:SAME is -- The arc sine of self*pi in the range [0.0, pi]
       -- (1/pi) times the arc cosine of self.  Result is in the range [0,1]
       return self.fltd.acospi.flt;
   end;

   atanpi:SAME is
       -- (1/pi) times the arc tangent of self.
       -- Result is in the range [-1/2,1/2]
       return self.fltd.atanpi.flt;
   end;

   atan2pi(f:SAME):SAME is
      -- (1/pi) times the arc tangent of self divided by arg.
      -- Result in the range [-1, 1].
      -- It chooses the quadrant specified by (self, arg).
      return self.fltd.atan2pi(f.fltd).flt;
   end;

   -- Miscellaneous

   abs:SAME is -- The absolute value of self.
       return self.fltd.abs.flt;
   end;

   signum:SAME is
      if self<0.0 then return -1.0;
      elsif self>0.0 then return 1.0;
--    else return self;                                                         -- NLP
      end; return self;                                                         -- NLP
--    end;                                                                      -- NLP
   end;

   log_gamma:SAME is
      -- log gamma function.  x.ln_gamma=x.gamma.abs.log
      return self.fltd.log_gamma.flt;
   end;

   gamma:SAME is
      -- gamma function.
      return self.fltd.gamma.flt;
   end;

   sqrt:SAME is -- The square root of self.
      return self.fltd.sqrt.flt;
   end;

   square:SAME is -- The square of self.
      return self*self;
   end;

   cube_root:SAME is -- The square root of self.
      return self.fltd.cube_root.flt;
   end;

   cube:SAME is -- The cube of self.
      return self*self*self;
   end;

   max(arg:SAME):SAME is -- The larger of self and arg.
      -- Doesn't work properly is an argument is NaN.
--    if self<arg then return arg; else return self; end;                       -- NLP
      if self<arg then return arg; end; return self;                            -- NLP
   end;

   min(arg:SAME):SAME is -- The smaller of self and arg.
      -- Doesn't work properly is an argument is NaN.
--    if self<arg then return self; else return arg; end;                       -- NLP
      if self<arg then return self; end; return arg;                            -- NLP
   end;

   -- Conversions.
   private shared fbuf: FSTR;

   str:STR is
      -- A string version of self.
      if ((void(fbuf)) or (fbuf.size < 30)) then fbuf := #FSTR(30) end;
      fstr ::= str_in(fbuf);
      return(fstr.str); end;

   str(prec:INT):STR is
      -- A string version of self with arg digits of precision.
      des_sz ::= prec+10;
      if ((void(fbuf)) or (fbuf.size < des_sz)) then fbuf:=#FSTR(des_sz) end;
      fstr ::= str_in(fbuf,prec);
      return(fstr.str);  end;

   str_in(arg:FSTR):FSTR is
      store_in: FSTR;
      if (arg.size >= 30) then store_in := arg;
      else store_in := #FSTR(30) end;
      sz ::= C_FLT::c_flt_str_in(self, store_in);
      store_in.loc := sz;
      return(store_in); end;

   str_in(arg: FSTR, prec: INT): FSTR is
      store_in: FSTR;
      des_sz ::= prec+10;
      if (arg.size >= des_sz) then store_in := arg;
      else store_in := #FSTR(des_sz) end;
      sz ::= C_FLT::c_flt_str_in_prec(self, prec, store_in);
      store_in.loc := sz;
      return(store_in);  end;

   create (s: STR): SAME is
      return C_FLTD::atof(s).flt
   end;

   int:INT is
      -- INT version of self.  It is an error if self is not integral.
      -- Use truncate, floor, ceiling, or round to achieve this.
      return self.fltd.int;
   end;

   fltd:FLTD is
      -- An FLTD version of self.  Built-in.
      raise "FLT::fltd:FLTD undefined." end;

   fltx:FLTX is
      -- An extended floating point version of self. It is an
      -- error if the value cannot be held in a FLTX. Built-in.
      raise "FLT::fltx:FLTX undefined." end;

   fltdx:FLTDX is
      -- An extended floating point version of self.  Built-in.
      raise "FLT::fltdx:FLTDX undefined." end;

   truncate:SAME is -- The nearest integer toward zero.
      return self.fltd.truncate.flt;
   end;

   floor:SAME is -- The largest integer not greater than self.
      return self.fltd.floor.flt;
   end;

   ceiling:SAME is -- The smallest integer not less than self.
      return self.fltd.ceiling.flt;
   end;

   round:SAME is -- The closest integer to self.
      return self.fltd.round.flt;
   end;

   -- Special values.

   -- An approximation of the mathematical value "pi".
   const pi:SAME:=FLTD::pi.flt;

   -- An approximation of the base of the natural logarithms "e".
   const e:SAME:=FLTD::e.flt;

   const sqrt_2:SAME:=FLTD::sqrt_2.flt;                 -- Approximation of 2.sqrt.
   const log_2:SAME:=FLTD::log_2.flt;                   -- Approximation of 2.log.
   const log2_e:SAME:=FLTD::log2_e.flt;                 -- Approximation of e.log2.
   const log10_e:SAME:=FLTD::log10_e.flt;               -- Approximation of e.log10.
   const log_10:SAME:=FLTD::log_10.flt;                 -- Approximation of 10.log.
   const half_pi:SAME:=FLTD::half_pi.flt;               -- Approximation of pi/2.
   const quarter_pi:SAME:=FLTD::quarter_pi.flt;         -- Approximation of pi/4.
   const inv_sqrt_2:SAME:=FLTD::inv_sqrt_2.flt;         -- Approximation of 1/(2.sqrt).
   const inv_pi:SAME:=FLTD::inv_pi.flt;                 -- Approximation of 1/pi.
   const double_inv_pi:SAME:=FLTD::double_inv_pi.flt;   -- Approximation of 2/pi.
   const double_sqrt_pi:SAME:=FLTD::double_sqrt_pi.flt; -- Approximation of 2*(pi.sqrt).

   -- The value to be used to represent no element in sets.
   const nil:SAME:=signaling_NaN(0);

   signaling_NaN(sig:INT):SAME is
       -- IEEE signalling NaN.  `sig' is the significand (presently unused).
       return C_FLT::r_signaling_nan_(sig);
   end;

   quiet_NaN(sig:INT):SAME is
       -- IEEE quiet NaN.  `sig' is the significand (presently unused).
       return C_FLT::r_quiet_nan_(sig);
   end;

   infinity:SAME is -- IEEE Infinity.
      return C_FLT::r_infinity_;
   end;

   const epsilon:SAME:=1.19209290e-07; -- The minimum x>0.0 such that 1.0+x/=x.

   const digits:INT:=6; -- The number of decimal digits of precision.

   -- The number of bits in the significand, including an implied bit.
   const mantissa_bits:INT:=24;

   -- The smallest normalized positive number.
   const min_normal:SAME:=1.17549435e-38;

   -- The largest normalized positive number.
   const max_normal:SAME:=3.40282347e38;

   min_subnormal:SAME is -- The smallest subnormal positive number.
      return C_FLT::r_min_subnormal_;
   end;

   max_subnormal:SAME is -- The largest subnormal positive number.
      return C_FLT::r_max_subnormal_;
   end;

   -- The minimum negative integer x such that b^(x-1) is in the range
   -- of normalized floating point numbers.
   const min_exp:INT:=-125;

   -- The minimum x such that 10^x is in the range of normalized
   -- floating point numbers.
   const min_exp10:INT:=-37;

   const max_exp:INT:=128; -- The maximum allowable exponent.

   const max_exp10:INT:=38; -- The maximum x such that 10^x is within range.

   -- Useful iters

   sum!(i:SAME!):SAME is
      -- Yields the sum of all previous values of `i'.
      r::=0.0; loop r:=r+i; yield r end end;

   product!(i:SAME!):SAME is
      -- Yields the product of all previous values of `i'.
      r::=1.0; loop r:=r*i; yield r end end;

end; -- class FLT

-------------------------------------------------------------------
external class C_FLT is
    -- This corresponds to the standard math functions linkable with "-lm".

    ir_finite_(a:FLT):BOOL;
    ir_isinf_(a:FLT):BOOL;
    ir_isnan_(a:FLT):BOOL;
    ir_isnormal_(a:FLT):BOOL;
    ir_issubnormal_(a:FLT):BOOL;
    ir_signbit_(a:FLT):BOOL;
    ir_iszero_(a:FLT):BOOL;
    ir_ilogb_(a:FLT):INT;
    r_copysign_(a:FLT,b:FLT):FLT;
    r_nextafter_(a:FLT,b:FLT):FLT;
    r_remainder_(a:FLT,b:FLT):FLT;
    r_fmod_(a:FLT,b:FLT):FLT;
    r_scalbn_(a:FLT,b:INT):FLT;
    --r_c_flt_sincos_(a:FLT):TUP{FLT,FLT};
    --r_c_flt_sincospi_(a:FLT):TUP{FLT,FLT};
    r_signaling_nan_(a:INT):FLT;
    r_quiet_nan_(a:INT):FLT;
    r_infinity_:FLT;
    r_min_subnormal_:FLT;
    r_max_subnormal_:FLT;
   c_flt_str_in(f: FLT, s: FSTR): INT;
   c_flt_str_in_prec(f: FLT, prec: INT, s: FSTR): INT;

end;
-------------------------------------------------------------------