--  -----------------------------------------------------------------------
--  Title:       double_primitive_functions
--  Last Mod:    Fri Mar 29 16:08:12 1991
--  Author:      Vincent Broman <broman@nosc.mil>
--   Copyright 1990 Vincent Broman
--     Permission granted to copy, modify, or compile this software for
--   one's own use, provided that this copyright notice is preserved intact.
--     Permission granted to distribute compiled binary copies of this
--   software which are linked in with some other application.
--     Permission granted to distribute other copies of this software,
--   provided that (1) any copy which is not source code, i.e. not in the
--   form in which the software is usually maintained, must be accompanied
--   by a copy of the source code from which it was compiled, and (2) the
--   one distributing it must refrain from imposing on the recipient
--   further restrictions on the distribution of this software.
--   
--  Visibility:  
--  Description: 
--		 functions on floating point types whose implementation
--		 depends on access to the mantissa/exponent representation
--		 of the floating point number.  this includes
--		 integer/fraction operations.
--		 
--		 you want to optimize this code so that constant-folding 
--		 and dead-code elimination can win big.
--		 
--  Exceptions:  numeric_error upon overflow of function scale.
--  -----------------------------------------------------------------------
with float_representation;

use float_representation;

package body double_primitive_functions is
-- 
   use double_rep;			-- inside float_representation

   overflow_error: exception renames numeric_error;
   -- AI-00387 specifies constraint_error for floating overflows.

   -- values computed during package elaboration, constant thereafter
   bias: integer;
   maxexp: integer;
   half_machine_epsilon: constant double := 0.5 ** double'machine_mantissa;
   machine_small: double;
   machine_normal_small: double;	-- smallest nondenormal positive


   function exponent( x: double) return integer is
--
-- return the exponent k such that 1/2 <= x/(2**k) < 1,
-- or zero for x = 0.0 .
-- 
   begin
      if x = 0.0 then
      	 return 0;
      end if;
      -- nonzero x
      -- essentially, just extract the exponent field and subtract the bias
      declare
	 ff: float_fields := split_fields( x);
	 ex: integer := integer( ff.exp);
      begin
	 if gradual_underflow and then ex = 0 then
	    -- nonzero, denormalized.  normalize with a multiplication.
	    ff := split_fields( two_to_mant * x);
	    return integer( ff.exp) - bias - mant;
	 else
	    return ex - bias;
	 end if;
      end;
   end exponent;


   function mantissa (x: double) return double is
--
-- return scale( x, - exponent( x)) if x is nonzero,  0.0 otherwise.
--
   begin
      if x = 0.0 then
      	 return 0.0;
      end if;
      -- nonzero x
      -- just set the exponent field to be 0+bias
      declare
	 ff: float_fields := split_fields( x);
      begin
	 if gradual_underflow and then ff.exp = 0 then
	    -- nonzero, denormalized.  normalize with a multiplication.
	    ff := split_fields( two_to_mant * x);
	 end if;
	 ff.exp := biased_exp( bias);
	 return join_fields( ff);
      end;
   end mantissa;
      

   function scale (x: double;
		   k: integer) return double is
--
-- return x * 2**k quickly, or quietly underflow to zero,
-- or raise an exception on overflow.
-- 
   begin
      if x = 0.0 then
      	 return 0.0;
      end if;
      -- nonzero x.
      -- essentially, just increment the exponent,
      -- checking the hairy special cases.
      declare
	 ff: float_fields := split_fields( x);
	 ex: integer := integer( ff.exp);
	 incr: integer := k;
	 newex: integer;
      begin
	 if gradual_underflow and then ex = 0 then
	    -- nonzero, denormalized.  normalize with a multiplication.
	    ff := split_fields( two_to_mant * x);
	    incr := incr - mant;
	    ex := integer( ff.exp);
	 end if;

	 newex := ex + incr;

	 if newex <= 0 then
	    if gradual_underflow and then newex >= 2 - mant then
	       -- must produce denormal by division
	       ff.exp := biased_exp( newex + mant);
	       return join_fields( ff) / two_to_mant;
	    else
	       return 0.0;
	    end if;
	 elsif newex > maxexp then
	    raise overflow_error;
	 else
	    ff.exp := biased_exp( newex);
	    return join_fields( ff);
	 end if;
      end;
   end scale;


-- 
-- the leading_part function reduces precision to a specified
-- number of bits by truncation, not rounding.
-- 

   function leading_part( x: double;
      	       	     	  k: integer) return double is
