--  -----------------------------------------------------------------------
--  Title:       longest_float_support_functions
--  Last Mod:    Thu Apr  4 07:35:09 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.
--   
--  Description: 
--
--	a longest-float-precision-only implementation of functions
--	supporting an implementation of generic_elementary_functions.
-- 
--  Exceptions:  No checking done.  An accidental overflow could be caused
--		 by an argument way outside the spec domain of a function.
--  -----------------------------------------------------------------------
-- 
with math_constants,
     longest_float_primitive_functions;

use math_constants,
    longest_float_primitive_functions;

package body longest_float_support_functions is

   arctan_crit: constant longest_float
      	       	:= 1.5 * longest_float( longest_float'machine_radix)
		       	 ** (1 - longest_float'machine_mantissa);
   sqrt_eps: constant longest_float := longest_float( 1.0 / sqrt_2)
		      	       	       ** longest_float'machine_mantissa;
   cubert_eps: constant longest_float := longest_float( 1.0 / cbrt_2)
		      	       	       	 ** longest_float'machine_mantissa;
   exp_crit: constant longest_float := cbrt_3 * cubert_eps;
   log_crit: constant longest_float := sqrt_5 * sqrt_eps;
   -- literal following approximates cubert( 4959) / 2.
   tan_crit: constant longest_float := 8.52644 * cubert_eps;


   function reduced_logdel_2 (y: in longest_float) return longest_float is
-- 
-- return the base-2 logarithm of (1+y)/(1-y),
-- for y in the interval [ -3+sqrt(8) .. 3-sqrt(8) ].
-- actually, this implementation supports any y smaller than 1/2.
-- the relative error is less than two roundoffs in the final result.
-- note: reduced_logdel_2( 0.0) = 0.0 exactly.
-- 
      y2: constant longest_float := y * y;
      two_over_ln2: constant := 2.0 / nat_log_2;
      n: positive;
      term: longest_float;

   begin				-- reduced_logdel_2
      if y2 < log_crit then
      	 return two_over_ln2 * y * (1.0 + y2 / 3.0);
      else
      	 n := (3 - longest_float'machine_mantissa) / exponent( y2);
	 -- n >= 3   and
	 -- y**(2*n+2)/(2*n+3) <= (1/2)**machine_mantissa
	 term := 1.0 / longest_float( n + n + 1);
	 while n > 1 loop
	    n := n - 1;
	    term := term * y2 + 1.0 / longest_float( n + n + 1);
	 end loop;
	 return two_over_ln2 * y * (1.0 + term * y2);
      end if;
   end reduced_logdel_2;


   function reduced_log_2 (x: in longest_float) return longest_float is
-- 
-- return the base-2 logarithm of x, x in the interval [sqrt(1/2),sqrt(2)].
-- 
   begin
      return reduced_logdel_2( (x - 1.0) / (x + 1.0));
   end reduced_log_2;


   procedure log_2_parts (int_part: out longest_float;
			  frac_part: out longest_float;
			  x: in longest_float) is
-- 
-- return the integer and fractional part of the
-- base-2 logarithm of x.
-- int_part is a floating pt integer, and abs( frac_part) <= 1/2.
-- 
      ex: integer := exponent( x);
      ma: longest_float := mantissa( x);
      sqrt_half: constant := sqrt_2 / 2.0;

   begin
      if ma < sqrt_half then
      	 ma := 2.0 * ma;
	 ex := ex - 1;
      end if;
      --     sqrt_half <= ma < sqrt_2 
      -- and x = ma * 2 ** ex
      int_part := longest_float( ex);
      frac_part := reduced_log_2( ma);
   end log_2_parts;

   function log_2 (x: in longest_float) return longest_float is
-- 
-- return the base-2 logarithm of x, assuming x > 0.0
-- 
-- the value is exact for x an integer power of 2.
-- 
      n, f: longest_float;

   begin				-- log_2
      log_2_parts( n, f, x);
      return n + f;
   end log_2;


   function reduced_two_to_the (x: in longest_float) return longest_float is
-- 
-- return 2.0**x for x in the interval [-0.72, 0.72].
-- the interval would be [-1/2, 1/2] but we allow some slop.
-- the error is less than two roundoffs in the result.
-- 
      y: constant longest_float := x * nat_log_2;
      n: positive;
      term: longest_float;

   begin				-- reduced_two_to_the
      if abs( y) < exp_crit then
      	 return 1.0 + y * (1.0 + 0.5 * y);
      else
      	 n := (longest_float'machine_mantissa + 4) / (2 - exponent( y));
	 -- n >= 4   and
	 -- abs( y)**(n+1)/factorial( n+1) <= (1/2)**(machine_mantissa+1)
	 term := y / longest_float( n);
	 while n > 1 loop
	    n := n - 1;
	    term := y / longest_float( n) * (1.0 + term);
	 end loop;
	 return 1.0 + term;
      end if;
   end reduced_two_to_the;

   
   function two_to_the (x: in longest_float) return longest_float is
-- 
-- return 2.0**x, underflowing to zero for very negative x,
-- and overflowing for very positive x.
-- 
      ix: longest_float := floor( x + 0.5);
      frac: longest_float := x - ix;

   begin
      return scale( reduced_two_to_the( frac), integer( ix));
   end two_to_the;


   function reduced_tan( x: in longest_float) return longest_float is
--
-- return 2 * tan( x / 2) for x in [-1.4,1.4],
-- with a relative error less than two roundoffs.
-- the continued fraction for tangent is evaluated with a sufficient
-- number of terms to guarantee the truncation error be insignificant.
-- small angles are handled as a separate case.
--
      y2: constant longest_float := 0.25 * x * x;
      n: positive;
      term: longest_float;

   begin
      if y2 < tan_crit then
      	 return x / (1.0 - y2 / (3.0 - y2 / 5.0));
	 -- if y2 <= 1/2 then this errs in x / reduced_tan x by <= 4/4959*y2**3.
      else
      	 n := 1 + (longest_float'machine_mantissa + 4) / (4 - exponent( y2));
	 -- n >= 4   and
	 -- error bound y**(2*n)/factorial( 2*n) <= (1/2)**(machine_mantissa+1)
	 term := longest_float( n + n - 1);
	 while n > 1 loop
	    n := n - 1;
	    term := longest_float( n + n - 1) - y2 / term;
	 end loop;
	 return x / term;
      end if;
   end reduced_tan;


   function reduced_sin (x: in longest_float) return longest_float is
-- 
-- return sin x for x in the interval [-pi/4,pi/4].
-- relative error in the final result is at most three roundoffs.
-- 
      t: constant longest_float := reduced_tan( x);

   begin
      return t / (1.0 + 0.25 * t * t);
   end reduced_sin;


   function reduced_cos (x: in longest_float) return longest_float is
-- 
-- return cos x for x in the interval [-pi/4,pi/4].
-- relative error in the final result is at most three roundoffs.
-- 
      t: constant longest_float := reduced_tan( x);
      t2: constant longest_float := t * t;

   begin
      return 1.0 - 0.5 * t2 / (1.0 + 0.25 * t2);
   end reduced_cos;


   function reduced_arctan (x: in longest_float) return longest_float is
-- 
-- returns arctan x in radians from [-pi/4,pi/4] for x in [-1,1].
-- the relative error is at most two roundoffs.
-- 
      x2: constant longest_float := x * x;
      n: positive;
      term: longest_float;

   begin
      if x2 < arctan_crit then
      	 return x;
      else
	 -- evaluate the continued fraction for x / atan x.
	 -- n terms give a relative error bounded by (1/3*x**2)**n,
	 -- which must be <= machine_epsilon / 2.
      	 n := 1 + (longest_float'machine_mantissa - 1)
	       	/ (exponent( 3.0 / x2) - 1);
	 term := longest_float( n + n - 1);
	 while n > 1 loop
	    n := n - 1;
	    term := longest_float( n + n - 1)
	       	  + longest_float( n * n) * x2 / term;
	 end loop;
	 return x / term;
      end if;
   end reduced_arctan;


end longest_float_support_functions;

-- $Header: l_support_functions_b.a,v 3.24 91/04/04 13:22:39 broman Exp $