--  -----------------------------------------------------------------------
--  Title:       generic_algebraic_functions
--  Last Mod:    Fri Jun  7 13:16:31 1991
--  Author:      Vincent Broman <broman@nosc.mil>
--  Visibility:  general purpose utility
--  Description: 
--		 This package provides certain algebraic functions
--		 of general utility on floating point types, viz.
--		 square, cube, square_root, cube_root, nth_root,
--		 hypot, complement, sec_to_tan, tan_to_sec, min, and max,
--		 as these are individually defined below.
--		 
--		 The rules for accuracy, exceptions, etc are analogous
--		 to those for generic_elementary_functions, except as
--		 deviations are noted below.  Unless otherwise specified,
--		 the relative accuracy of the function results must be
--		 2.0 * float_type'base'epsilon or better.
--		 
--  Exceptions:  Argument_error is raised if and only if the argument to a
--		 function falls outside the mathematical domain of the function.
--		 The exception which signals floating point overflow
--		 is raised only if some of the values permitted by the
--		 relative accuracy requirements for a function result
--		 exceed float_type'safe_large.
--		 Only square and hypot may overflow.
--		 
--		 Otherwise, constraint_error is raised if and only if the
--		 value returned for a function violates a range constraint
--		 associated with float_type.
--  -----------------------------------------------------------------------
with float_types,
     single_algebraic_functions,
     double_algebraic_functions;

use float_types;

package body generic_algebraic_functions is
-- 

   package sgl_fun renames single_algebraic_functions;
   package dbl_fun renames double_algebraic_functions;
   
   
   function square( x: in float_type) return float_type is
-- 
-- Returns x * x, 
-- to an accuracy satisfying the Ada model of floating point numbers.
-- 
   begin
      return x * x;			-- no type conversions needed
   end square;

   
   function cube( x: in float_type) return float_type is
-- 
-- Returns x * x * x, 
-- to an accuracy satisfying the Ada model of floating point numbers.
-- 
   begin
      return x * x * x;
   end cube;


   function square_root( x: in float_type) return float_type is
-- 
-- Returns the nonnegative square root of x,
-- or raises argument_error for negative x.
-- square_root( 0.0) = 0.0
-- 
   begin
      if float_type'digits <= single'base'digits then
	 return float_type( sgl_fun.square_root( single( x)));
      else
	 return float_type( dbl_fun.square_root( double( x)));
      end if;
   end square_root;


   function cube_root( x: in float_type;
		       y: in float_type := 0.0) return float_type is
--
-- Returns (exactly) the greatest real root, t, of the cubic equation:
-- t**3 = xx + yy*t, for some xx and yy which approximate x and y
-- to a relative accuracy of TBD * float_type'base'epsilon.
-- cube_root( 0.0, y) = 0.0 for y <= 0.0 .
-- 
-- Function cube_root is discontinuous near the curve
-- (x/2)**2 = (y/3)**3, x <= 0.
-- 
-- It is conjectured that the cube_root function and the rational
-- operations in package standard are sufficient to accurately compute
-- real roots of any real polynomial of degree up to four in a finite number
-- of steps, the number of steps being independent of the accuracy required.
--
   begin
      if float_type'digits <= single'base'digits then
	 return float_type( sgl_fun.cube_root( single( x),
	       	     	      	       	       single( y)));
      else
	 return float_type( dbl_fun.cube_root( double( x),
	       	     	      	       	       double( y)));
      end if;
   end cube_root;


   function nth_root( x: in float_type;
		      n: in positive) return float_type is
--
-- Returns the nth root of x,
-- the positive root being chosen in case n is even.
-- nth_root( 0.0, n) = 0.0 and nth_root( 1.0, n) = 1.0, for all n.
-- Argument_error is raised iff n is even and x negative.
--
   begin
      if float_type'digits <= single'base'digits then
	 return float_type( sgl_fun.nth_root( single( x),        n));
      else
	 return float_type( dbl_fun.nth_root( double( x),        n));
      end if;
   end nth_root;


   function hypot( x, y: in float_type) return float_type is