-- 
   begin
      if k < 1 then
	 return 0.0;
      elsif k >= double'machine_mantissa then
	 return x;
      end if;
      declare
	 clear: integer := mant - k;
	 -- the number of least significant bits to clear to zero
	 ff: float_fields := split_fields( x);
      begin
	 -- now,  1 <= clear <= mant - 1
	 if clear < mantissa_word_size then
	    ff.mant_4 := ff.mant_4 and truncate_mask_wd( clear);
	    -- mant_3, mant_2, and mant_1 unchanged
	 elsif clear < 2 * mantissa_word_size then
	    ff.mant_4 := zero_word;
	    ff.mant_3 := ff.mant_3 and
	       	     	 truncate_mask_wd( clear - mantissa_word_size);
	    -- mant_2 and mant_1 unchanged
	 elsif clear < 3 * mantissa_word_size then
	    ff.mant_4 := zero_word;
	    ff.mant_3 := zero_word;
	    ff.mant_2 := ff.mant_2 and
	       	     	 truncate_mask_wd( clear - 2 * mantissa_word_size);
	    -- mant_1 unchanged
	 elsif clear < mant - 1 then
	    ff.mant_4 := zero_word;
	    ff.mant_3 := zero_word;
	    ff.mant_2 := zero_word;
	    ff.mant_1 := ff.mant_1 and
	       	     	 truncate_mask_hd( clear - 3 * mantissa_word_size);
	 else				-- clear = mant - 1
	    ff.mant_4 := zero_word;
	    ff.mant_3 := zero_word;
	    ff.mant_2 := zero_word;
	    ff.mant_1 := zero_head;
	 end if;
	 return join_fields( ff);
      end;
   end leading_part;


--
-- oddness is checked by examining one bit of the base 2
-- representation of x.  since leading zeroes, one leading one,
-- and trailing zeroes of the number are not present explicitly
-- in floating point, they are checked for with preliminary if tests.
--

   function odd( x: double) return boolean is
-- 
      ax: constant double := abs( x);

   begin
      if ax < 2.0 then
      	 return ax >= 1.0;
      else
	 declare
	    ff: float_fields := split_fields( ax);
	    ex: natural := natural( ff.exp) - bias - 2;
	    -- bigendian bit nbr in mantissa for 2**0
	    exoff: natural;
	    bitoff: mantissa_word_bitno;
	 begin
	    if ax >= two_to_mant then
	       return false;
	    -- now,  2 <= exponent( x) <= mant
	    elsif ex < mantissa_head_size then
	       return boolean( ff.mant_1( bitno_hd( ex)));
	    else
	       exoff := ex - mantissa_head_size;
	       bitoff := bitno_wd( exoff mod mantissa_word_size);
	       case exoff / mantissa_word_size is
		  when 0 =>
		     return boolean( ff.mant_2( bitoff));
		  when 1 =>
		     return boolean( ff.mant_3( bitoff));
		  when 2 =>
		     return boolean( ff.mant_4( bitoff));
		  when others =>
		     null;		-- no value causes program_error
	       end case;
	    end if;
	 end;
      end if;
   end odd;


--
-- truncation requires making all fractional bits zero in a signed
-- magnitude representation of floating point numbers.
-- for large numbers, this has no effect, as all large floats are integers.
-- for abs( x) < 1.0 this produces zero.
-- for middling numbers, we AND the mantissa bit array with a mask
-- corresponding to the number of bits we want cleared.
-- This is accomplished inside the leading_part function.
--

   function truncate (x: double) return double is
--
-- truncate x to the nearest integer value with absolute value
-- not exceeding abs( x).  No conversion to an integer type
-- is expected, so truncate cannot overflow for large arguments.
--
   begin
      return leading_part( x, exponent( x));
   end truncate;


--
-- all the other integer/fraction functions are based on truncate.
--

   function floor (x: double) return double is
--
-- return as a double the greatest integer value <= x.
--
   begin
      declare
	 tx: constant double := truncate( x); -- no constraint_error
      begin
	 if x >= 0.0 or x = tx then
	    return tx;
	 else
	    return tx - 1.0;		-- exactly
	 end if;
      end;
   end floor;


   function ceiling (x: double) return double is
--
-- return as a double the least integer value >= x.
--
   begin
      declare
	 tx: constant double := truncate( x); -- no constraint_error
      begin
	 if x <= 0.0 or x = tx then
	    return tx;
	 else
	    return tx + 1.0;		-- exactly
	 end if;
      end;
   end ceiling;


   function round (x: double) return double is
--
-- return as a double the integer value nearest x.
-- in case of a tie, prefer the even value.
--
   begin
      declare
	 tx: double := truncate( x);
	 diff: double := x - tx;	-- computed exactly
      begin
	 if abs( diff) < 0.5 then
	    return tx;
	 elsif diff > 0.5 then
	    return tx + 1.0;
	 elsif diff < -0.5 then
	    return tx - 1.0;
	 elsif diff = 0.5 then
	    if odd( tx) then
	       return tx + 1.0;
	    else
	       return tx;
	    end if;
	 else				-- diff = -0.5
	    if odd( tx) then
	       return tx - 1.0;
	    else
	       return tx;
	    end if;
	 end if;
      end;
   end round;


   function "rem"( x, y: double) return double is
