--  -----------------------------------------------------------------------
--  Title:       single_algebraic_functions
--  Last Mod:    Fri Jun  7 12:58:17 1991
--  Author:      Vincent Broman <broman@nosc.mil>
--  Visibility:  withed by generic_algebraic_functions
--		 and single_elementary_functions.
--  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 * single'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 single'safe_large.
--		 Only square and hypot may overflow.
--  -----------------------------------------------------------------------
with math_constants,
     single_primitive_functions,
     single_support_functions;

use math_constants,
    single_primitive_functions,
    single_support_functions;

package body single_algebraic_functions is
-- 

   mant: constant := single'machine_mantissa;
   -- the following thresholds may err on the small and large sides, resp.
   sqrt_2_epsilon: constant := 0.5 ** ((mant + 1) / 2);
   recip_sqrt_2_epsilon: constant := 2.0 ** ((mant + 1) / 2);

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

   
   function cube( x: in single) return single 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 single) return single is
-- 
-- Returns the nonnegative square root of x,
-- or raises argument_error for negative x.
-- square_root( 0.0) = 0.0
-- 
   begin
      if x > 0.0 then
	 declare
	    expo: constant integer := exponent( x);
	    scaling: integer;		-- for scaling of result
	    reduced_arg, estimate: single;
	    --
	    -- minimax approximation constants from Hart p192, SQRT 0290
	    -- adjusted for an interval of [1/2,2] instead of [1/4,1].
	    p1: constant :=  3.09033297;
	    p0: constant :=  1.0000233;
	    q0: constant :=  3.0903708;
	 begin
	    if expo >= 0 then
	       scaling := expo / 2;
	    else
	       scaling := (expo - 1) / 2;
	    end if;
	    reduced_arg := scale( x, - 2 * scaling);
	    -- 0.5 <= reduced_arg < 2.0
	    -- 
	    -- rational approx of order 1,1 to square_root( reduced_arg)
	    estimate := (p1 * reduced_arg + p0) / (reduced_arg + q0);
	    -- which has 2**-8.6 max relative error.
	    -- need only 2**-5.3 to support iteration to follow.
	    -- Now apply Heron's iteration twice
	    estimate := 0.5 * (estimate + reduced_arg / estimate);
	    estimate := 0.5 * (estimate + reduced_arg / estimate);
	    return scale( estimate, scaling);
	 end;
      elsif x = 0.0 then
      	 return x;
      else
	 raise argument_error;
      end if;
   end square_root;


   function cubert( x: in single) return single is
--
-- Returns the cube root of x.
-- cubert( 0.0) = 0.0
--
   begin
      if x = 0.0 then
	 return x;
      else
	 declare
	    expo: constant integer := exponent( x);
	    scaling: integer;		-- for scaling of result
	    reduced_arg, estimate: single;
	    --
	    -- minimax approximation constants from Hart p199, CBRT 0760
	    -- adjusted for an interval of [1/2,4] instead of [1/8,1].
	    p1: constant :=  2.3452798;
	    p0: constant :=  1.5873712;
	    q0: constant :=  2.95493372;
	 begin
	    if expo >= 0 then
	       scaling := expo / 3;
	    else
	       scaling := (expo - 2) / 3;
	    end if;
	    reduced_arg := abs( scale( x, - 3 * scaling));
	    -- 0.5 <= reduced_arg < 4.0
	    -- 
	    -- rational approx of order 1,1 to cubert( reduced_arg)
	    estimate := (p1 * reduced_arg + p0) / (reduced_arg + q0);
	    -- which has 2**-7.3 max relative error.
	    -- need only 2**-6.0 to support iteration to follow.
	    -- Now apply a Newton update twice
	    estimate := (2.0 * estimate +
	       	     	 reduced_arg / (estimate * estimate)) / 3.0;
	    estimate := (2.0 * estimate +
	       	     	 reduced_arg / (estimate * estimate)) / 3.0;
	    return scale( copy_sign( estimate, x), scaling);
	 end;
      end if;
   end cubert;


   function cube_root( x: in single;
		       y: in single := 0.0) return single 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 * single'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 y = 0.0 then
      	 return cubert( x);
      elsif x = 0.0 then
      	 if y > 0.0 then
	    return square_root( y);
	 else
	    return 0.0;
	 end if;
      else
      	 raise program_error;		-- unimplemented
      end if;
   end cube_root;


   function log_2 (x: in single) return single is