-- 
-- Returns the nonnegative square root of x**2 + y**2,
-- overflowing only if the final result makes it necessary.
-- hypot( x, 0.0) = abs( x)
-- hypot( 0.0, y) = abs( y)
-- 
   begin
      if float_type'digits <= single'base'digits then
	 return float_type( sgl_fun.hypot( single( x),
	       	     	      	       	   single( y)));
      else
	 return float_type( dbl_fun.hypot( double( x),
	       	     	      	       	   double( y)));
      end if;
   end hypot;


   function hypot( x, y, z: in float_type) return float_type is
--
-- Returns the nonnegative square root of x**2 + y**2 + z**2,
-- overflowing only if the final result makes it necessary.
-- hypot( x, 0.0, 0.0) = abs( x)
-- hypot( 0.0, y, 0.0) = abs( y)
-- hypot( 0.0, 0.0, z) = abs( z)
-- 
   begin
      if float_type'digits <= single'base'digits then
	 return float_type( sgl_fun.hypot( single( x),
	       	     	      	       	   single( y),
					   single( z)));
      else
	 return float_type( dbl_fun.hypot( double( x),
	       	     	      	       	   double( y),
					   double( z)));
      end if;
   end hypot;
   

   function complement( x: in float_type) return float_type is
--
-- Returns the nonnegative square root of (1 - x**2) if abs( x) <= 1.0,
-- raising argument_error otherwise.
-- complement(  1.0) = 0.0
-- complement( -1.0) = 0.0
-- complement(  0.0) = 1.0
--
   begin
      if float_type'digits <= single'base'digits then
	 return float_type( sgl_fun.complement( single( x)));
      else
	 return float_type( dbl_fun.complement( double( x)));
      end if;
   end complement;


   function sec_to_tan( x: in float_type) return float_type is
--
-- Returns the nonnegative square root of (x**2 - 1) if abs( x) >= 1.0,
-- otherwise raises argument_error.
-- Must never overflow.
-- sec_to_tan(  1.0) = 0.0
-- sec_to_tan( -1.0) = 0.0
--
   begin
      if float_type'digits <= single'base'digits then
	 return float_type( sgl_fun.sec_to_tan( single( x)));
      else
	 return float_type( dbl_fun.sec_to_tan( double( x)));
      end if;
   end sec_to_tan;


   function tan_to_sec( x: in float_type) return float_type is
--
-- Returns the positive square root of (x**2 + 1).
-- Must never overflow.
-- tan_to_sec( 0.0) = 1.0
-- 
   begin
      if float_type'digits <= single'base'digits then
	 return float_type( sgl_fun.tan_to_sec( single( x)));
      else
	 return float_type( dbl_fun.tan_to_sec( double( x)));
      end if;
   end tan_to_sec;


   function min( x, y: in float_type) return float_type is
--
-- Returns either x or y, choosing the minimum of the two,
-- where minimum means that some point in the safe interval
-- for the number chosen is less than or equal to some point
-- in the safe interval for the other number.
-- If the intervals overlap, the choice is implementation dependent.
--
   begin
      if x <= y then
      	 return x;
      else
      	 return y;
      end if;
   end min;


   function min( x, y, z: in float_type) return float_type is
--
-- Returns either x or y or z, choosing the minimum of the three,
-- where minimum means that some point in the safe interval
-- for the number chosen is less than or equal to some point
-- in the safe intervals for each of the other two numbers.
-- If the intervals overlap, the choice is implementation dependent.
--
   begin
      return min( min( x, y), z);
   end min;


   function max( x, y: in float_type) return float_type is
--
-- Returns either x or y, choosing the maximum of the two,
-- where maximum means that some point in the safe interval
-- for the number chosen is greater than or equal to some point
-- in the safe interval for the other number.
-- If the intervals overlap, the choice is implementation dependent.
--
   begin
      if x >= y then
      	 return x;
      else
      	 return y;
      end if;
   end max;


   function max( x, y, z: in float_type) return float_type is
--
-- Returns either x or y or z, choosing the maximum of the three,
-- where maximum means that some point in the safe interval
-- for the number chosen is greater than or equal to some point
-- in the safe intervals for each of the other two numbers.
-- If the intervals overlap, the choice is implementation dependent.
--
   begin
      return max( max( x, y), z);
   end max;


end generic_algebraic_functions;

-- $Header: g2_algebraic_functions_b.a,v 3.25 91/06/07 13:17:44 broman Rel $