--
-- returns the machine-representable value closest to
-- x - y * round( x/ y), evaluated with exact arithmetic.
-- The result is exact, except for y near zero,
-- and only if the implementation does not handle denormalized numbers.
-- numeric_error is raised if y = 0.0 .
-- 
   begin
      if y = 0.0 then
	 raise numeric_error;
      elsif x = 0.0 then
      	 return x;
      else
	 declare
	    ax: double := abs( x);
	    ay: double := abs( y);
	    delexp: integer := exponent( ax) - exponent( ay);
	    sy: double := scale( ay, delexp);
	    -- ax & sy now have same exponent
	    odd_mult: boolean := false;

	 begin
	    -- We repeatedly subtract (ay times a power of 2) from ax.
	    while delexp > 0 loop
	       if ax >= sy then
		  -- It is assumed that the floating point subtraction of
		  -- two positives a and b is exact if 1/2 < a/b < 2.
		  -- This condition holds for ax and sy.
		  ax := ax - sy;
	       end if;
	       delexp := delexp - 1;
	       sy := sy * 0.5;	-- might use scale, for speed or exactness
	    end loop;
	    
	    -- unroll last two iterations to handle the requirements
	    -- dealing with evenness and size
	    if delexp = 0 then
	       -- sy = ay
	       if ax >= sy then
	       	  ax := ax - sy;
		  odd_mult := true;
	       end if;
	       delexp := -1;
	       sy := sy * 0.5;
	    end if;

	    if delexp = -1 then
	       -- 0 <= ax < 2 * sy = ay
	       if odd_mult and ax = sy then
		  ax := - sy;
	       elsif ax > sy then
	       	  ax := ax - ay;
	       end if;
	    end if;

	    if ax = 0.0 then
	       return copy_sign( 0.0, x);
	    elsif x > 0.0 then
	       return ax;
	    else
	       return - ax;
	    end if;
	 end;
      end if;
   end "rem";


   function adjacent( x: double;
		      towards: double) return double is
--
-- returns x if towards = x, otherwise returns the float_point number
-- nearest to x in the direction toward towards.
--
   begin
      if towards = x then
      	 return x;
      elsif towards > x then
      	 return successor( x);
      else
      	 return predecessor( x);
      end if;
   end adjacent;


   function successor( x: double) return double is
--
-- return the next floating point number greater than x,
-- or overflow if x = double'base'last.
--
   begin
      if - machine_normal_small <= x and x < machine_normal_small then
      	 return x + machine_small;
      elsif x < 0.0 and mantissa( x) = 0.5 then
      	 return x + scale( half_machine_epsilon, exponent( x) - 1);
      elsif double'machine_overflows or else x < double'base'last then
      	 return x + scale( half_machine_epsilon, exponent( x));
      else
      	 raise overflow_error;
      end if;
   end successor;


   function predecessor( x: double) return double is
--
-- return the next floating point number less than x,
-- or overflow if x = double'base'first.
--
   begin
      if - machine_normal_small < x and x <= machine_normal_small then
      	 return x - machine_small;
      elsif x > 0.0 and mantissa( x) = 0.5 then
      	 return x - scale( half_machine_epsilon, exponent( x) - 1);
      elsif double'machine_overflows or else x > double'base'first then
      	 return x - scale( half_machine_epsilon, exponent( x));
      else
      	 raise overflow_error;
      end if;
   end predecessor;


   function copy_sign( value: double;
		       sign: double) return double is
-- 
-- returns abs( value) with the sign of sign, i.e.
-- when sign > 0.0 or sign = +0.0   returns abs( value),
-- when sign < 0.0 or sign = -0.0   returns - abs( value).
-- 
      vf: float_fields := split_fields( value);
      sf: constant float_fields := split_fields( sign);

   begin
      -- what happens if signed zeroes are constructible but forbidden?
      vf.sign := sf.sign;
      return join_fields( vf);
   end copy_sign;


   function multiply_sign( value: double;
      	       	     	   sign: double) return double is
-- 
-- returns value multiplied by the sign of sign, i.e.
-- when sign > 0.0 or sign = +0.0   returns value,
-- when sign < 0.0 or sign = -0.0   returns - value.
-- 
      vf: float_fields := split_fields( value);
      sf: constant float_fields := split_fields( sign);

   begin
      -- what happens if signed zeroes are constructible but forbidden?
      vf.sign := bit( vf.sign /= sf.sign);
      return join_fields( vf);
   end multiply_sign;


begin					-- double_primitive_functions
   -- initialize constants requiring some computation
   declare
      ff: constant float_fields := split_fields( 0.6);
   begin
      bias := integer( ff.exp);
   end;
   maxexp := double'machine_emax + bias;
   if gradual_underflow then
      -- depends on interpretation of 'machine_emin
      machine_normal_small := 2.0 ** (double'machine_emin - 1);
      machine_small := machine_normal_small *
      	       	       2.0 ** (1 - double'machine_mantissa);
   else
      machine_normal_small := 2.0 ** (double'machine_emin - 1);
      machine_small := machine_normal_small;
   end if;
end double_primitive_functions;

-- $Header: d_primitive_functions_b.a,v 3.23 91/03/29 16:08:36 broman Stab $