-- 
-- return the approximate base-2 logarithm of x, assuming x > 0.0
-- the value is exact for x an integer power of 2.
-- the approximant has a continuous first derivative everywhere.
-- the absolute approximation error is less than 20/27-log2(5/3), or 0.00378 .
-- 
      -- the following coefficients are 21 bit model numbers approximating
      -- p1= 4 -9/5 log2(e), p2= -7 +9/2 log2(e), p3= 4 -27/10 log2(e).
      p1: constant :=   367827.0 / 2.0 ** 18; -- approx  1.4031487
      p2: constant := -1065085.0 / 2.0 ** 21; -- approx -0.50787210
      p3: constant :=   219621.0 / 2.0 ** 21; -- approx  0.10472345
      -- 
      ex: integer := exponent( x) - 1;
      sma: single := scale( x, - ex) - 1.0;
      -- 0 <= sma < 1

   begin				-- log_2
      return single( ex) + sma * (p1 + sma * (p2 + sma * p3));
   end log_2;


   function two_to_the (x: in single) return single is
-- 
-- return 2.0**x, underflowing to zero for very negative x,
-- and overflowing for very positive x.
-- the value is exact for model integer x.
-- the approximant has a continuous first derivative everywhere.
-- the relative approximation error is less than 1-17/27 cubert4, or 0.000525 .
-- 
      -- the following coefficients are 21 bit model numbers approximating
      -- p1= 2 -17/9 ln2, p2= -5 +68/9 ln2, p3= 4 -17/3 ln2.
      p1: constant := 1448549.0 / 2.0 ** 21; -- approx 0.69072199
      p2: constant :=  124315.0 / 2.0 ** 19; -- approx 0.23711205
      p3: constant :=  151343.0 / 2.0 ** 21; -- approx 0.072165966
      -- 
      ix: single := floor( x);
      fr: single := x - ix;

   begin
      return scale( 1.0 + fr * (p1 + fr * (p2 + fr * p3)), integer( ix));
   end two_to_the;


   function nth_root( x: in single;
		      n: in positive) return single 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				-- nth_root
      if n = 1 or x = 0.0 or x = 1.0 then
      	 return x;
      elsif n = 2 then
      	 return square_root( x);
      elsif n = 3 then
      	 return cubert( x);
      elsif x < 0.0 and n / 2 * 2 = n then
	 raise argument_error;
      else
      	 declare
	    fn: constant single := single( n);
	    ex: integer := exponent( x);
	    ma: single := mantissa( abs( x));
	    sqrt_half: constant := sqrt_2 / 2.0;
	    l2ma, rm: single;
	    rs: integer;
	 begin
	    if ma < sqrt_half then
	       ex := ex - 1;
	       ma := ma * 2.0;
	    end if;
	    -- abs( x) = ma * 2.0 ** ex  and  sqrt_half <= ma < sqrt_2.
	    l2ma := reduced_log_2( ma);
	    rs := integer( (l2ma + single( ex)) / fn);
	    -- rs is the integer nearest log_2( the nth root)
	    rm := (l2ma + single( ex - n * rs)) / fn;
	    -- abs( rm) <= 1/2 plus some roundoff
	    return copy_sign( scale( reduced_two_to_the( rm), rs), x);
	 end;
      end if;
   end nth_root;


   function hypot( x, y: in single) return single 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)
-- 
      ax: single := abs( x);
      ay: single := abs( y);

   begin
      if ax = 0.0 then
      	 return ay;
      else
      	 if ax < ay then
	    declare
	       swi: constant single := ax;
	    begin
	       ax := ay;
	       ay := swi;
	    end;
	 end if;
	 -- now ax >= ay >= 0.0 and ax > 0.0
	 if ay / ax < sqrt_2_epsilon then
	    return ax;
	 else
	    declare
	       xoy: constant single := ax / ay;
	    begin
	       return ax + ay / (xoy + tan_to_sec( xoy));
	    end;
	 end if;
      end if;
   end hypot;


   function hypot( x, y, z: in single) return single 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)
-- 
      ax: single := abs( x);
      ay: single := abs( y);
      az: single := abs( z);

      procedure sort_two( a, b: in out single) is
