home *** CD-ROM | disk | FTP | other *** search
/ ADA Programming Guide / ADA Programming Guide.iso / ada_gwu / adaed / math / math.adb < prev    next >
Encoding:
Text File  |  1996-01-30  |  61.3 KB  |  1,930 lines

  1. --  -----------------------------------------------------------------------
  2. --  Title:       Float_elementary_functions
  3. --  Last Mod:    Thu Apr 19 11:26:27 1990
  4. --  Author:      Vincent Broman <broman@nosc.mil>
  5. --   Copyright 1990 Vincent Broman
  6. --     Permission granted to copy, modify, or compile this software for
  7. --   one's own use, provided that this copyright notice is preserved intact.
  8. --     Permission granted to distribute compiled binary copies of this
  9. --   software which are linked in with some other application.
  10. --     Permission granted to distribute other copies of this software,
  11. --   provided that (1) any copy which is not source code, i.e. not in the
  12. --   form in which the software is usually maintained, must be accompanied
  13. --   by a copy of the source code from which it was compiled, and (2) the
  14. --   one distributing it must refrain from imposing on the recipient
  15. --   further restrictions on the distribution of this software.
  16. --   
  17. --  Description: 
  18. --
  19. --    a Float-precision-only implementation of 
  20. --    the proposed standard generic package of elementary functions, 
  21. --    from the ISO-IEC/JTC1/SC22/WG9 (Ada)
  22. --    Numerics Rapporteur Group proposal, Draft 1.1.
  23. --    
  24. --     This will evolve to conform to later drafts.
  25. -- 
  26. --  Exceptions:  argument_error, numeric_error are raised.
  27. --  -----------------------------------------------------------------------
  28. -- 
  29. with math_constants;  use math_constants;
  30. with primitive_functions;
  31. use  primitive_functions;
  32.  
  33. package body Math is
  34.    
  35.  
  36.    overflow_error: exception renames numeric_error;
  37. -- AI-00387 specifies constraint_error for floating overflows.
  38.  
  39.    fminexp: constant Float := Float( Float'machine_emin - 1);
  40.    fmaxexp: constant Float := Float( Float'machine_emax);
  41.  
  42.    mant: constant := Float'machine_mantissa;
  43.    recip_sqrt_epsilon: constant Float := 2.0 ** ((mant - 1) / 2);
  44.  
  45.    two_pi: constant := 2.0 * pi;
  46.    half_pi: constant := pi / 2.0;
  47.  
  48.    subtype octant_number is integer range 0 .. 7;
  49.  
  50.  
  51.    function sqr (x: Float) return Float is
  52. -- 
  53.    begin
  54.       return x * x;
  55.    end sqr;
  56.    pragma inline (sqr);
  57.    
  58.    
  59. ----------------------------------------------------------------------
  60. -- high precision argument reduction routines
  61.  
  62.  
  63.    procedure octant_reduce( x: in Float;
  64.                 r: out Float;
  65.                 nmod8: out octant_number) is
  66. --
  67. -- Octant_reduce reduces x to the interval [0, pi/4] with the high 
  68. -- absolute accuracy needed by trigonometric functions computed in radians
  69. -- which are subject to stringent relative accuracy requirements.
  70. -- 
  71. -- R is in 0.0 .. pi/4 and, to within the accuracy described below,
  72. --   either  (abs( x) - r)          is equal to pi/4 times an even integer
  73. --   or else (abs( x) - (pi/4 - r)) is equal to pi/4 times an odd integer.
  74. -- The angles in odd octants are complemented relative to pi/4 so that
  75. -- small sines and cosines may be computed from small angles,
  76. -- i.e. r is small anywhere near the x and y coordinate axes.
  77. -- Further, r=0 iff x=0, because of the way huge arguments get reduced.
  78. -- 
  79. -- NMOD8 is the octant number of x, in 0 .. 7, i.e.
  80. -- for x >= 0 nmod8 is     truncate( x / (pi/4)) mod 8,
  81. -- for x <  0 nmod8 is 7 - truncate( - x / (pi/4)) mod 8.
  82. -- 
  83. -- The relative accuracy of the reduced angle corresponding to (r,nmod8)
  84. -- is better than 2 ulp if abs( x) < 6433 < pi/4*2**13.
  85. -- 
  86. -- The inexactness of (r,nmod8) is caused entirely by
  87. --   (1) inexact reduction of huge arguments into the interval [-6433,6433].
  88. --   (2) the difference of 1.01e-21 between pi/4 and sum( pterm), 
  89. --   (3) the two roundings in producing r itself.
  90. -- If abs(x) < 6433, the absolute error in r caused by (1) and (2) 
  91. -- is less than  (sum(pterm)-pi/4)*abs(x)/(pi/4) < 8.3e-18 .
  92. -- 
  93. -- On (0, 22000) the maximum value of abs( x / (sin( x)*cos( x))) 
  94. -- for x a Float precision floating point number is less than 6.1e10,
  95. -- so the maximum relative error in r before rounding in (3) is less than
  96. --  6.1e10*1.01e-21/(pi/4) < 7.85e-11, which is .00066 ulp.
  97. -- 
  98. -- For arguments abs( x) > 6433, x is reduced and rounded more than
  99. -- once, so an extra absolute error is introduced, bounded by:
  100. --    max( 2**(exponent(2pi)-1-24), abs(x)*2**(-13+1-24)),
  101. -- where the former bound applies if reducing by whole periods produces a
  102. -- remainder small than 2pi, the latter bound otherwise.
  103. -- 
  104. -- Test your own system with 161*pi/2 to see how well argument reduction
  105. -- is done.  You should get cos(252.8982086181640625) =~ -4.1857068922e-9 .
  106. -- 
  107.       nbits: constant := 13;        -- significant bits in n
  108.       oneminus: constant := 1.0 - 3.0 * Float'epsilon;    -- 3 minimal?
  109.       toobig: constant := pi / 4.0 * (2.0 ** nbits) * oneminus; -- 6433.3+
  110.       p: constant Float := pi / 4.0;
  111.       ax: Float := abs( x);
  112.       n: Float;
  113.       odd_n: boolean;
  114.       
  115.       function truncate_nbits (x: Float) return Float is
  116. --
  117. -- truncate to an integer with nbits or fewer significant bits
  118. -- 
  119.       begin                -- truncate_nbits
  120.            if exponent( x) < nbits then
  121.         return truncate( x);
  122.      else
  123.         return shorten( x, nbits);
  124.      end if;
  125.       end truncate_nbits;
  126.       
  127.       function reduce_by_n (x: Float;
  128.                 n: Float) return Float is
  129. -- 
  130. -- compute x - n * sum(pterm(k),k=1..7) 
  131. -- erring by less than two roundoffs in the final result.
  132. -- x and n are assumed nonnegative. n has an nbit-length integer value.
  133. -- 
  134.            type softreal is array( 1 .. 7) of Float;
  135. --
  136. -- The high precision eighth-period:
  137. -- 
  138. -- sum(pterm(k),k=1..7) = pi/4 rounded to 68 bits, in 7 8-bit pieces.
  139.      pterm: constant softreal :=
  140.          (16#0.c9000000000000000#,
  141.           16#0.000fd000000000000#,
  142.           16#0.00000aa0000000000#,
  143.           16#0.00000002200000000#,
  144.           16#0.00000000016800000#,
  145.           16#0.000000000000c2000#,
  146.           16#0.0000000000000034c#);
  147. --
  148. -- each product of a 13 bit integer and an 8 bit pterm(k) computed below
  149. -- lands exactly on a 21 bit model number, no roundoff.
  150. --
  151.            diff: Float := x;
  152.      sumtail: Float := 0.0;
  153.      prod, nextdiff: Float;
  154.  
  155.       begin                -- reduce_by_n
  156.      for k in softreal'range loop
  157.         -- combine largest terms while the subtraction is still exact,
  158.         -- then sum tail to combine with head all at once.
  159.         prod := n * pterm( k);
  160.         nextdiff := diff - prod;
  161.         if abs( nextdiff) > prod then
  162.            -- nextdiff might be inexact
  163.            for m in reverse k + 1 .. softreal'last loop
  164.           sumtail := sumtail + n * pterm( m);
  165.            end loop;
  166.            return diff - (prod + sumtail); -- two roundings
  167.         end if;
  168.         diff := nextdiff;        -- difference is still exact
  169.      end loop;
  170.      return diff;
  171.       end reduce_by_n;
  172.  
  173.       function mod8 (x: Float) return octant_number is
  174. -- 
  175. -- return x mod 8 for x having a nonnegative integer value.
  176. -- 
  177.      xo8: Float := x * (1.0 / 8.0);
  178.  
  179.       begin                -- mod8
  180.      return octant_number( (xo8 - truncate( xo8)) * 8.0);
  181.       end mod8;
  182.  
  183.  
  184.    begin                -- octant_reduce
  185.       -- if ax too big to handle accurately then reduce by whole periods first
  186.       while ax > toobig loop
  187.      n := 8.0 * truncate_nbits( ax / (8.0 * p) * oneminus); -- n < ax / p
  188.      ax := reduce_by_n( ax, n);
  189.       end loop;
  190.       --
  191.       -- now ax < pi/4*2**nbits
  192.       -- 
  193.       n := truncate( ax / p);        -- 13 or fewer bits
  194.       odd_n := odd( n);
  195.       
  196.       if odd_n then
  197.      -- odd octant, complement angle wrt pi/4
  198.      -- p - (ax - n*p) = -ax + (n+1)*p
  199.      ax := - reduce_by_n( ax, n + 1.0);
  200.       elsif n > 0.0 then
  201.      -- for even nbr of octants, ax - n*p
  202.      ax := reduce_by_n( ax, n);
  203.       else                -- n = 0
  204.            null;
  205.       end if;
  206.       
  207.       if ax < 0.0 then
  208.            -- initial n was off by 1 because of roundoff in ax/p
  209.      r := - ax;
  210.      if odd_n then
  211.         n := n + 1.0;
  212.      else
  213.         n := n - 1.0;
  214.      end if;
  215.       elsif ax > p then
  216.            -- initial n was off by 1 because of roundoff in ax/p
  217.            -- want 2 * p - ax
  218.            r := - reduce_by_n( ax, 2.0);
  219.      if odd_n then
  220.         n := n - 1.0;
  221.      else
  222.         n := n + 1.0;
  223.      end if;
  224.       else
  225.            -- initial n correct, vast majority of cases
  226.            r := ax;
  227.       end if;
  228.       
  229.       if x >= 0.0 then
  230.      nmod8 := mod8( n);
  231.       else
  232.      nmod8 := octant_number'last - mod8( n);
  233.       end if;
  234.  
  235.    end octant_reduce;
  236.  
  237.  
  238.    procedure period_reduce( x: in Float;
  239.                 p: in Float;
  240.                 r: out Float;
  241.                 n: out octant_number) is
  242. -- 
  243. -- Period_reduce reduces x to the interval [0, p/8] with the high 
  244. -- absolute accuracy needed by trigonometric functions
  245. -- which are subject to stringent relative accuracy requirements.
  246. -- 
  247. -- X is the input angle, in units such that p is the cyclic period.
  248. -- Execution is erroneous if p is not positive.
  249. -- R is in 0.0 .. p/8 and, to within the accuracy described below,
  250. --   either  (abs( x) - r)          is equal to p/8 times an even integer
  251. --   or else (abs( x) - (p/8 - r))  is equal to p/8 times an odd integer.
  252. -- The angles in odd octants are complemented relative to p/8 so that
  253. -- small sines and cosines may be computed from small angles,
  254. -- i.e. r is small anywhere near the x and y coordinate axes.
  255. -- 
  256. -- N is the octant number of x, in 0 .. 7, i.e.
  257. -- for x >= 0 n is     truncate( x / (p/8)) mod 8,
  258. -- for x <  0 n is 7 - truncate( - x / (p/8)) mod 8.
  259. -- 
  260. -- The reduced angle corresponding to (r,n) is rounded and 
  261. -- correct up to the relative machine accuracy of r.
  262. -- All floating point arithmetic used is exact except possibly for 
  263. --  (1) complementing r in an odd octant, computing (p/8 - r1),
  264. --  (2) the last reductions, if p/8 loses bits by denormalization.
  265. -- 
  266.       r1: Float := abs( x);
  267.       exponent_diff: constant integer := exponent( r1) - exponent( p);
  268.       p1: Float := scale( p, exponent_diff); -- p1 & r1 now have same exponent
  269.       s: integer := exponent_diff + 3;
  270.  
  271.       octant: octant_number := 0;
  272.  
  273.    begin                -- period reduce
  274.       if s >= 0 then
  275.      -- We repeatedly subtract (p times a power of 2) from r.
  276.      loop
  277.         if r1 >= p1 then
  278.            -- It is assumed that the floating point subtraction of
  279.            -- two positives (a-b) is exact if either
  280.            -- (1) exponent(a) = exponent(b) or
  281.            -- (2) exponent(a)-1 = exponent(b) >= exponent(a-b).
  282.            -- This is the case for (r1-p1).
  283.            r1 := r1 - p1;
  284.            -- 
  285.            case s is
  286.           when 0 =>
  287.              r1 := p1 - r1;    -- odd octant, complement
  288.              octant := octant + 1;
  289.           when 1 =>
  290.              octant := octant + 2;
  291.           when 2 =>
  292.              octant := octant + 4;
  293.           when others =>
  294.              null;
  295.            end case;
  296.         end if; -- if r1 >= p1
  297.  
  298.      exit when s = 0;
  299.         s := s - 1;
  300.         p1 := p1 * 0.5;        -- could scale if that's faster
  301.  
  302.      end loop; -- down s
  303.  
  304.      r := r1;
  305.  
  306.       else
  307.      r := r1;            -- small angle, no reduction
  308.  
  309.       end if; -- if s nonnegative
  310.  
  311.       if x >= 0.0 then
  312.      n := octant;
  313.       else
  314.      n := octant_number'last - octant; -- no effect on odd_octant
  315.       end if;
  316.  
  317.    end period_reduce;
  318.  
  319.  
  320.    procedure ln2_reduce( x: in Float;
  321.              r: out Float;
  322.              n: out integer) is
  323. --
  324. -- Ln2_reduce reduces x to an interval of width approximately ln(2),
  325. -- with the high absolute accuracy needed by an exponential function which is
  326. -- subject to stringent relative accuracy requirements.
  327. -- 
  328. -- This routine requires abs( x) < 89.0 for spec accuracy.
  329. -- 
  330. -- Upon return,        x = r + n * ln(2),
  331. -- to within an absolute accuracy of 1 * epsilon (using the Ada model), 
  332. -- or more realistically to within 6e-8 .
  333. -- Also, abs( r) <= nat_log_2 / 2 + 256 * epsilon, 
  334. -- and n is in the range -128 .. 128,
  335. -- 
  336. -- The inexactness of the resultant (r,n) is caused entirely by 
  337. --   (1) the negligible difference of about 5.5e-14 between ln(2) and ln2, 
  338. --   (2) the 1 or 2 roundings in producing the result r itself.
  339. -- 
  340.       ln2_a: constant := 16#0.b1700000000#;
  341.       ln2_b: constant := 16#0.000217f0000#;
  342.       ln2_c: constant := 16#0.00000007d1c#;
  343.       -- ln(2) = ln2_a+ln2_b+ln2_c accurate to 43 bits, in 3 14-bit pieces.
  344.  
  345.       nrd: Float := floor( x / nat_log_2 + 0.5);
  346.  
  347.    begin                -- ln2_reduce
  348.       r := ((x - ln2_a * nrd) - ln2_b * nrd) - ln2_c * nrd;
  349.       -- The multiplications did happen exactly, using model numbers.
  350.       -- The a term subtracts exactly.  The c term is so small
  351.       -- that a bound on the absolute error can depend only on the size
  352.       -- of the result of the second subtraction.  So, two roundoffs
  353.       -- make the absolute error < 1*epsilon, pessimistically.
  354.       n := integer( nrd);
  355.       -- if we wanted to eliminate slop in the output interval,
  356.       -- then just test r, adjust n, and repeat the reduction.
  357.    end ln2_reduce;
  358.  
  359.  
  360.    function complement (x: in Float) return Float is
  361. -- compute sqrt(1-x**2)
  362.       ax: constant Float := abs( x);
  363.    begin
  364.       if ax >= 0.5 then
  365.            declare
  366.         t: constant Float := 1.0 - ax;
  367.      begin
  368.         return sqrt( 2.0 * t - t * t);
  369.      end;
  370.       else
  371.            return sqrt( (1.0 + x) * (1.0 - x));
  372.       end if;
  373.    end complement;
  374.  
  375.  
  376. ----------------------------------------------------------------------
  377. -- functions on reduced or simplified domains
  378.  
  379.  
  380.    function reduced_logdel_2 (y: in Float) return Float is
  381. -- 
  382. -- return the base-2 logarithm of (1+y)/(1-y),
  383. -- y in the interval [ -3+sqrt(8) .. 3-sqrt(8) ].
  384. -- the coefficients are from a rational (even) 1,1 approximation to
  385. -- log_2( (1+y)/(1-y) )/y  with relative error (omitting roundoff) < 1.91e-8.
  386. -- this delivers high relative accuracy for small y.
  387. -- Note: reduced_logdel_2( 0.0) = 0.0 exactly.
  388. -- 
  389.       p0: constant := -4.7845025;
  390.       p1: constant :=  1.2906108;
  391.       q0: constant := -1.6581823;
  392.       y2: constant Float := y * y;
  393.  
  394.    begin
  395.       return y * (p0 + p1 * y2) / (q0 + y2);
  396.    end reduced_logdel_2;
  397.  
  398.  
  399.    function reduced_log_2 (x: in Float) return Float is
  400. -- 
  401. -- return the base-2 logarithm of x, x in the interval [sqrt(1/2),sqrt(2)].
  402. -- 
  403.    begin
  404.       return reduced_logdel_2( (x - 1.0) / (x + 1.0));
  405.    end reduced_log_2;
  406.  
  407.    procedure log_2_parts (int_part: out Float;
  408.               frac_part: out Float;
  409.               x: in Float) is
  410. -- 
  411. -- return the integer and fractional part of the
  412. -- base-2 logarithm of x.
  413. -- int_part is a floating pt integer, and abs( frac_part) <= 1/2.
  414. -- 
  415.       ex: integer := exponent( x);
  416.       ma: Float := mantissa( x);
  417.       sqrt_half: constant := sqrt_2 / 2.0;
  418.  
  419.    begin
  420.       if ma < sqrt_half then
  421.            ma := 2.0 * ma;
  422.      ex := ex - 1;
  423.       end if;
  424.       --     sqrt_half <= ma < sqrt_2 
  425.       -- and x = ma * 2 ** ex
  426.       int_part := Float( ex);
  427.       frac_part := reduced_log_2( ma);
  428.    end log_2_parts;
  429.  
  430.    function log_2 (x: in Float) return Float is
  431. -- 
  432. -- return the base-2 logarithm of x, assuming x > 0.0
  433. -- 
  434. -- the value is exact for x an integer power of 2.
  435. -- 
  436.       n, f: Float;
  437.  
  438.    begin                -- log_2
  439.       log_2_parts( n, f, x);
  440.       return n + f;
  441.    end log_2;
  442.  
  443.  
  444.    function reduced_two_to_the (x: in Float) return Float is
  445. -- 
  446. -- return 2.0**x for x in the interval [-0.524, 0.524].
  447. -- the coefficients are a Remez-optimized rational(3,2) approximation,
  448. -- with relative error (omitting roundoff) < 1e-8.
  449. -- The interval would be [-1/2, 1/2] but we allow some slop.
  450. -- 
  451.       p0: constant :=  41.694187;
  452.       p1: constant :=  17.342325;
  453.       p2: constant :=  3.0047081;
  454.       p3: constant :=  0.23083161;
  455.       q0: constant :=  p0;              -- assures 2**0 = 1
  456.       q1: constant := -11.5578823;
  457.  
  458.    begin
  459.       return (p0 + x * (p1 + x * (p2 + x * p3)))
  460.        / (q0 + x * (q1 + x));
  461.    end reduced_two_to_the;
  462.  
  463.    
  464.    function two_to_the (x: in Float) return Float is
  465. -- 
  466. -- return 2.0**x, underflowing to zero for very negative x,
  467. -- and overflowing for very positive x.
  468. -- 
  469.       ix: Float := floor( x + 0.5);
  470.       frac: Float := x - ix;
  471.  
  472.    begin
  473.       return scale( reduced_two_to_the( frac), integer( ix));
  474.    end two_to_the;
  475.  
  476.  
  477.    function reduced_sin (x: in Float) return Float is
  478. -- 
  479. -- return sin x for x in the interval [0,pi/4].
  480. -- the coefficients are a chebychev-economized taylor series for
  481. -- sin((1-x)*pi/8)/((1-x)*pi/8) on [-1,1],
  482. -- with absolute error (omitting roundoff) < 2e-8.
  483. -- this gives high relative accuracy for small x.
  484. -- 
  485.       -- scale factor for size of interval [0,pi/4] already reckoned in below
  486.       a0: constant :=  9.74495337094254e-01;
  487.       a1: constant :=  1.28892138844516e-01;
  488.       a2: constant := -1.59024045807998e-01;
  489.       a3: constant := -1.28509163870500e-02;
  490.       a4: constant :=  7.83587003622433e-03;
  491.       a5: constant :=  4.55929708079455e-04;
  492.  
  493.       pio8: constant := pi / 8.0;
  494.       y: constant Float := pio8 - x;
  495.  
  496.    begin
  497.       return x * (a0 + y * (a1
  498.              + y * (a2
  499.              + y * (a3
  500.              + y * (a4
  501.              + y *  a5)))));
  502.    end reduced_sin;
  503.  
  504.  
  505.    function reduced_cos (x: in Float) return Float is
  506. -- 
  507. -- return cos x for x in the interval [0,pi/4].
  508. -- the coefficients are a chebychev-economized taylor series for
  509. -- cos((1-x)*pi/8) on [-1,1],
  510. -- with absolute error (omitting roundoff) < 1.8e-9
  511. -- 
  512.       -- scale factor for size of interval [0,pi/4] already reckoned in below
  513.       a0: constant :=  9.23879532410434e-01;
  514.       a1: constant :=  3.82683402033700e-01;
  515.       a2: constant := -4.61939745322990e-01;
  516.       a3: constant := -6.37789978530365e-02;
  517.       a4: constant :=  3.84943015536361e-02;
  518.       a5: constant :=  3.16859376215302e-03;
  519.       a6: constant := -1.27611586733778e-03;
  520.  
  521.       pio8: constant := pi / 8.0;
  522.       y: constant Float := pio8 - x;
  523.  
  524.    begin
  525.       return a0 + y * (a1
  526.         + y * (a2
  527.         + y * (a3
  528.         + y * (a4
  529.         + y * (a5
  530.         + y *  a6)))));
  531.    end reduced_cos;
  532.  
  533.  
  534.    function reduced_arctan (x: in Float) return Float is
  535. -- 
  536. -- returns arctan x in radians from [-pi/4,pi/4] for x in [-1,1].
  537. -- coefficents from Abramowitz and Stegun p81, #4.4.49 .
  538. -- relative error < 2e-8.
  539. -- 
  540.       x2: constant Float := x * x;
  541.       --
  542.       a2: constant :=  -0.3333314528;
  543.       a4: constant :=   0.1999355085;
  544.       a6: constant :=  -0.1420889944;
  545.       a8: constant :=   0.1065626393;
  546.       a10: constant := -0.0752896400;
  547.       a12: constant :=  0.0429096138;
  548.       a14: constant := -0.0161657367;
  549.       a16: constant :=  0.0028662257;
  550.  
  551.    begin
  552.       return x * (1.0 + x2 * ( a2 +
  553.             x2 * ( a4 +
  554.             x2 * ( a6 +
  555.             x2 * ( a8 +
  556.             x2 * (a10 +
  557.             x2 * (a12 +
  558.             x2 * (a14 +
  559.             x2 *  a16))))))));
  560.    end reduced_arctan;
  561.  
  562.  
  563. ----------------------------------------------------------------------
  564. -- implementations of user-callable fns
  565.  
  566.  
  567.    function sqrt (x : Float) return Float is
  568. -- 
  569. -- Declaration:
  570. --   function SQRT (X : FLOAT_TYPE) return FLOAT_TYPE;
  571. -- Definition:
  572. --   SQRT(X) ~ sqrt(X)
  573. -- Usage:
  574. --   Z := SQRT(X);
  575. -- Domain:
  576. --   X >= 0.0
  577. -- Range:
  578. --   SQRT(X) >= 0.0
  579. -- Accuracy:
  580. --   (a) Maximum relative error = 2.0 * FLOAT_TYPE'BASE'EPSILON
  581. --   (b) SQRT(0.0) = 0.0
  582. -- Notes:
  583. --   (a) The upper bound of the reachable range of SQRT is approximately given by
  584. --   SQRT(X)<=sqrt(FLOAT_TYPE'SAFE_LARGE)
  585. --   (b) Other standards might impose additional constraints on SQRT. For
  586. --   example, the IEEE standards for binary and radix-independent
  587. --   floating-point arithmetic require greater accuracy in the result of SQRT
  588. --   than this standard requires, and they require that SQRT(-0.0)=-0.0.
  589. --   An implementation of SQRT in GENERIC_ELEMENTARY_FUNCTIONS that
  590. --   conforms to this standard will conform to those other standards if it
  591. --   satisfies their additional constraints.
  592. -- 
  593.    begin
  594.       if x > 0.0 then
  595.      declare
  596.         expo: integer := exponent( x);
  597.         scaling: integer;        -- for scaling of result
  598.         reduced_arg, estimate: Float;
  599.         --
  600.         -- minimax approximation constants from Hart p192, SQRT 0290
  601.         -- adjusted for an interval of [1/2,2] instead of [1/4,1].
  602.         p1: constant :=  3.09033297;
  603.         p0: constant :=  1.0000233;
  604.         q0: constant :=  3.0903708;
  605.      begin
  606.         if expo >= 0 then
  607.            scaling := expo / 2;
  608.         else
  609.            scaling := (expo - 1) / 2;
  610.         end if;
  611.         reduced_arg := scale( x, - 2 * scaling);
  612.         -- 0.5 <= reduced_arg < 2.0
  613.         -- 
  614.         -- rational approx of order 1,1 to sqrt( reduced_arg)
  615.         estimate := (p1 * reduced_arg + p0) / (reduced_arg + q0);
  616.         -- which has 2**-8.6 max relative error.
  617.         -- need only 2**-5.3 to support iteration to follow.
  618.         -- Now apply Heron's iteration twice
  619.         estimate := 0.5 * (estimate + reduced_arg / estimate);
  620.         estimate := 0.5 * (estimate + reduced_arg / estimate);
  621.         return scale( estimate, scaling);
  622.      end;
  623.       elsif x = 0.0 then
  624.            return 0.0;
  625.       else
  626.      raise argument_error;
  627.       end if;
  628.    end sqrt;
  629.    
  630.    
  631.    function log (x : Float) return Float is
  632. --
  633. -- Declaration:
  634. --   function LOG (X : FLOAT_TYPE) return FLOAT_TYPE;
  635. -- Definition:
  636. --   LOG(X) ~ ln(X)
  637. -- Usage:
  638. --   Z := LOG(X);   -- natural logarithm
  639. -- Domain:
  640. --   X > 0.0
  641. -- Range:
  642. --   Mathematically unbounded
  643. -- Accuracy:
  644. --   (a) Maximum relative error = 4.0 * FLOAT_TYPE'BASE'EPSILON 
  645. --   (b) LOG(1.0) = 0.0
  646. -- Notes:
  647. --   The reachable range of LOG is approximately given by
  648. --   ln(FLOAT_TYPE'SAFE_SMALL) <= LOG(X) <= ln(FLOAT_TYPE'SAFE_LARGE)
  649. --
  650.    begin
  651.       if x > 0.0 then
  652.      -- log_2( 1.0) = 0.0
  653.      return log_2( x) * nat_log_2;
  654.       else
  655.      raise argument_error;
  656.       end if;
  657.    end log;
  658.    
  659.    
  660.    function log (x, base : Float) return Float is
  661. --
  662. -- Declaration:
  663. --   function LOG (X, BASE : FLOAT_TYPE) return FLOAT_TYPE;
  664. -- Definition:
  665. --   LOG(X,BASE) ~ log to base BASE of X
  666. -- Usage:
  667. --   Z := LOG(X, 10.0);   -- base 10 logarithm
  668. --   Z := LOG(X, 2.0);    -- base 2 logarithm
  669. --   Z := LOG(X, BASE);   -- base BASE logarithm
  670. -- Domain:
  671. --   (a) X > 0.0
  672. --   (b) BASE > 0.0
  673. --   (c) BASE /= 1.0
  674. -- Range:
  675. --   Mathematically unbounded
  676. -- Accuracy:
  677. --   (a) Maximum relative error = 4.0 * FLOAT_TYPE'BASE'EPSILON 
  678. --   (b) LOG(1.0,BASE) = 0.0
  679. -- Notes:
  680. --   (a) When BASE > 1.0, the reachable range of LOG is approximately given by
  681. --       log to base BASE of FLOAT_TYPE'SAFE_SMALL <= LOG(X,BASE) <=
  682. --       log to base BASE of FLOAT_TYPE'SAFE_LARGE
  683. --   (b) When 0.0 < BASE < 1.0, the reachable range of LOG is approximately given by
  684. --       log to base BASE of FLOAT_TYPE'SAFE_LARGE <= LOG(X,BASE) <=
  685. --       log to base BASE of FLOAT_TYPE'SAFE_SMALL
  686. --
  687.    begin
  688.       if x > 0.0 and base > 0.0 and base /= 1.0 then
  689.      -- log_2( 1.0) = 0.0
  690.      return log_2( x) / log_2( base);
  691.       else
  692.      raise argument_error;
  693.       end if;
  694.    end log;
  695.    
  696.  
  697.    function exp (x : Float) return Float is
  698. --
  699. -- Declaration:
  700. --   function EXP (X : FLOAT_TYPE) return FLOAT_TYPE;
  701. -- Definition:
  702. --   EXP(X) ~ e raised to the X power
  703. -- Usage:
  704. --   Z := EXP(X);   -- e raised to the power X
  705. -- Domain:
  706. --   Mathematically unbounded
  707. -- Range:
  708. --   EXP(X) >= 0.0
  709. -- Accuracy:
  710. --   (a) Maximum relative error = 4.0 * FLOAT_TYPE'BASE'EPSILON 
  711. --   (b) EXP(0.0) = 1.0
  712. -- Notes:
  713. --   The usable domain of EXP is approximately given by
  714. --   X <= ln(FLOAT_TYPE'SAFE_LARGE)
  715. --
  716.       minarg: constant Float := fminexp * nat_log_2;
  717.       maxarg: constant Float := fmaxexp * nat_log_2;
  718.  
  719.    begin
  720.       if x >= maxarg then
  721.            raise overflow_error;
  722.       elsif x < minarg then
  723.            return 0.0;
  724.       else
  725.            declare
  726.         r: Float;
  727.         n: integer;
  728.      begin
  729.         ln2_reduce( x, r, n);
  730.         return scale( reduced_two_to_the( r / nat_log_2), n);
  731.      end;
  732.       end if;
  733.    end exp;
  734.    
  735.    
  736.    function "**" (x, y : Float) return Float is
  737. --
  738. -- Declaration:
  739. --   function "**" (X, Y : FLOAT_TYPE) return FLOAT_TYPE;
  740. -- Definition:
  741. --   X ** Y ~ X raised to the power Y
  742. -- Usage:
  743. --   Z := X ** Y;   -- X raised to the power Y
  744. -- Domain:
  745. --   (a) X >= 0.0
  746. --   (b) Y > 0.0 when X = 0.0 
  747. -- Range:
  748. --   X ** Y >= 0.0
  749. -- Accuracy:
  750. --   (a) Maximum relative error (when X > 0.0) =
  751. --   (4.0+|Y*ln(X)|/32.0) * FLOAT_TYPE'BASE'EPSILON
  752. --   (b) X ** 0.0 = 1.0 when X > 0.0 
  753. --   (c) 0.0 ** Y = 0.0 when Y > 0.0
  754. --   (d) X ** 1.0 = X
  755. --   (e) 1.0 ** Y =1.0
  756. -- Notes:
  757. --   The usable domain of "**", when X > 0.0, is approximately the set of
  758. --   values for X and Y satisfying
  759. --   Y*ln(X) <= ln(FLOAT_TYPE'SAFE_LARGE)
  760. --   This imposes a positive upper bound on Y (as a function of X) when 
  761. --   X > 1.0, and a negative lower bound on Y (as a function of X) when 
  762. --   0.0 < X < 1.0.
  763. --
  764.       function half_precision (r: Float) return Float is
  765.      -- truncate to half as many bits in the mantissa
  766.       begin
  767.      return shorten( r, Float'machine_mantissa / 2);
  768.       end half_precision;
  769.       pragma inline( half_precision);
  770.  
  771.    begin                -- function "**"
  772.       if x = 1.0 then
  773.      return 1.0;            -- 1**y
  774.       elsif x > 0.0 then
  775.      if y = 0.0 then
  776.         return 1.0;            -- x**0
  777.      elsif y = 1.0 then
  778.         return x;            -- x**1
  779.      else
  780.         declare            -- general case of 2**( y*log_2( x))
  781.            sy: Float := y;
  782.            fex, ma, logx: Float;
  783.         begin
  784.            log_2_parts( fex, ma, x);
  785.            logx := fex + ma;
  786.            -- If y*log_2( x) overflows positively, we don't care
  787.            -- at what point the overflow exception happens.
  788.            -- Large negative values of y*log_2( x) are detected by
  789.            -- the following conservative test, which is equivalent to
  790.            -- sy * logx < 2 * fminexp, but will not overflow,
  791.            -- because abs( logx) < 2 * abs( fminexp).
  792.            if (sy / (2.0 * fminexp)) * logx > 1.0 then
  793.                  return 0.0;        -- underflow
  794.            else
  795.                  declare
  796.              -- y*log_2( x) must have 7 extra bits of accuracy 
  797.              -- to get full precision in the mantissa of the result.
  798.              -- so we split up the factors in two parts
  799.              y1: constant Float := half_precision( sy);
  800.              y2: constant Float := sy - y1;
  801.              -- sy = y1 + y2 exactly
  802.              w1: constant Float := half_precision( logx);
  803.              w2: constant Float := (fex - w1) + ma;
  804.              -- fex + ma = w1 + w2 exactly
  805.              -- 
  806.              -- what we need is (y1+y2)*(w1+w2) very accurately
  807.              y1w1: constant Float := y1*w1; -- exactly
  808.              n: constant Float := floor( y1w1 + 0.5);
  809.              frac: constant Float := (y1w1 - n) +
  810.                                             (sy * w2 + y2 * w1);
  811.              -- rounded to within 2**-Float'mantissa,
  812.              -- frac is nearly confined to [-.5,.5] .
  813.           begin
  814.              return scale( reduced_two_to_the( frac), integer( n));
  815.           end;
  816.            end if;
  817.         end;
  818.      end if;
  819.       elsif x = 0.0 and y > 0.0 then
  820.      return 0.0;            -- 0**y
  821.       else
  822.      raise argument_error;        -- x < 0 or x = 0 and y <= 0
  823.       end if;
  824.    end "**";
  825.  
  826.  
  827.    function sin (x : Float) return Float is
  828. --
  829. -- Declaration:
  830. --   function SIN (X : FLOAT_TYPE) return FLOAT_TYPE;
  831. -- Definition:
  832. --   SIN(X) ~ sin(X)
  833. -- Usage:
  834. --   Z := SIN(X);   -- X in radians
  835. -- Domain:
  836. --   Mathematically unbounded
  837. -- Range:
  838. --   |SIN(X)| <= 1.0
  839. -- Accuracy:
  840. --   (a) Maximum relative error = 2.0 * FLOAT_TYPE'BASE'EPSILON 
  841. --   when |X| is less than or equal to some documented implementation-dependent
  842. --   threshold, which must not be less than 
  843. --     FLOAT_TYPE'MACHINE_RADIX ** (FLOAT_TYPE'MACHINE_MANTISSA/2)
  844. --   For larger values of |X|, degraded accuracy is allowed. An implementation
  845. --   must document its behavior for large |X|. 
  846. --   (b) SIN(0.0) = 0.0
  847. --
  848.       r: Float;
  849.       oct: octant_number;
  850.  
  851.    begin
  852.       octant_reduce( x, r, oct);
  853.       case oct is
  854.            when 0| 3 =>
  855.         return reduced_sin( r);
  856.      when 1| 2 =>
  857.         return reduced_cos( r);
  858.      when 4| 7 =>
  859.         return - reduced_sin( r);
  860.      when 5| 6 =>
  861.         return - reduced_cos( r);
  862.       end case;
  863.    end sin;
  864.    
  865.    
  866.    function sin (x, cycle : Float) return Float is
  867. --
  868. -- Declaration:
  869. --   function SIN (X, CYCLE : FLOAT_TYPE) return FLOAT_TYPE;
  870. -- Definition:
  871. --   SIN(X,CYCLE) ~ sin(2Pi * X/CYCLE)
  872. -- Usage:
  873. --   Z := SIN(X, 360.0);   -- X in degrees
  874. --   Z := SIN(X, 1.0);     -- X in bams (binary angular measure)
  875. --   Z := SIN(X, CYCLE);   -- X in units such that one complete
  876. --                         -- cycle of rotation corresponds to
  877. --                         -- X = CYCLE
  878. -- Domain:
  879. --   (a) X mathematically unbounded
  880. --   (b) CYCLE > 0.0
  881. -- Range:
  882. --   |SIN(X,CYCLE)| <= 1.0
  883. -- Accuracy:
  884. --   (a) Maximum relative error = 2.0 * FLOAT_TYPE'BASE'EPSILON 
  885. --   (b) For integer k, SIN(X,CYCLE)= 0.0 when X=k*CYCLE/2.0
  886. --                                    1.0 when X=(4k+1)*CYCLE/4.0
  887. --                                    -1.0 when X=(4k+3)*CYCLE/4.0
  888. --
  889.       r: Float;
  890.       oct: octant_number;
  891.  
  892.    begin
  893.       if cycle > 0.0 then
  894.      period_reduce( x, cycle, r, oct);
  895.      r := r / cycle * two_pi;    -- underflow denormal allowed
  896.      case oct is
  897.         when 0| 3 =>
  898.            return reduced_sin( r);
  899.         when 1| 2 =>
  900.            return reduced_cos( r);
  901.         when 4| 7 =>
  902.            return - reduced_sin( r);
  903.         when 5| 6 =>
  904.            return - reduced_cos( r);
  905.      end case;
  906.       else
  907.      raise argument_error;
  908.       end if;
  909.    end sin;
  910.    
  911.    
  912.    function cos (x : Float) return Float is
  913. --
  914. -- Declaration:
  915. --   function COS (X : FLOAT_TYPE) return FLOAT_TYPE;
  916. -- Definition:
  917. --   COS(X) ~ cos(X)
  918. -- Usage:
  919. --   Z := COS(X);   -- X in radians
  920. -- Domain:
  921. --   Mathematically unbounded
  922. -- Range:
  923. --   |COS(X)| <= 1.0
  924. -- Accuracy:
  925. --   (a) Maximum relative error = 2.0 * FLOAT_TYPE'BASE'EPSILON 
  926. --   when |X| is less than or equal to some documented implementation-dependent
  927. --   threshold, which must not be less than 
  928. --     FLOAT_TYPE'MACHINE_RADIX ** (FLOAT_TYPE'MACHINE_MANTISSA/2)
  929. --   For larger values of |X|, degraded accuracy is allowed. An implementation
  930. --   must document its behavior for large |X|. 
  931. --   (b) COS(0.0) = 1.0
  932. --
  933.       r: Float;
  934.       oct: octant_number;
  935.  
  936.    begin
  937.       octant_reduce( x, r, oct);
  938.       case oct is
  939.      when 0| 7 =>
  940.         return reduced_cos( r);
  941.            when 1| 6 =>
  942.         return reduced_sin( r);
  943.      when 2| 5 =>
  944.         return - reduced_sin( r);
  945.      when 3| 4 =>
  946.         return - reduced_cos( r);
  947.       end case;
  948.    end cos;
  949.    
  950.    
  951.    function cos (x, cycle : Float) return Float is
  952. --
  953. -- Declaration:
  954. --   function COS (X, CYCLE : FLOAT_TYPE) return FLOAT_TYPE;
  955. -- Definition:
  956. --   COS(X,CYCLE) ~ cos(2Pi*X/CYCLE)
  957. -- Usage:
  958. --   Z := COS(X, 360.0);   -- X in degrees
  959. --   Z := COS(X, 1.0);     -- X in bams
  960. --   Z := COS(X, CYCLE);   -- X in units such that one complete
  961. --                         -- cycle of rotation corresponds to
  962. --                         -- X = CYCLE
  963. -- Domain:
  964. --   (a) X mathematically unbounded
  965. --   (b) CYCLE > 0.0
  966. -- Range:
  967. --   |COS(X,CYCLE)| <= 1.0
  968. -- Accuracy:
  969. --   (a) Maximum relative error = 2.0 * FLOAT_TYPE'BASE'EPSILON 
  970. --   (b) For integer k, COS(X,CYCLE) = 1.0 when X=k*CYCLE
  971. --                                     0.0 when X=(2k+1)*CYCLE/4.0
  972. --                                     -1.0 when X=(2k+1)*CYCLE/2.0 
  973. --
  974.       r: Float;
  975.       oct: octant_number;
  976.  
  977.    begin
  978.       if cycle > 0.0 then
  979.      period_reduce( x, cycle, r, oct);
  980.      r := r / cycle * two_pi;    -- underflow denormal allowed
  981.      case oct is
  982.         when 0| 7 =>
  983.            return reduced_cos( r);
  984.         when 1| 6 =>
  985.            return reduced_sin( r);
  986.         when 2| 5 =>
  987.            return - reduced_sin( r);
  988.         when 3| 4 =>
  989.            return - reduced_cos( r);
  990.      end case;
  991.       else
  992.      raise argument_error;
  993.       end if;
  994.    end cos;
  995.    
  996.    
  997.    function tan (x : Float) return Float is
  998. --
  999. -- Declaration:
  1000. --   function TAN (X : FLOAT_TYPE) return FLOAT_TYPE;
  1001. -- Definition:
  1002. --   TAN(X) ~ tan(X)
  1003. -- Usage:
  1004. --   Z := TAN(X);   -- X in radians
  1005. -- Domain:
  1006. --   Mathematically unbounded
  1007. -- Range:
  1008. --   Mathematically unbounded
  1009. -- Accuracy:
  1010. --   (a) Maximum relative error = 4.0 * FLOAT_TYPE'BASE'EPSILON 
  1011. --   when |X| is less than or equal to some documented implementation-dependent
  1012. --   threshold, which must not be less than 
  1013. --     FLOAT_TYPE'MACHINE_RADIX ** (FLOAT_TYPE'MACHINE_MANTISSA/2)
  1014. --   For larger values of |X|, degraded accuracy is allowed. An implementation
  1015. --   must document its behavior for large |X|. 
  1016. --   (b) TAN(0.0) = 0.0
  1017. --
  1018.       r: Float;
  1019.       oct: octant_number;
  1020.  
  1021.    begin
  1022.       octant_reduce( x, r, oct);
  1023.       case oct is
  1024.      when 0| 4 =>
  1025.         return reduced_sin( r) / reduced_cos( r);
  1026.            when 1| 5 =>
  1027.         if r /= 0.0 then
  1028.            return reduced_cos( r) / reduced_sin( r);
  1029.         else            -- x must have been huge
  1030.            return 0.0;
  1031.         end if;
  1032.      when 2| 6 =>
  1033.         if r /= 0.0 then
  1034.            return - reduced_cos( r) / reduced_sin( r);
  1035.         else            -- x must have been huge
  1036.            return 0.0;
  1037.         end if;
  1038.      when 3| 7 =>
  1039.         return - reduced_sin( r) / reduced_cos( r);
  1040.       end case;
  1041.    end tan;
  1042.    
  1043.    
  1044.    function tan (x, cycle : Float) return Float is
  1045. --
  1046. -- Declaration:
  1047. --   function TAN (X, CYCLE : FLOAT_TYPE) return FLOAT_TYPE;
  1048. -- Definition:
  1049. --   TAN(X,CYCLE) ~ tan(2Pi*X/CYCLE)
  1050. -- Usage:
  1051. --   Z := TAN(X, 360.0);   -- X in degrees
  1052. --   Z := TAN(X, 1.0);     -- X in bams
  1053. --   Z := TAN(X, CYCLE);   -- X in units such that one complete
  1054. --                         -- cycle of rotation corresponds to
  1055. --                         -- X = CYCLE
  1056. -- Domain:
  1057. --   (a) X /= (2k+1)*CYCLE/4.0, for integer k 
  1058. --   (b) CYCLE > 0.0
  1059. -- Range:
  1060. --   Mathematically unbounded
  1061. -- Accuracy:
  1062. --   (a) Maximum relative error = 4.0 * FLOAT_TYPE'BASE'EPSILON 
  1063. --   (b) TAN(X,CYCLE) = 0.0 when X=k*CYCLE/2.0, for integer k 
  1064. --
  1065.       r: Float;
  1066.       oct: octant_number;
  1067.  
  1068.    begin
  1069.       if cycle > 0.0 then
  1070.      period_reduce( x, cycle, r, oct);
  1071.      r := r / cycle * two_pi;    -- underflow denormal allowed
  1072.      case oct is
  1073.         when 0| 4 =>
  1074.            return reduced_sin( r) / reduced_cos( r);
  1075.         when 1| 5 =>
  1076.            if r /= 0.0 then
  1077.           return reduced_cos( r) / reduced_sin( r);
  1078.            else
  1079.           raise argument_error;
  1080.            end if;
  1081.         when 2| 6 =>
  1082.            if r /= 0.0 then
  1083.           return - reduced_cos( r) / reduced_sin( r);
  1084.            else
  1085.           raise argument_error;
  1086.            end if;
  1087.         when 3| 7 =>
  1088.            return - reduced_sin( r) / reduced_cos( r);
  1089.      end case;
  1090.       else
  1091.      raise argument_error;
  1092.       end if;
  1093.    end tan;
  1094.    
  1095.    
  1096.    function cot (x : Float) return Float is
  1097. --
  1098. -- Declaration:
  1099. --   function COT (X : FLOAT_TYPE) return FLOAT_TYPE;
  1100. -- Definition:
  1101. --   COT(X) ~ cot(X)
  1102. -- Usage:
  1103. --   Z := COT(X);   -- X in radians
  1104. -- Domain:
  1105. --   X /= 0.0
  1106. -- Range:
  1107. --   Mathematically unbounded
  1108. -- Accuracy:
  1109. --   (a) Maximum relative error = 4.0 * FLOAT_TYPE'BASE'EPSILON 
  1110. --   when |X| is less than or equal to some documented implementation-dependent
  1111. --   threshold, which must not be less than 
  1112. --     FLOAT_TYPE'MACHINE_RADIX ** (FLOAT_TYPE'MACHINE_MANTISSA/2)
  1113. --   For larger values of |X|, degraded accuracy is allowed. An implementation
  1114. --   must document its behavior for large |X|. 
  1115. --
  1116.       r: Float;
  1117.       oct: octant_number;
  1118.  
  1119.    begin
  1120.       if x /= 0.0 then
  1121.      octant_reduce( x, r, oct);
  1122.      case oct is
  1123.         when 0| 4 =>
  1124.            if r /= 0.0 then
  1125.           return reduced_cos( r) / reduced_sin( r);
  1126.            else            -- x must have been huge
  1127.           return 0.0;
  1128.            end if;
  1129.         when 1| 5 =>
  1130.            return reduced_sin( r) / reduced_cos( r);
  1131.         when 2| 6 =>
  1132.            return - reduced_sin( r) / reduced_cos( r);
  1133.         when 3| 7 =>
  1134.            if r /= 0.0 then
  1135.           return - reduced_cos( r) / reduced_sin( r);
  1136.            else            -- x must have been huge
  1137.           return 0.0;
  1138.            end if;
  1139.      end case;
  1140.       else
  1141.      raise argument_error;
  1142.       end if;
  1143.    end cot;
  1144.    
  1145.    
  1146.    function cot (x, cycle : Float) return Float is
  1147. --
  1148. -- Declaration:
  1149. --   function COT (X, CYCLE : FLOAT_TYPE) return FLOAT_TYPE;
  1150. -- Definition:
  1151. --   COT(X,CYCLE) ~ cot(2Pi*X/CYCLE)
  1152. -- Usage:
  1153. --   Z := COT(X, 360.0);   -- X in degrees
  1154. --   Z := COT(X, 1.0);     -- X in bams
  1155. --   Z := COT(X, CYCLE);   -- X in units such that one complete
  1156. --                         -- cycle of rotation corresponds to
  1157. --                         -- X = CYCLE
  1158. -- Domain:
  1159. --   (a) X /= k*CYCLE/2.0, for integer k 
  1160. --   (b) CYCLE > 0.0
  1161. -- Range:
  1162. --   Mathematically unbounded
  1163. -- Accuracy:
  1164. --   (a) Maximum relative error = 4.0 * FLOAT_TYPE'BASE'EPSILON 
  1165. --   (b) COT(X,CYCLE) = 0.0 when X=(2k+1)*CYCLE/4.0, for integer k
  1166. --
  1167.       r: Float;
  1168.       oct: octant_number;
  1169.  
  1170.    begin
  1171.       if cycle > 0.0 then
  1172.      period_reduce( x, cycle, r, oct);
  1173.      r := r / cycle * two_pi;    -- underflow denormal allowed
  1174.      case oct is
  1175.         when 0| 4 =>
  1176.            if r /= 0.0 then
  1177.           return reduced_cos( r) / reduced_sin( r);
  1178.            else
  1179.           raise argument_error;
  1180.            end if;
  1181.         when 1| 5 =>
  1182.            return reduced_sin( r) / reduced_cos( r);
  1183.         when 2| 6 =>
  1184.            return - reduced_sin( r) / reduced_cos( r);
  1185.         when 3| 7 =>
  1186.            if r /= 0.0 then
  1187.           return - reduced_cos( r) / reduced_sin( r);
  1188.            else
  1189.           raise argument_error;
  1190.            end if;
  1191.      end case;
  1192.       else
  1193.      raise argument_error;
  1194.       end if;
  1195.    end cot;
  1196.    
  1197.    
  1198.    function arcsin (x : Float) return Float is
  1199. --
  1200. -- Declaration:
  1201. --   function ARCSIN (X : FLOAT_TYPE) return FLOAT_TYPE;
  1202. -- Definition:
  1203. --   ARCSIN(X) ~ arcsine(X)
  1204. -- Usage:
  1205. --   Z := ARCSIN(X);   -- Z in radians
  1206. -- Domain:
  1207. --   |X| <= 1.0
  1208. -- Range:
  1209. --   |ARCSIN(X)| <= Pi/2
  1210. -- Accuracy:
  1211. --   (a) Maximum relative error = 4.0 * FLOAT_TYPE'BASE'EPSILON 
  1212. --   (b) ARCSIN(0.0) = 0.0
  1213. --   (c) ARCSIN(1.0) = Pi/2
  1214. --   (d) ARCSIN(-1.0) = -Pi/2
  1215. -- Notes:
  1216. --   - Pi/2 and Pi/2 are not safe numbers of FLOAT_TYPE. Accordingly,
  1217. --   an implementation may exceed the range limits, but only slightly;
  1218. --   cf. Section 9 for a precise statement of the requirements. Similarly,
  1219. --   when accuracy requirement (c) or (d) applies, an implementation may
  1220. --   approximate the prescribed result, but only within narrow limits;
  1221. --   cf. Section 10 for a precise statement of the requirements.
  1222. --
  1223.    begin
  1224.       if abs( x) <= 1.0 then
  1225.            return arctan( y => x, x => complement( x));
  1226.       else
  1227.            raise argument_error;
  1228.       end if;
  1229.    end arcsin;
  1230.    
  1231.    
  1232.    function arcsin (x, cycle : Float) return Float is
  1233. --
  1234. -- Declaration:
  1235. --   function ARCSIN (X, CYCLE : FLOAT_TYPE) return FLOAT_TYPE;
  1236. -- Definition:
  1237. --   ARCSIN(X,CYCLE) ~ arcsin(X)*CYCLE/2Pi
  1238. -- Usage:
  1239. --   Z := ARCSIN(X, 360.0);   -- Z in degrees
  1240. --   Z := ARCSIN(X, 1.0);     -- Z in bams
  1241. --   Z := ARCSIN(X, CYCLE);   -- Z in units such that one complete
  1242. --                            -- cycle of rotation corresponds to
  1243. --                            -- Z = CYCLE
  1244. -- Domain:
  1245. --   (a) |X| <= 1.0
  1246. --   (b) CYCLE > 0.0
  1247. -- Range:
  1248. --   |ARCSIN(X,CYCLE) <= CYCLE/4.0
  1249. -- Accuracy:
  1250. --   (a) Maximum relative error = 4.0 * FLOAT_TYPE'BASE'EPSILON 
  1251. --   (b) ARCSIN(0.0,CYCLE) = 0.0
  1252. --   (c) ARCSIN(1.0,CYCLE) = CYCLE/4.0
  1253. --   (d) ARCSIN(-1.0,CYCLE) = -CYCLE/4.0
  1254. -- Notes:
  1255. --   - CYCLE/4.0 and CYCLE/4.0 
  1256. --   might not be safe numbers of FLOAT_TYPE. Accordingly,
  1257. --   an implementation may exceed the range limits, but only slightly;
  1258. --   cf. Section 9 for a precise statement of the requirements. Similarly,
  1259. --   when accuracy requirement (c) or (d) applies, an implementation may
  1260. --   approximate the prescribed result, but only within narrow limits;
  1261. --   cf. Section 10 for a precise statement of the requirements.
  1262. --
  1263.    begin
  1264.       if abs( x) <= 1.0 and cycle > 0.0 then
  1265.            return arctan( y => x,
  1266.                    x => complement( x),
  1267.             cycle => cycle);
  1268.       else
  1269.            raise argument_error;
  1270.       end if;
  1271.    end arcsin;
  1272.    
  1273.    
  1274.    function arccos (x : Float) return Float is
  1275. --
  1276. -- Declaration:
  1277. --   function ARCCOS (X : FLOAT_TYPE) return FLOAT_TYPE;
  1278. -- Definition:
  1279. --   ARCCOS(X) ~ arccos(X)
  1280. -- Usage:
  1281. --   Z := ARCCOS(X);   -- Z in radians
  1282. -- Domain:
  1283. --   |X| <= 1.0
  1284. -- Range:
  1285. --   0.0 <= ARCCOS(X) <= Pi
  1286. -- Accuracy:
  1287. --   (a)  Maximum relative error = 4.0 * FLOAT_TYPE'BASE'EPSILON 
  1288. --   (b) ARCCOS(1.0) = 0.0
  1289. --   (c) ARCCOS(0.0) = Pi/2
  1290. --   (d) ARCCOS(-1.0) = Pi
  1291. -- Notes:
  1292. --   Pi/2 and Pi are not safe numbers of FLOAT_TYPE. Accordingly,
  1293. --   an implementation may exceed the range limits, but only slightly;
  1294. --   cf. Section 9 for a precise statement of the requirements. Similarly,
  1295. --   when accuracy requirement (c) or (d) applies, an implementation may
  1296. --   approximate the prescribed result, but only within narrow limits;
  1297. --   cf. Section 10 for a precise statement of the requirements.
  1298. --
  1299.    begin
  1300.       if abs( x) <= 1.0 then
  1301.            return arctan( y => complement( x), x => x);
  1302.       else
  1303.            raise argument_error;
  1304.       end if;
  1305.    end arccos;
  1306.    
  1307.    
  1308.    function arccos (x, cycle : Float) return Float is
  1309. --
  1310. -- Declaration:
  1311. --   function ARCCOS (X, CYCLE : FLOAT_TYPE) return FLOAT_TYPE;
  1312. -- Definition:
  1313. --   ARCCOS(X,CYCLE) ~ arccos(X)*CYCLE/2Pi
  1314. -- Usage:
  1315. --   Z := ARCCOS(X, 360.0);   -- Z in degrees
  1316. --   Z := ARCCOS(X, 1.0);     -- Z in bams
  1317. --   Z := ARCCOS(X, CYCLE);   -- Z in units such that one complete
  1318. --                            -- cycle of rotation corresponds to
  1319. --                            -- Z = CYCLE
  1320. -- Domain:
  1321. --   (a) |X| <= 1.0
  1322. --   (b) CYCLE > 0.0
  1323. -- Range:
  1324. --   0.0 <= ARCCOS(X,CYCLE) <= CYCLE/2.0
  1325. -- Accuracy:
  1326. --   (a) Maximum relative error = 4.0 * FLOAT_TYPE'BASE'EPSILON 
  1327. --   (b) ARCCOS(1.0,CYCLE) = 0.0
  1328. --   (c) ARCCOS(0.0,CYCLE) = CYCLE/4.0
  1329. --   (d) ARCCOS(-1.0,CYCLE) = CYCLE/2.0
  1330. -- Notes:
  1331. --   CYCLE/4.0 and CYCLE/2.0 
  1332. --   might not be safe numbers of FLOAT_TYPE. Accordingly,
  1333. --   an implementation may exceed the range limits, but only slightly;
  1334. --   cf. Section 9 for a precise statement of the requirements. Similarly,
  1335. --   when accuracy requirement (c) or (d) applies, an implementation may
  1336. --   approximate the prescribed result, but only within narrow limits;
  1337. --   cf. Section 10 for a precise statement of the requirements.
  1338. --
  1339.    begin
  1340.       if abs( x) <= 1.0 and cycle > 0.0 then
  1341.            return arctan( y => complement( x),
  1342.                         x => x,
  1343.             cycle => cycle);
  1344.       else
  1345.            raise argument_error;
  1346.       end if;
  1347.    end arccos;
  1348.    
  1349.  
  1350.    function arctan (y : Float;
  1351.                          x : Float := 1.0) return Float is
  1352. --
  1353. -- Declaration:
  1354. --   function ARCTAN (Y : FLOAT_TYPE;
  1355. --                    X : FLOAT_TYPE := 1.0) return FLOAT_TYPE;
  1356. -- Definition:
  1357. --   (a) ARCTAN(Y) ~ arctan(Y)
  1358. --   (b) ARCTAN(Y,X) ~ arctan(Y/X) when X >= 0.0
  1359. --                     arctan(Y/X)+Pi when X < 0.0 and Y >= 0.0
  1360. --                     arctan(Y/X)-Pi when X < 0.0 and Y < 0.0
  1361. -- Usage:
  1362. --   Z := ARCTAN(Y);      -- Z, in radians, is the angle (in the
  1363. --                        -- quadrant containing the point (1.0,Y))
  1364. --                        -- whose tangent is Y
  1365. --   Z := ARCTAN(Y, X);   -- Z, in radians, is the angle (in the
  1366. --                        -- quadrant containing the point (X,Y))
  1367. --                        -- whose tangent is Y/X
  1368. -- Domain:
  1369. --   X /= 0.0 when Y = 0.0
  1370. -- Range:
  1371. --     (a) |ARCTAN(Y)| <= Pi/2
  1372. --     (b) 0.0 < ARCTAN(Y,X) <= Pi when Y >= 0.0
  1373. --     (c) -Pi <= ARCTAN(Y,X) <= 0.0 when Y < 0.0
  1374. -- Accuracy:
  1375. --   (a) Maximum relative error = 4.0 * FLOAT_TYPE'BASE'EPSILON 
  1376. --   (b) ARCTAN(0.0) = 0.0
  1377. --   (c) ARCTAN((0.0,X) = 0.0 when X > 0.0
  1378. --                        Pi when X < 0.0
  1379. --   (d) ARCTAN(Y,0.0) = Pi/2 when Y > 0.0
  1380. --                       -Pi/2 when Y < 0.0
  1381. -- Notes:
  1382. --   -Pi,-Pi/2,Pi/2 and Pi are not safe numbers of FLOAT_TYPE. Accordingly,
  1383. --   an implementation may exceed the range limits, but only slightly;
  1384. --   cf. Section 9 for a precise statement of the requirements. Similarly,
  1385. --   when accuracy requirement (c) or (d) applies, an implementation may
  1386. --   approximate the prescribed result, but only within narrow limits;
  1387. --   cf. Section 10 for a precise statement of the requirements.
  1388. --
  1389.    begin
  1390.       if x = 0.0 and y = 0.0 then
  1391.      raise argument_error;
  1392.       elsif x >= y then
  1393.      if x >= - y then
  1394.         return reduced_arctan( y / x); -- near +x axis
  1395.      else
  1396.         return - half_pi - reduced_arctan( x / y); -- near -y axis
  1397.      end if;
  1398.       elsif x >= - y then
  1399.      return half_pi - reduced_arctan( x / y); -- near +y axis
  1400.       elsif y >= 0.0 then
  1401.      return pi + reduced_arctan( y / x); -- just above -x axis
  1402.       else
  1403.      return - pi + reduced_arctan( y / x); -- just below -x axis
  1404.       end if;
  1405.    end arctan;
  1406.    
  1407.    
  1408.    function arctan (y     : Float;
  1409.                          x     : Float := 1.0;
  1410.             cycle : Float) return Float is
  1411. --
  1412. -- Declaration:
  1413. --   function ARCTAN (Y     : FLOAT_TYPE;
  1414. --                    X     : FLOAT_TYPE := 1.0;
  1415. --                    CYCLE : FLOAT_TYPE)        return FLOAT_TYPE;
  1416. -- Definition:
  1417. --   (a) ARCTAN(Y,CYCLE) ~ arctan(Y)*CYCLE/2Pi
  1418. --   (b) ARCTAN(Y,X,CYCLE) ~ arctan(Y/X)*CYCLE/2Pi when X >= 0.0
  1419. --                           (arctan(Y/X)+Pi)*CYCLE/2Pi when X < 0.0 and Y >= 0.0
  1420. --                           (arctan(Y/X)-Pi)*CYCLE/2Pi when X < 0.0 and Y < 0.0
  1421. -- Usage:       
  1422. --   Z := ARCTAN(Y, CYCLE => 360.0);   -- Z, in degrees, is the
  1423. --                                     -- angle (in the quadrant
  1424. --                                     -- containing the point
  1425. --                                     -- (1.0,Y)) whose tangent is Y
  1426. --   Z := ARCTAN(Y, CYCLE => 1.0);     -- Z, in bams, is the
  1427. --                                     -- angle (in the quadrant
  1428. --                                     -- containing the point
  1429. --                                     -- (1.0,Y)) whose tangent is Y
  1430. --   Z := ARCTAN(Y, CYCLE => CYCLE);   -- Z, in units such that one
  1431. --                                     -- complete cycle of rotation
  1432. --                                     -- corresponds to Z = CYCLE,
  1433. --                                     -- is the angle (in the
  1434. --                                     -- quadrant containing the
  1435. --                                     -- point (1.0,Y)) whose
  1436. --                                     -- tangent is Y
  1437. --   Z := ARCTAN(Y, X, 360.0);         -- Z, in degrees, is the
  1438. --                                     -- angle (in the quadrant
  1439. --                                     -- containing the point (X,Y))
  1440. --                                     -- whose tangent is Y/X
  1441. --   Z := ARCTAN(Y, X, 1.0);           -- Z, in bams, is the
  1442. --                                     -- angle (in the quadrant
  1443. --                                     -- containing the point (X,Y))
  1444. --                                     -- whose tangent is Y/X
  1445. --   Z := ARCTAN(Y, X, CYCLE);         -- Z, in units such that one
  1446. --                                     -- complete cycle of rotation
  1447. --                                     -- corresponds to Z = CYCLE,
  1448. --                                     -- is the angle (in the
  1449. --                                     -- quadrant containing the
  1450. --                                     -- point (X,Y)) whose
  1451. --                                     -- tangent is Y/X
  1452. -- Domain:
  1453. --   (a) X /= 0.0 when Y = 0.0
  1454. --   (b) CYCLE > 0.0
  1455. -- Range:
  1456. --   (a) |ARCTAN(Y,CYCLE)| <= CYCLE/4.0
  1457. --   (b) 0.0 <= ARCTAN(Y,X,CYCLE) <= CYCLE/2.0 when Y >= 0.0
  1458. --   (c) -CYCLE/2.0 <= ARCTAN(Y,X,CYCLE) <= 0.0 when Y < 0.0
  1459. -- Accuracy:
  1460. --   (a) Maximum relative error = 4.0 * FLOAT_TYPE'BASE'EPSILON 
  1461. --   (b) ARCTAN(0.0,CYCLE) = 0.0
  1462. --   (c) ARCTAN(0.0,X,CYCLE) = 0.0 when X > 0.0
  1463. --                             CYCLE/2.0 when X < 0.0
  1464. --   (d) ARCTAN(Y,0.0,CYCLE) = CYCLE/4.0 when Y > 0.0
  1465. --                             -CYCLE/4.0 when Y < 0.0
  1466. -- Notes:
  1467. --   -CYCLE/2.0, -CYCLE/4.0, CYCLE/4.0 and CYCLE/2.0 
  1468. --   might not be safe numbers of FLOAT_TYPE. Accordingly,
  1469. --   an implementation may exceed the range limits, but only slightly;
  1470. --   cf. Section 9 for a precise statement of the requirements. Similarly,
  1471. --   when accuracy requirement (c) or (d) applies, an implementation may
  1472. --   approximate the prescribed result, but only within narrow limits;
  1473. --   cf. Section 10 for a precise statement of the requirements.
  1474. --
  1475.       scale: constant Float := cycle / two_pi;
  1476.  
  1477.    begin
  1478.       -- what happens if scale underflows?  cycle/4 underflows?
  1479.       if cycle <= 0.0 or else (x = 0.0 and y = 0.0) then
  1480.      raise argument_error;
  1481.       elsif x >= y then
  1482.      if x >= - y then
  1483.         -- near +x axis
  1484.         return scale * reduced_arctan( y / x);
  1485.      else
  1486.         -- near -y axis
  1487.         return - cycle / 4.0 - scale * reduced_arctan( x / y);
  1488.      end if;
  1489.       elsif x >= - y then
  1490.      -- near +y axis
  1491.      return cycle / 4.0 - scale * reduced_arctan( x / y);
  1492.       elsif y >= 0.0 then
  1493.      -- just above -x axis
  1494.      return cycle / 2.0 + scale * reduced_arctan( y / x);
  1495.       else
  1496.      -- just below -x axis
  1497.      return - cycle / 2.0 + scale * reduced_arctan( y / x);
  1498.       end if;
  1499.    end arctan;
  1500.    
  1501.    
  1502.    function arccot (x : Float;
  1503.                          y : Float := 1.0) return Float is
  1504. --
  1505. -- Declaration:
  1506. --   function ARCCOT (X : FLOAT_TYPE;
  1507. --                    Y : FLOAT_TYPE := 1.0) return FLOAT_TYPE;
  1508. -- Definition:
  1509. --   (a) ARCCOT(X) ~ arccot(X)
  1510. --   (b) ARCCOT(X,Y) ~ arccot(X/Y) when Y >= 0.0
  1511. --                     arccot(X/Y)-Pi when Y < 0.0
  1512. -- Usage:
  1513. --   Z := ARCCOT(X);      -- Z, in radians, is the angle (in the
  1514. --                        -- quadrant containing the point (X,1.0)
  1515. --                        -- whose cotangent is X
  1516. --   Z := ARCCOT(X, Y);   -- Z, in radians, is the angle (in the
  1517. --                        -- quadrant containing the point (X,Y))
  1518. --                        -- whose cotangent is X/Y
  1519. -- Domain:
  1520. --   Y /= 0.0 when X = 0.0
  1521. -- Range:
  1522. --   (a) 0.0 <= ARCCOT(X) <= Pi
  1523. --   (b) 0.0 <= ARCCOT(X,Y) <= Pi when Y >= 0.0
  1524. --   (c) -Pi <= ARCCOT(X,Y) <= 0.0 when Y < 0.0
  1525. -- Accuracy:
  1526. --   (a) Maximum relative error = 4.0 * FLOAT_TYPE'BASE'EPSILON 
  1527. --   (b) ARCCOT(0.0) = Pi/2
  1528. --   (c) ARCCOT(0.0,Y) = Pi/2 when Y > 0.0
  1529. --                       -Pi/2 when Y < 0.0
  1530. --   (d) ARCCOT(X,0.0) = 0.0 when X > 0.0
  1531. --                       Pi when X < 0.0
  1532. -- Notes:
  1533. --   -Pi,-Pi/2,Pi/2 and Pi are not safe numbers of FLOAT_TYPE. Accordingly,
  1534. --   an implementation may exceed the range limits, but only slightly;
  1535. --   cf. Section 9 for a precise statement of the requirements. Similarly,
  1536. --   when accuracy requirement (b), (c), or (d) applies, an implementation may
  1537. --   approximate the prescribed result, but only within narrow limits;
  1538. --   cf. Section 10 for a precise statement of the requirements.
  1539. --
  1540.    begin
  1541.       return arctan( y => y, x => x);
  1542.    end arccot;
  1543.    
  1544.    
  1545.    function arccot (x     : Float;
  1546.                          y     : Float := 1.0;
  1547.             cycle : Float) return Float is
  1548. --
  1549. -- Declaration:
  1550. --   function ARCCOT (X     : FLOAT_TYPE;
  1551. --                    Y     : FLOAT_TYPE := 1.0;
  1552. --                    CYCLE : FLOAT_TYPE)        return FLOAT_TYPE;
  1553. -- Definition:
  1554. --   (a) ARCCOT(X,CYCLE) ~ arccot(X)*CYCLE/2Pi
  1555. --   (b) ARCCOT(X,Y) ~ arccot(X/Y)*CYCLE/2Pi when Y >= 0.0
  1556. --                     (arccot(X/Y)-Pi)*CYCLE/2Pi
  1557. -- Usage:
  1558. --   Z := ARCCOT(X, CYCLE => 360.0);   -- Z, in degrees, is the
  1559. --                                     -- angle (in the quadrant
  1560. --                                     -- containing the point
  1561. --                                     -- (X,1.0)) whose cotangent
  1562. --                                     -- is X
  1563. --   Z := ARCCOT(X, CYCLE => 1.0);     -- Z, in bams, is the
  1564. --                                     -- angle (in the quadrant
  1565. --                                     -- containing the point
  1566. --                                     -- (X,1.0)) whose cotangent
  1567. --                                     -- is X
  1568. --   Z := ARCCOT(X, CYCLE => CYCLE);   -- Z, in units such that one
  1569. --                                     -- complete cycle of rotation
  1570. --                                     -- corresponds to Z = CYCLE,
  1571. --                                     -- is the angle (in the
  1572. --                                     -- quadrant containing the
  1573. --                                     -- point (X,1.0)) whose
  1574. --                                     -- cotangent is X
  1575. --   Z := ARCCOT(X, Y, 360.0);         -- Z, in degrees, is the
  1576. --                                     -- angle (in the quadrant
  1577. --                                     -- containing the point (X,Y))
  1578. --                                     -- whose cotangent is X/Y
  1579. --   Z := ARCCOT(X, Y, 1.0);           -- Z, in bams, is the
  1580. --                                     -- angle (in the quadrant
  1581. --                                     -- containing the point (X,Y)
  1582. --                                     -- whose cotangent is X/Y
  1583. --   Z := ARCCOT(X, Y, CYCLE);         -- Z, in units such that one
  1584. --                                     -- complete cycle of rotation
  1585. --                                     -- corresponds to Z = CYCLE
  1586. --                                     -- is the angle (in the
  1587. --                                     -- quadrant containing the
  1588. --                                     -- point (X,Y)) whose
  1589. --                                     -- cotangent is X/Y
  1590. -- Domain:
  1591. --   (a) Y /= 0.0 when X = 0.0
  1592. --   (b) CYCLE > 0.0
  1593. -- Range:
  1594. --   (a) 0.0 <= ARCCOT((X,CYCLE) <= CYCLE/2.0
  1595. --   (b) 0.0 <= ARCCOT(X,Y,CYCLE) <= CYCLE/2.0 when Y >= 0.0
  1596. --   (c) -CYCLE/2.0 <= ARCCOT(X,Y,CYCLE) <= 0.0 when Y < 0.0
  1597. -- Accuracy:
  1598. --   (a) Maximum relative error = 4.0 * FLOAT_TYPE'BASE'EPSILON 
  1599. --   (b) ARCCOT(0.0,CYCLE) = CYCLE/4.0
  1600. --   (c) ARCCOT(0.0,Y,CYCLE) = CYCLE/4.0 when Y > 0.0
  1601. --                             -CYCLE/4.0 when Y < 0.0
  1602. --   (d) ARCCOT(X,0.0,CYCLE) = 0.0 when X > 0.0
  1603. --                             CYCLE/2.0 when X < 0.0
  1604. -- Notes:
  1605. --   - CYCLE/2.0, - CYCLE/4.0, CYCLE/4.0 and CYCLE/2.0
  1606. --   might not be safe numbers of FLOAT_TYPE. Accordingly,
  1607. --   an implementation may exceed the range limits, but only slightly;
  1608. --   cf. Section 9 for a precise statement of the requirements. Similarly,
  1609. --   when accuracy requirement (b), (c), or (d) applies, an implementation may
  1610. --   approximate the prescribed result, but only within narrow limits;
  1611. --   cf. Section 10 for a precise statement of the requirements.
  1612. --
  1613.    begin
  1614.       return arctan( y => y, x => x, cycle => cycle);
  1615.    end arccot;
  1616.    
  1617.    
  1618.    function sinh (x : Float) return Float is
  1619. --
  1620. -- Declaration:
  1621. --   function SINH (X : FLOAT_TYPE) return FLOAT_TYPE;
  1622. -- Definition:
  1623. --   SINH(X) ~ sinh X
  1624. -- Usage:
  1625. --   Z := SINH(X);
  1626. -- Domain:
  1627. --   Mathematically unbounded
  1628. -- Range:
  1629. --   Mathematically unbounded
  1630. -- Accuracy:
  1631. --   (a) Maximum relative error = 8.0 * FLOAT_TYPE'BASE'EPSILON 
  1632. --   (b) SINH(0.0) = 0.0
  1633. -- Notes:
  1634. --   The usable domain of SINH is approximately given by
  1635. --   |X| <= ln(FLOAT_TYPE'SAFE_LARGE)+ln(2.0)
  1636. --
  1637.       ax: constant Float := abs( x);
  1638.       meet_exp: constant := (mant + 1) * nat_log_2 / 2.0;
  1639.  
  1640.    begin
  1641.       if ax > meet_exp then
  1642.      -- essentially +-exp(ax)/2, but must round and overflow properly
  1643.      declare
  1644.         ex: constant Float := (exp_1 / 2.0) * exp( ax - 1.0);
  1645.      begin
  1646.         if x >= 0.0 then
  1647.            return ex;
  1648.         else
  1649.            return - ex;
  1650.         end if;
  1651.      end; 
  1652.       elsif abs( x) <= 0.25 then
  1653.            -- small args need an economized taylor series, relerr<1.5e-9.
  1654.      declare
  1655.         a1: constant := 1.00000000151607e+00;
  1656.         a3: constant := 1.66666230112174e-01;
  1657.         a5: constant := 8.35195337329340e-03;
  1658.         x2: constant Float := x * x;
  1659.      begin
  1660.         return x * (a1 + x2 * (a3 + x2 *  a5));
  1661.      end;
  1662.       else                -- 1/4 < abs x <= meet_exp
  1663.            declare
  1664.         ex: constant Float := exp( x);
  1665.      begin
  1666.         return 0.5 * (ex - 1.0 / ex);
  1667.      end;
  1668.       end if;
  1669.    end sinh;
  1670.    
  1671.    
  1672.    function cosh (x : Float) return Float is
  1673. --
  1674. -- Declaration:
  1675. --   function COSH (X : FLOAT_TYPE) return FLOAT_TYPE;
  1676. -- Definition:
  1677. --   COSH(X) ~ cosh X
  1678. -- Usage:
  1679. --   Z := COSH(X);
  1680. -- Domain:
  1681. --   Mathematically unbounded
  1682. -- Range:
  1683. --   COSH(X) >= 1.0
  1684. -- Accuracy:
  1685. --   (a) Maximum relative error = 8.0 * FLOAT_TYPE'BASE'EPSILON 
  1686. --   (b) COSH(0.0) = 1.0
  1687. -- Notes:
  1688. --   The usable domain of COSH is approximately given by
  1689. --   |X| <= ln(FLOAT_TYPE'SAFE_LARGE)+ln(2.0)
  1690. --
  1691.       ax: constant Float := abs( x);
  1692.       meet_exp: constant := (mant + 1) * nat_log_2 / 2.0;
  1693.  
  1694.    begin
  1695.       if ax > meet_exp then
  1696.      -- essentially exp(ax)/2, but must round and overflow properly
  1697.      return (exp_1 / 2.0) * exp( ax - 1.0);
  1698.       else
  1699.      declare
  1700.               ex: constant Float := exp( x); -- moderate size
  1701.      begin
  1702.         return 0.5 * (ex + 1.0 / ex);
  1703.      end;
  1704.       end if;
  1705.    end cosh;
  1706.    
  1707.    
  1708.    function tanh (x : Float) return Float is
  1709. --
  1710. -- Declaration:
  1711. --   function TANH (X : FLOAT_TYPE) return FLOAT_TYPE;
  1712. -- Definition:
  1713. --   TANH(X) ~ tanh X
  1714. -- Usage:
  1715. --   Z := TANH(X);
  1716. -- Domain:
  1717. --   Mathematically unbounded
  1718. -- Range:
  1719. --   |TANH(X)| <= 1.0
  1720. -- Accuracy:
  1721. --   (a) Maximum relative error = 8.0 * FLOAT_TYPE'BASE'EPSILON 
  1722. --   (b) TANH(0.0) = 0.0
  1723. --
  1724.       ax: constant Float := abs( x);
  1725.       meet_exp: constant := (mant + 1) * nat_log_2 / 2.0;
  1726.  
  1727.    begin
  1728.       -- check whether exp(-ax) negligible wrt exp(ax)
  1729.       if ax > meet_exp then
  1730.      if x >= 0.0 then
  1731.         return 1.0;
  1732.      else
  1733.         return -1.0;
  1734.      end if;
  1735.       elsif ax <= 0.25 then
  1736.            -- small args need an economized taylor series, relerr<2.5e-9.
  1737.      declare
  1738.         x2: constant Float := x * x;
  1739.         a0: constant :=  9.99999997535969e-01;
  1740.         a2: constant := -3.33332067277568e-01;
  1741.         a4: constant :=  1.33231448939255e-01;
  1742.         a6: constant := -5.13290121422199e-02;
  1743.      begin
  1744.         return x * (a0 + x2 * (a2
  1745.                + x2 * (a4
  1746.                + x2 *  a6)));
  1747.      end;
  1748.       else                -- 1/4 < abs x <= meet_exp
  1749.      declare
  1750.               ex: Float := exp( x);
  1751.      begin
  1752.         return (ex - 1.0 / ex) / (ex + 1.0 / ex);
  1753.      end;
  1754.       end if;
  1755.    end tanh;
  1756.    
  1757.    
  1758.    function coth (x : Float) return Float is
  1759. --
  1760. -- Declaration:
  1761. --   function COTH (X : FLOAT_TYPE) return FLOAT_TYPE;
  1762. -- Definition:
  1763. --   COTH(X) ~ coth X
  1764. -- Usage:
  1765. --   Z := COTH(X);
  1766. -- Domain:
  1767. --   X /= 0.0
  1768. -- Range:
  1769. --   |COTH(X)| >= 1.0
  1770. -- Accuracy:
  1771. --   Maximum relative error = 8.0 * FLOAT_TYPE'BASE'EPSILON 
  1772. --
  1773.    begin
  1774.       if x /= 0.0 then
  1775.      return 1.0 / tanh( x);
  1776.       else
  1777.      raise argument_error;
  1778.       end if;
  1779.    end coth;
  1780.    
  1781.    
  1782.    function arcsinh (x : Float) return Float is
  1783. --
  1784. -- Declaration:
  1785. --   function ARCSINH (X : FLOAT_TYPE) return FLOAT_TYPE;
  1786. -- Definition:
  1787. --   ARCSINH(X) ~ arcsinh X
  1788. -- Usage:
  1789. --   Z := ARCSINH(X);
  1790. -- Domain:
  1791. --   Mathematically unbounded
  1792. -- Range:
  1793. --   Mathematically unbounded
  1794. -- Accuracy:
  1795. --   (a) Maximum relative error = 8.0 * FLOAT_TYPE'BASE'EPSILON 
  1796. --   (b) ARCSINH(0.0) = 0.0
  1797. -- Notes:
  1798. --   The reachable range of ARCSINH is approximately given by
  1799. --     |ARCSINH(X)| <= ln(FLOAT_TYPE'SAFE_LARGE)+ln(2.0)
  1800. --
  1801.       ax: Float := abs( x);
  1802.       as: Float;
  1803.       reduce_edge: constant := sqrt_2 / 4.0;
  1804.       -- see reduced_logdel_2 for its domain definition
  1805.  
  1806.    begin
  1807.       if ax >= recip_sqrt_epsilon then
  1808.      as := 1.0 + log_2( ax);
  1809.       elsif ax >= reduce_edge then
  1810.      as := log_2( ax + sqrt( 1.0 + sqr(ax)));
  1811.       elsif ax > Float'base'epsilon then
  1812.                     -- ax < sqrt(1/8)
  1813.      --       z     = x+sqrt(1+x**2)
  1814.      -- (z-1)/(z+1) = x/(1+sqrt(1+x**2))
  1815.      as := reduced_logdel_2( ax / (1.0 + sqrt( 1.0 + sqr( ax))));
  1816.       else
  1817.      return x;            -- sign already correct for small args
  1818.       end if;
  1819.       if x >= 0.0 then
  1820.            return nat_log_2 * as;
  1821.       else
  1822.            return - nat_log_2 * as;
  1823.       end if;
  1824.    end arcsinh;
  1825.    
  1826.    
  1827.    function arccosh (x : Float) return Float is
  1828. --
  1829. -- Declaration:
  1830. --   function ARCCOSH (X : FLOAT_TYPE) return FLOAT_TYPE;
  1831. -- Definition:
  1832. --   ARCCOSH(X) ~ arccosh X
  1833. -- Usage:
  1834. --   Z := ARCCOSH(X);
  1835. -- Domain:
  1836. --   X >= 1.0
  1837. -- Range:
  1838. --   ARCCOSH(X) >= 0.0
  1839. -- Accuracy:
  1840. --   (a) Maximum relative error = 8.0 * FLOAT_TYPE'BASE'EPSILON 
  1841. --   (b) ARCCOSH(1.0) = 0.0
  1842. -- Notes:
  1843. --   The upper bound of the reachable range of ARCCOSH is approximately given
  1844. --   by  ARCCOSH(X) <= ln(FLOAT_TYPE'SAFE_LARGE)+ln(2.0)
  1845. --
  1846.       reduce_edge: constant := 0.75 * sqrt_2;
  1847.       -- see reduced_logdel_2 for its domain definition
  1848.  
  1849.    begin
  1850.       if x >= recip_sqrt_epsilon then
  1851.      return nat_log_2 * (1.0 + log_2( x));
  1852.       elsif x >= reduce_edge then
  1853.      return nat_log_2 * log_2( x + sqrt( sqr( x) - 1.0));
  1854.       elsif x >= 1.0 then
  1855.      --    z        = x+sqrt(x**2-1)
  1856.      -- (z-1)/(z+1) = sqrt((x-1)/(x+1))
  1857.      return nat_log_2 * reduced_logdel_2( sqrt( (x - 1.0) / (x + 1.0) ));
  1858.       else
  1859.      raise argument_error;
  1860.       end if;
  1861.    end arccosh;
  1862.    
  1863.    
  1864.    function arctanh (x : Float) return Float is
  1865. --
  1866. -- Declaration:
  1867. --   function ARCTANH (X : FLOAT_TYPE) return FLOAT_TYPE;
  1868. -- Definition:
  1869. --   ARCTANH(X) ~ arctanh X
  1870. -- Usage:
  1871. --   Z := ARCTANH(X);
  1872. -- Domain:
  1873. --   |X| < 1.0
  1874. -- Range:
  1875. --   Mathematically unbounded
  1876. -- Accuracy:
  1877. --   (a) Maximum relative error = 8.0 * FLOAT_TYPE'BASE'EPSILON 
  1878. --   (b) ARCTANH(0.0) = 0.0
  1879. --
  1880.       ax: constant Float := abs( x);
  1881.       half_ln_2: constant := nat_log_2 / 2.0;
  1882.       reduce_edge: constant := 1.0 / (3.0 + 2.0 * sqrt_2);
  1883.       -- see reduced_logdel_2 for its domain definition
  1884.  
  1885.    begin
  1886.       if ax >= 1.0 then
  1887.      raise argument_error;
  1888.       elsif ax < reduce_edge then
  1889.      return reduced_logdel_2( x) * half_ln_2;
  1890.       else
  1891.      return log_2( (1.0 + x) / (1.0 - x)) * half_ln_2;
  1892.       end if;
  1893.    end arctanh;
  1894.    
  1895.    
  1896.    function arccoth (x : Float) return Float is
  1897. --
  1898. -- Declaration:
  1899. --   function ARCCOTH (X : FLOAT_TYPE) return FLOAT_TYPE;
  1900. -- Definition:
  1901. --   ARCCOTH(X) ~ arccoth X
  1902. -- Usage:
  1903. --   Z := ARCCOTH(X);
  1904. -- Domain:
  1905. --   |X| > 1.0
  1906. -- Range:
  1907. --   Mathematically unbounded
  1908. -- Accuracy:
  1909. --   Maximum relative error = 8.0 * FLOAT_TYPE'BASE'EPSILON 
  1910. --
  1911.       ax: constant Float := abs( x);
  1912.       half_ln_2: constant := nat_log_2 / 2.0;
  1913.       reduce_edge: constant := 3.0 + 2.0 * sqrt_2;
  1914.       -- see reduced_logdel_2 for its domain definition
  1915.  
  1916.    begin
  1917.       if ax <= 1.0 then
  1918.      raise argument_error;
  1919.       elsif ax > reduce_edge then
  1920.      return reduced_logdel_2( 1.0 / x) * half_ln_2;
  1921.       else
  1922.            -- 1 < ax <= 3+2*sqrt2
  1923.      return log_2( (x + 1.0) / (x - 1.0)) * half_ln_2;
  1924.       end if;
  1925.    end arccoth;
  1926.  
  1927. end Math;
  1928.  
  1929. -- $Header: s_elementary_functions_b.ada,v 3.13 90/04/19 12:34:31 broman Rel $
  1930.