-- sort a, b so that  a >= b.
      begin
      	 if a < b then
	    declare
	       swi: constant single := a;
	    begin
	       a := b;
	       b := swi;
	    end;
	 end if;
      end sort_two;
      pragma inline( sort_two);

   begin				-- hypot
      sort_two( ax, ay);
      sort_two( ax, az);
      return hypot( ax, hypot( ay, az));
   end hypot;
   

   function complement( x: in single) return single 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
--
      ax: constant single := abs( x);

   begin
      if ax < 0.5 then
	 declare
	    -- solve 1 - ax**2 = (1 - dels)**2 by Newton's method.
	    -- s=1-dels .  s is initialized high and always decreases, in fact,
	    -- (s - exact_s) = (olds - exact_s)**2 / (2 * olds) .
	    dels: single := 0.5 * ax * ax; -- init with >=6.6 bits of accuracy
	    -- The initial value of s=1-dels has relative error < .5**t
	    -- if ax**4/(1-ax**2) < .5**(t-3) approximately.
	    -- Each Newton step increases t to 2t+1.
	    -- One step  attains t=25 for ax < 0.2079
	    -- Two steps attains t=25 for ax < 0.5841
	 begin
	    dels := (ax - dels) * (ax + dels) / (2.0 * (1.0 - dels));
	    dels := (ax - dels) * (ax + dels) / (2.0 * (1.0 - dels));
	    return 1.0 - dels;
	 end;
      elsif ax <= 1.0 then
	 declare
	    omx: constant single := 1.0 - ax;
	 begin
	    return square_root( 2.0 * omx - omx * omx);
	 end;
      else
      	 raise argument_error;
      end if;
   end complement;


   function sec_to_tan( x: in single) return single 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
--
      ax: single := abs( x);

   begin
      if ax > 2.0 then
      	 return ax * complement( 1.0 / ax);
      elsif ax >= 1.0 then
      	 return square_root( (ax - 1.0) * (ax + 1.0));
      else
      	 raise argument_error;
      end if;
   end sec_to_tan;


   function tan_to_sec( x: in single) return single is
--
-- Returns the positive square root of (x**2 + 1).
-- Must never overflow.
-- tan_to_sec( 0.0) = 1.0
-- 
      ax: single := abs( x);

   begin
      if ax > recip_sqrt_2_epsilon then
	 -- this case separated out to avoid overflow in evaluation of
	 -- single'safe_large + 0.5/single'safe_large .
      	 return ax;
      else
      	 declare
	    base, compl, dels, size: single;
	    -- Solve (base + dels)**2 = x**2 + 1  for dels by Newton's method.
	    --   s = base + dels
	    --   base = max( 1, ax)
	    --   initial dels = 1/(2*ax)  for ax => 1
	    --                  ax**2 / 2 for ax <  1
	    --   compl = x**2 + 1 - base**2
	    -- Convergence is quadratic. In fact,
	    -- (s - exact_s) = (olds - exact_s)**2 / (2 * olds)
	    -- This implies that convergence is monotonic decreasing.

	    -- The initial value of s=base+dels has relative error < .5**t
	    -- if size = either ax or 1/ax 
	    -- satisfies size**4/(1+size**2) < .5**(t-3) approximately.
	    -- Each Newton step increases t to better than 2t+1.
	    -- One step  attains t=25 for ax or 1/ax < 0.2125
	    -- Two steps attains t=25 for ax or 1/ax < 0.7197
	    -- Three steps attains t=25 everywhere.
	 begin
	    if ax >= 1.0 then
	       base := ax;
	       compl := 1.0;
	       dels := 0.5 / ax;
	       size := 1.0 / ax;
	    else
	       base := 1.0;
	       compl := ax * ax;
	       dels := 0.5 * compl;
	       size := ax;
	    end if;
	    dels := (dels * dels + compl) / (2.0 * (base + dels));
	    dels := (dels * dels + compl) / (2.0 * (base + dels));
	    if size > 0.7197 then
	       dels := (dels * dels + compl) / (2.0 * (base + dels));
	    end if;
	    return base + dels;
	 end;
      end if;
   end tan_to_sec;


   function min( x, y: in single) return single 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 single) return single 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 single) return single 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 single) return single 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 single_algebraic_functions;

-- $Header: s_algebraic_functions_b.a,v 3.25 91/06/07 13:15:05 broman Rel $