The Unsorted BBS Collection/thegreatunsorted.tar/thegreatunsorted/programming/misc_programming/MATH/GEFPC.A

Text File | 1990-08-14 | 42.0 KB | 1,319 lines

-- This body is specifically for Meridian Ada using their MATH_LIB
-- This is not efficient. There are some inaccuracies. Just a sample to
-- get the idea of one method of building the body based on an
-- existing set SQRT, EXP, LN, SIN, COS, and ATAN

with MATH_LIB ;

package body GENERIC_ELEMENTARY_FUNCTIONS is

  PI      : constant := 3.14159_26535_89793_23846_26433_83279_50288_41972 ;
  LOG_TWO : constant := 0.69314_71805_59945_30941_72321_21458_17656_80755 ;
  TWO_PI : constant := 2.0 * PI ;
  HALF_PI : constant := PI / 2.0 ;
  FOURTH_PI : constant := PI / 4.0 ;
  EPSILON : FLOAT_TYPE ;
  SQUARE_ROOT_EPSILON : FLOAT_TYPE ;
  HALF_LOG_EPSILON : FLOAT_TYPE ;
  PREVIOUS : FLOAT_TYPE ;
  LOG_LAST : FLOAT_TYPE ;
  LOG_INVERSE_EPSILON : FLOAT_TYPE ;

  function EXACT_REMAINDER ( X , Y : FLOAT_TYPE ) return FLOAT_TYPE is
    DENOMINATOR : FLOAT_TYPE := abs X ;
    DIVISOR : FLOAT_TYPE := abs Y ;
    REDUCER : FLOAT_TYPE ;
    SIGN : FLOAT_TYPE := 1.0 ;
  begin
    if Y = 0.0 then
      raise CONSTRAINT_ERROR ;
    elsif X = 0.0 then
      return 0.0 ;
    elsif X = Y then
      return 0.0 ;
    elsif DENOMINATOR < DIVISOR then
      return X ;
    end if ;
    while DENOMINATOR >= DIVISOR loop
      -- put divisors mantissa with denominators exponent to make reducer
      REDUCER := DIVISOR ;
      begin
        while REDUCER * 1_048_576.0 < DENOMINATOR loop
          REDUCER := REDUCER * 1_048_576.0 ;
        end loop ;
      exception when others=> null ;
      end;
    begin
      while REDUCER * 1_024.0 < DENOMINATOR loop
        REDUCER := REDUCER * 1_024.0 ;
      end loop ;
    exception when others=> null ;
    end;
    begin
      while REDUCER * 2.0 < DENOMINATOR loop
          REDUCER := REDUCER * 2.0 ;
        end loop ;
      exception
        when others=>
          null ;
      end;
      DENOMINATOR := DENOMINATOR - REDUCER ;
    end loop ;
    if X < 0.0 then
      return - DENOMINATOR ;
    else
      return DENOMINATOR ;
    end if ;
  end EXACT_REMAINDER ;


  function LOCAL_ATAN ( Y : FLOAT_TYPE ;
                        X : FLOAT_TYPE := 1.0 ) return FLOAT_TYPE is
    Z : FLOAT_TYPE ;
    RAW_ATAN : FLOAT_TYPE ;
  begin
    if abs Y > abs X then
      Z := abs ( X / Y ) ;
    else
      Z := abs ( Y / X ) ;
    end if;
    if Z = 0.0 then
      RAW_ATAN := 0.0 ;
    elsif Z = 1.0 then
      RAW_ATAN := FOURTH_PI ;
    elsif Z < SQUARE_ROOT_EPSILON then
      RAW_ATAN := Z ;
    else
      RAW_ATAN := FLOAT_TYPE ( MATH_LIB.ATAN ( FLOAT ( Z ) ) ) ;
    end if ;
    if abs Y > abs X then
      RAW_ATAN := HALF_PI - RAW_ATAN ;
    end if ;
    if X > 0.0 then
      if Y > 0.0 then
        return RAW_ATAN ;
      else -- Y < 0.0
        return -RAW_ATAN ;
      end if ;
    else -- X < 0.0
      if Y > 0.0 then
        return PI - RAW_ATAN ;
      else -- Y < 0.0
        return - ( PI - RAW_ATAN ) ;
      end if ;
    end if ;
  end LOCAL_ATAN ;




-- SQRT
-- Declaration:
--   function SQRT (X : FLOAT_TYPE) return FLOAT_TYPE;
-- Definition:
--   SQRT(X) ~ sqrt(X)
-- Usage:
--   Z := SQRT(X);
-- Domain:
--   X >= 0.0
-- Range:
--   SQRT(X) >= 0.0
-- Accuracy:
--   (a) Maximum relative error = 2.0 * FLOAT_TYPE'BASE'EPSILON
--   (b) SQRT(0.0) = 0.0
-- Notes:
--   (a) The upper bound of the reachable range of SQRT is approximately given
--   by SQRT(X)<=sqrt(FLOAT_TYPE'SAFE_LARGE)
--   (b) Other standards might impose additional constraints on SQRT. For
--   example, the IEEE standards for binary and radix-independent
--   floating-point arithmetic require greater accuracy in the result of SQRT
--   than this standard requires, and they require that SQRT(-0.0)=-0.0.
--   An implementation of SQRT in GENERIC_ELEMENTARY_FUNCTIONS that
--   conforms to this standard will conform to those other standards if it
--   satisfies their additional constraints.

  function SQRT ( X : FLOAT_TYPE ) return FLOAT_TYPE is
  begin
    if X < 0.0 then
      raise ARGUMENT_ERROR ;
    elsif X = 0.0 then
      return X ; -- may be +0.0 or -0.0
    elsif X = 1.0 then
      return 1.0 ; -- needs to be exact for good COMPLEX accuracy
    end if;
    return FLOAT_TYPE ( MATH_LIB.SQRT ( FLOAT ( X ) ) ) ;
  end SQRT ;

-- LOG (natural base)
-- Declaration:
--   function LOG (X : FLOAT_TYPE) return FLOAT_TYPE;
-- Definition:
--   LOG(X) ~ ln(X)
-- Usage:
--   Z := LOG(X);   -- natural logarithm
-- Domain:
--   X > 0.0
-- Range:
--   Mathematically unbounded
-- Accuracy:
--   (a) Maximum relative error = 4.0 * FLOAT_TYPE'BASE'EPSILON
--   (b) LOG(1.0) = 0.0
-- Notes:
--   The reachable range of LOG is approximately given by
--   ln(FLOAT_TYPE'SAFE_SMALL) <= LOG(X) <= ln(FLOAT_TYPE'SAFE_LARGE)

  function LOG ( X : FLOAT_TYPE ) return FLOAT_TYPE is
  begin
    if X <= 0.0 then
      raise ARGUMENT_ERROR ;
    elsif X = 1.0 then
      return 0.0 ;
    end if ;
    return FLOAT_TYPE ( MATH_LIB.LN ( FLOAT ( X ) ) ) ;
  end LOG ;


-- LOG (arbitrary base)
-- Declaration:
--   function LOG (X, BASE : FLOAT_TYPE) return FLOAT_TYPE;
-- Definition:
--   LOG(X,BASE) ~ log to base BASE of X
-- Usage:
--   Z := LOG(X, 10.0);   -- base 10 logarithm
--   Z := LOG(X, 2.0);    -- base 2 logarithm
--   Z := LOG(X, BASE);   -- base BASE logarithm
-- Domain:
--   (a) X > 0.0
--   (b) BASE > 0.0
--   (c) BASE /= 1.0
-- Range:
--   Mathematically unbounded
-- Accuracy:
--   (a) Maximum relative error = 4.0 * FLOAT_TYPE'BASE'EPSILON
--   (b) LOG(1.0,BASE) = 0.0
-- Notes:
--   (a) When BASE > 1.0, the reachable range of LOG is approximately given by
--       log to base BASE of FLOAT_TYPE'SAFE_SMALL <= LOG(X,BASE) <=
--       log to base BASE of FLOAT_TYPE'SAFE_LARGE
--   (b) When 0.0 < BASE < 1.0, the reachable range of LOG is approximately given by
--       log to base BASE of FLOAT_TYPE'SAFE_LARGE <= LOG(X,BASE) <=
--       log to base BASE of FLOAT_TYPE'SAFE_SMALL

  function LOG ( X , BASE : FLOAT_TYPE ) return FLOAT_TYPE is
  begin
    if X <= 0.0 then
      raise ARGUMENT_ERROR ;
    elsif BASE <= 0.0 or else BASE = 1.0 then
      raise ARGUMENT_ERROR ;
    elsif X = 1.0 then
      return 0.0 ;
    end if ;
    return FLOAT_TYPE ( MATH_LIB.LN ( FLOAT ( X ) )
                      / MATH_LIB.LN ( FLOAT ( BASE ) ) ) ;
  end LOG ;


-- EXP
-- Declaration:
--   function EXP (X : FLOAT_TYPE) return FLOAT_TYPE;
-- Definition:
--   EXP(X) ~ natural e raised to the X power
-- Usage:
--   Z := EXP(X);   -- e raised to the power X
-- Domain:
--   Mathematically unbounded
-- Range:
--   EXP(X) >= 0.0
-- Accuracy:
--   (a) Maximum relative error = 4.0 * FLOAT_TYPE'BASE'EPSILON
--   (b) EXP(0.0) = 1.0
-- Notes:
--   The usable domain of EXP is approximately given by
--   X <= ln(FLOAT_TYPE'SAFE_LARGE)

  function EXP ( X : FLOAT_TYPE ) return FLOAT_TYPE is
    RESULT : FLOAT_TYPE ;
  begin
    if X = 0.0 then
      return 1.0 ;
    elsif X > LOG_LAST then
      raise ARGUMENT_ERROR ;
    end if ;
    RESULT := FLOAT_TYPE ( MATH_LIB.EXP ( FLOAT ( X ) ) ) ;
    if RESULT <= 0.0 then
      return 0.0 ;
    end if ;
    return RESULT ;
  exception
    when others =>
      raise ARGUMENT_ERROR ;
  end EXP ;



-- "**"
-- Declaration:
--   function "**" (X, Y : FLOAT_TYPE) return FLOAT_TYPE;
-- Definition:
--   X ** Y ~ X raised to the power Y
-- Usage:
--   Z := X ** Y;   -- X raised to the power Y
-- Domain:
--   (a) X >= 0.0
--   (b) Y > 0.0 when X = 0.0
-- Range:
--   X ** Y >= 0.0
-- Accuracy:
--   (a) Maximum relative error (when X > 0.0) =
--   (4.0+|Y*ln(X)|/32.0) * FLOAT_TYPE'BASE'EPSILON
--   (b) X ** 0.0 = 1.0 when X > 0.0
--   (c) 0.0 ** Y = 0.0 when Y > 0.0
--   (d) X ** 1.0 = X
--   (e) 1.0 ** Y =1.0
-- Notes:
--   The usable domain of "**", when X > 0.0, is approximately the set of
--   values for X and Y satisfying
--   Y*ln(X) <= ln(FLOAT_TYPE'SAFE_LARGE)
--   This imposes a positive upper bound on Y (as a function of X) when
--   X > 1.0, and a negative lower bound on Y (as a function of X) when
--   0.0 < X < 1.0.

  function "**" ( X , Y : FLOAT_TYPE ) return FLOAT_TYPE is
    RESULT : FLOAT_TYPE ;
  begin
    if X = 0.0 and then Y = 0.0 then
      raise ARGUMENT_ERROR ;
    elsif X < 0.0 then
      raise ARGUMENT_ERROR ;
    elsif Y = 0.0 then
      return 1.0 ;
    elsif X = 0.0 then
      return 0.0 ;
    elsif X = 1.0 then
      return 1.0 ;
    elsif Y = 1.0 then
      return X ;
    elsif Y = 2.0 then
      return X * X ;
    end if ;
    RESULT := FLOAT_TYPE ( MATH_LIB.EXP ( FLOAT ( Y ) * MATH_LIB.LN (
        FLOAT ( X ) ) ) );
    if RESULT <= 0.0 then
      return 0.0 ;
    end if ;
    return RESULT ;
  exception
    when others =>
      raise ARGUMENT_ERROR ;
  end "**" ;


-- SIN (natural cycle)
-- Declaration:
--   function SIN (X : FLOAT_TYPE) return FLOAT_TYPE;
-- Definition:
--   SIN(X) ~ sin(X)
-- Usage:
--   Z := SIN(X);   -- X in radians
-- Domain:
--   Mathematically unbounded
-- Range:
--   |SIN(X)| <= 1.0
-- Accuracy:
--   (a) Maximum relative error = 2.0 * FLOAT_TYPE'BASE'EPSILON
--   when |X| is less than or equal to some documented implementation-dependent
--   threshold, which must not be less than
--     FLOAT_TYPE'MACHINE_RADIX ** FLOAT_TYPE'MACHINE_MANTISSA/2
--   For larger values of |X|, degraded accuracy is allowed. An implementation
--   must document its behavior for large |X|.
--   (b) SIN(0.0) = 0.0

  function SIN ( X : FLOAT_TYPE ) return FLOAT_TYPE is
  begin
    if X = 0.0 then
      return 0.0 ;
    elsif abs X < SQUARE_ROOT_EPSILON then
      return X ;
    end if ;
    return FLOAT_TYPE ( MATH_LIB.SIN ( FLOAT ( X ) ) ) ;
  end SIN ;

-- SIN (arbitrary cycle)
-- Declaration:
--   function SIN (X, CYCLE : FLOAT_TYPE) return FLOAT_TYPE;
-- Definition:
--   SIN(X,CYCLE) ~ sin(2Pi * X/CYCLE)
-- Usage:
--   Z := SIN(X, 360.0);   -- X in degrees
--   Z := SIN(X, 1.0);     -- X in bams (binary angular measure)
--   Z := SIN(X, CYCLE);   -- X in units such that one complete
--                         -- cycle of rotation corresponds to
--                         -- X = CYCLE
-- Domain:
--   (a) X mathematically unbounded
--   (b) CYCLE > 0.0
-- Range:
--   |SIN(X,CYCLE)| <= 1.0
-- Accuracy:
--   (a) Maximum relative error = 2.0 * FLOAT_TYPE'BASE'EPSILON
--   (b) For integer k, SIN(X,CYCLE)= 0.0 when X=k*CYCLE/2.0
--                                    1.0 when X=(4k+1)*CYCLE/4.0
--                                    -1.0 when X=(4k+3)*CYCLE/4.0

  function SIN ( X , CYCLE : FLOAT_TYPE ) return FLOAT_TYPE is
  begin
    if CYCLE <= 0.0 then
      raise ARGUMENT_ERROR ;
    elsif X = 0.0 then
      return 0.0 ;
    elsif X = CYCLE then
      return 0.0 ;
    end if ;
    if abs X > abs CYCLE then
      return FLOAT_TYPE ( MATH_LIB.SIN (
                    FLOAT(EXACT_REMAINDER(X,CYCLE)) * TWO_PI /
                                         FLOAT ( CYCLE ) ) ) ;
    else
      return FLOAT_TYPE ( MATH_LIB.SIN ( FLOAT( X ) * TWO_PI /
                                         FLOAT ( CYCLE ) ) ) ;
    end if ;
  end SIN ;



-- COS (natural cycle)
-- Declaration:
--   function COS (X : FLOAT_TYPE) return FLOAT_TYPE;
-- Definition:
--   COS(X) ~ cos(X)
-- Usage:
--   Z := COS(X);   -- X in radians
-- Domain:
--   Mathematically unbounded
-- Range:
--   |COS(X)| <= 1.0
-- Accuracy:
--   (a) Maximum relative error = 2.0 * FLOAT_TYPE'BASE'EPSILON
--   when |X| is less than or equal to some documented implementation-dependent
--   threshold, which must not be less than
--     FLOAT_TYPE'MACHINE_RADIX ** FLOAT_TYPE'MACHINE_MANTISSA/2
--   For larger values of |X|, degraded accuracy is allowed. An implementation
--   must document its behavior for large |X|.
--   (b) COS(0.0) = 1.0

  function COS ( X : FLOAT_TYPE ) return FLOAT_TYPE is
  begin
    if X = 0.0 then
      return 1.0 ;
    elsif abs X < SQUARE_ROOT_EPSILON then
      return 1.0 ;
    end if ;
    return FLOAT_TYPE ( MATH_LIB.COS ( FLOAT( X ) ) ) ;
  end COS ;



-- COS (arbitrary cycle)
-- Declaration:
--   function COS (X, CYCLE : FLOAT_TYPE) return FLOAT_TYPE;
-- Definition:
--   COS(X,CYCLE) ~ cos(2Pi*X/CYCLE)
-- Usage:
--   Z := COS(X, 360.0);   -- X in degrees
--   Z := COS(X, 1.0);     -- X in bams
--   Z := COS(X, CYCLE);   -- X in units such that one complete
--                         -- cycle of rotation corresponds to
--                         -- X = CYCLE
-- Domain:
--   (a) X mathematically unbounded
--   (b) CYCLE > 0.0
-- Range:
--   |COS(X,CYCLE)| <= 1.0
-- Accuracy:
--   (a) Maximum relative error = 2.0 * FLOAT_TYPE'BASE'EPSILON
--   (b) For integer k, COS(X,CYCLE) = 1.0 when X=k*CYCLE
--                                     0.0 when X=(2k+1)*CYCLE/4.0
--                                    -1.0 when X=(2k+1)*CYCLE/2.0

  function COS ( X , CYCLE : FLOAT_TYPE ) return FLOAT_TYPE is
  begin
    if CYCLE <= 0.0 then
      raise ARGUMENT_ERROR ;
    elsif X = 0.0 then
      return 1.0 ;
    elsif X = CYCLE then
      return 1.0 ;
    end if ;
    if abs X > abs CYCLE then
      return FLOAT_TYPE ( MATH_LIB.COS (
                    FLOAT(EXACT_REMAINDER(X,CYCLE)) * TWO_PI /
                                         FLOAT ( CYCLE ) ) ) ;
    else
      return FLOAT_TYPE ( MATH_LIB.COS ( FLOAT( X ) * TWO_PI /
                                         FLOAT ( CYCLE ) ) ) ;
    end if;
  end COS ;



-- TAN (natural cycle)
-- Declaration:
--   function TAN (X : FLOAT_TYPE) return FLOAT_TYPE;
-- Definition:
--   TAN(X) ~ tan(X)
-- Usage:
--   Z := TAN(X);   -- X in radians
-- Domain:
--   Mathematically unbounded
-- Range:
--   Mathematically unbounded
-- Accuracy:
--   (a) Maximum relative error = 4.0 * FLOAT_TYPE'BASE'EPSILON
--   when |X| is less than or equal to some documented implementation-dependent
--   threshold, which must not be less than
--     FLOAT_TYPE'MACHINE_RADIX ** FLOAT_TYPE'MACHINE_MANTISSA/2
--   For larger values of |X|, degraded accuracy is allowed. An implementation
--   must document its behavior for large |X|.
--   (b) TAN(0.0) = 0.0

  function TAN ( X : FLOAT_TYPE ) return FLOAT_TYPE is
  begin
    if X = 0.0 then
      return 0.0 ;
    elsif abs X < SQUARE_ROOT_EPSILON then
      return X ;
    end if ;
    return SIN ( X ) / COS ( X ) ;
  end TAN ;

-- TAN (arbitrary cycle)
-- Declaration:
--   function TAN (X, CYCLE : FLOAT_TYPE) return FLOAT_TYPE;
-- Definition:
--   TAN(X,CYCLE) ~ tan(2Pi*X/CYCLE)
-- Usage:
--   Z := TAN(X, 360.0);   -- X in degrees
--   Z := TAN(X, 1.0);     -- X in bams
--   Z := TAN(X, CYCLE);   -- X in units such that one complete
--                         -- cycle of rotation corresponds to
--                         -- X = CYCLE
-- Domain:
--   (a) X /= (2k+1)*CYCLE/4.0, for integer k
--   (b) CYCLE > 0.0
-- Range:
--   Mathematically unbounded
-- Accuracy:
--   (a) Maximum relative error = 4.0 * FLOAT_TYPE'BASE'EPSILON
--   (b) TAN(X,CYCLE) = 0.0 when X=k*CYCLE/2.0, for integer k

  function TAN ( X , CYCLE : FLOAT_TYPE ) return FLOAT_TYPE is
  begin
    if CYCLE <= 0.0 then
      raise ARGUMENT_ERROR ;
    elsif X = 0.0 then
      return 0.0 ;
    end if ;
    return SIN ( X , CYCLE ) / COS ( X , CYCLE ) ;
  end TAN ;



-- COT (natural cycle)
-- Declaration:
--   function COT (X : FLOAT_TYPE) return FLOAT_TYPE;
-- Definition:
--   COT(X) ~ cot(X)
-- Usage:
--   Z := COT(X);   -- X in radians
-- Domain:
--   X /= 0.0
-- Range:
--   Mathematically unbounded
-- Accuracy:
--   (a) Maximum relative error = 4.0 * FLOAT_TYPE'BASE'EPSILON
--   when |X| is less than or equal to some documented implementation-dependent
--   threshold, which must not be less than
--     FLOAT_TYPE'MACHINE_RADIX ** FLOAT_TYPE'MACHINE_MANTISSA/2
--   For larger values of |X|, degraded accuracy is allowed. An implementation
--   must document its behavior for large |X|.

  function COT ( X : FLOAT_TYPE ) return FLOAT_TYPE is
  begin
    if X = 0.0 then
      raise ARGUMENT_ERROR ;
    end if ;
    return COS ( X ) / SIN ( X ) ;
  end COT ;


-- COT (arbitrary cycle)
-- Declaration:
--   function COT (X, CYCLE : FLOAT_TYPE) return FLOAT_TYPE;
-- Definition:
--   COT(X,CYCLE) ~ cot(2Pi*X/CYCLE)
-- Usage:
--   Z := COT(X, 360.0);   -- X in degrees
--   Z := COT(X, 1.0);     -- X in bams
--   Z := COT(X, CYCLE);   -- X in units such that one complete
--                         -- cycle of rotation corresponds to
--                         -- X = CYCLE
-- Domain:
--   (a) X /= k*CYCLE/2.0, for integer k
--   (b) CYCLE > 0.0
-- Range:
--   Mathematically unbounded
-- Accuracy:
--   (a) Maximum relative error = 4.0 * FLOAT_TYPE'BASE'EPSILON
--   (b) COT(X,CYCLE) = 0.0 when X=(2k+1)*CYCLE/4.0, for integer k

  function COT ( X , CYCLE : FLOAT_TYPE ) return FLOAT_TYPE is
  begin
    if CYCLE <= 0.0 then
      raise ARGUMENT_ERROR ;
    elsif X = 0.0 then
      raise ARGUMENT_ERROR ;
    end if ;
    return  COS ( X , CYCLE ) / SIN ( X , CYCLE ) ;
  end COT ;



-- ARCSIN (natural cycle)
-- Declaration:
--   function ARCSIN (X : FLOAT_TYPE) return FLOAT_TYPE;
-- Definition:
--   ARCSIN(X) ~ arcsine(X)
-- Usage:
--   Z := ARCSIN(X);   -- Z in radians
-- Domain:
--   |X| <= 1.0
-- Range:
--   |ARCSIN(X)| <= Pi/2
-- Accuracy:
--   (a) Maximum relative error = 4.0 * FLOAT_TYPE'BASE'EPSILON
--   (b) ARCSIN(0.0) = 0.0
--   (c) ARCSIN(1.0) = Pi/2
--   (d) ARCSIN(-1.0) = -Pi/2
-- Notes:
--   - Pi/2 and Pi/2 are not safe numbers of FLOAT_TYPE. Accordingly,
--   an implementation may exceed the range limits, but only slightly;
--   cf.Section 9 for a precise statement of the requirements. Similarly,
--   when accuracy requirement (c) or (d) applies, an implementation may
--   approximate the prescribed result, but only within narrow limits;
--   cf.Section 10 for a precise statement of the requirements.

  function ARCSIN ( X : FLOAT_TYPE ) return FLOAT_TYPE is
  begin
    if abs X > 1.0 then
      raise ARGUMENT_ERROR ;
    elsif X = 0.0 then
      return 0.0 ;
    elsif abs X < SQUARE_ROOT_EPSILON then
      return X ;
    elsif X = 1.0 then
      return PI / 2.0 ;
    elsif X = -1.0 then
      return -PI / 2.0 ;
    end if ;
    return ARCTAN ( X / SQRT ( 1.0 - X * X ) ) ;
  end ARCSIN ;


-- ARCSIN (arbitrary cycle)
-- Declaration:
--   function ARCSIN (X, CYCLE : FLOAT_TYPE) return FLOAT_TYPE;
-- Definition:
--   ARCSIN(X,CYCLE) ~ arcsin(X)*CYCLE/2Pi
-- Usage:
--   Z := ARCSIN(X, 360.0);   -- Z in degrees
--   Z := ARCSIN(X, 1.0);     -- Z in bams
--   Z := ARCSIN(X, CYCLE);   -- Z in units such that one complete
--                            -- cycle of rotation corresponds to
--                            -- Z = CYCLE
-- Domain:
--   (a) |X| <= 1.0
--   (b) CYCLE > 0.0
-- Range:
--   |ARCSIN(X,CYCLE) <= CYCLE/4.0
-- Accuracy:
--   (a) Maximum relative error = 4.0 * FLOAT_TYPE'BASE'EPSILON
--   (b) ARCSIN(0.0,CYCLE) = 0.0
--   (c) ARCSIN(1.0,CYCLE) = CYCLE/4.0
--   (d) ARCSIN(-1.0,CYCLE) = -CYCLE/4.0
-- Notes:
--   - CYCLE/4.0 and CYCLE/4.0
--   might not be safe numbers of FLOAT_TYPE. Accordingly,
--   an implementation may exceed the range limits, but only slightly;
--   cf.Section 9 for a precise statement of the requirements. Similarly,
--   when accuracy requirement (c) or (d) applies, an implementation may
--   approximate the prescribed result, but only within narrow limits;
--   cf.Section 10 for a precise statement of the requirements.

  function ARCSIN ( X , CYCLE : FLOAT_TYPE ) return FLOAT_TYPE is
  begin
    if CYCLE <= 0.0 then
      raise ARGUMENT_ERROR ;
    elsif abs X > 1.0 then
      raise ARGUMENT_ERROR ;
    elsif X = 0.0 then
      return 0.0 ;
    elsif X = 1.0 then
      return CYCLE / 4.0 ;
    elsif X = -1.0 then
      return -CYCLE / 4.0 ;
    end if ;
    return ARCTAN ( X / SQRT ( 1.0 - X * X ) , 1.0 , CYCLE ) ;
  end ARCSIN ;


-- ARCCOS (natural cycle)
-- Declaration:
--   function ARCCOS (X : FLOAT_TYPE) return FLOAT_TYPE;
-- Definition:
--   ARCCOS(X) ~ arccos(X)
-- Usage:
--   Z := ARCCOS(X);   -- Z in radians
-- Domain:
--   |X| <= 1.0
-- Range:
--   0.0 <= ARCCOS(X) <= Pi
-- Accuracy:
--   (a)  Maximum relative error = 4.0 * FLOAT_TYPE'BASE'EPSILON
--   (b) ARCCOS(1.0) = 0.0
--   (c) ARCCOS(0.0) = Pi/2
--   (d) ARCCOS(-1.0) = Pi
-- Notes:
--   Pi/2 and Pi are not safe numbers of FLOAT_TYPE. Accordingly,
--   an implementation may exceed the range limits, but only slightly;
--   cf.Section 9 for a precise statement of the requirements. Similarly,
--   when accuracy requirement (c) or (d) applies, an implementation may
--   approximate the prescribed result, but only within narrow limits;
--   cf.Section 10 for a precise statement of the requirements.

  function ARCCOS ( X : FLOAT_TYPE ) return FLOAT_TYPE is
    TEMP : FLOAT_TYPE ;
  begin
    if abs X > 1.0 then
      raise ARGUMENT_ERROR ;
    elsif X = 0.0 then
      return PI / 2.0 ;
    elsif X = 1.0 then
      return 0.0 ;
    elsif X = -1.0 then
      return PI ;
    end if ;
    TEMP := ARCTAN ( SQRT ( 1.0 - X * X ) / X ) ;
    if TEMP < 0.0 then
      TEMP := PI + TEMP ;
    end if ;
    return TEMP ;
  end ARCCOS ;


-- ARCCOS (arbitrary cycle)
-- Declaration:
--   function ARCCOS (X, CYCLE : FLOAT_TYPE) return FLOAT_TYPE;
-- Definition:
--   ARCCOS(X,CYCLE) ~ arccos(X)*CYCLE/2Pi
-- Usage:
--   Z := ARCCOS(X, 360.0);   -- Z in degrees
--   Z := ARCCOS(X, 1.0);     -- Z in bams
--   Z := ARCCOS(X, CYCLE);   -- Z in units such that one complete
--                            -- cycle of rotation corresponds to
--                            -- Z = CYCLE
-- Domain:
--   (a) |X| <= 1.0
--   (b) CYCLE > 0.0
-- Range:
--   0.0 <= ARCCOS(X,CYCLE) <= CYCLE/2.0
-- Accuracy:
--   (a) Maximum relative error = 4.0 * FLOAT_TYPE'BASE'EPSILON
--   (b) ARCCOS(1.0,CYCLE) = 0.0
--   (c) ARCCOS(0.0,CYCLE) = CYCLE/4.0
--   (d) ARCCOS(-1.0,CYCLE) = CYCLE/2.0
-- Notes:
--   CYCLE/4.0 and CYCLE/2.0
--   might not be safe numbers of FLOAT_TYPE. Accordingly,
--   an implementation may exceed the range limits, but only slightly;
--   cf.Section 9 for a precise statement of the requirements. Similarly,
--   when accuracy requirement (c) or (d) applies, an implementation may
--   approximate the prescribed result, but only within narrow limits;
--   cf.Section 10 for a precise statement of the requirements.

  function ARCCOS ( X , CYCLE : FLOAT_TYPE ) return FLOAT_TYPE is
    TEMP : FLOAT_TYPE ;
  begin
  if CYCLE <= 0.0 then
      raise ARGUMENT_ERROR ;
    elsif abs X > 1.0 then
      raise ARGUMENT_ERROR ;
    elsif X = 0.0 then
      return CYCLE / 4.0 ;
    elsif X = 1.0 then
      return 0.0 ;
    elsif X = -1.0 then
      return CYCLE / 2.0 ;
    end if ;
    TEMP := ARCTAN ( SQRT ( 1.0 - X * X ) / X , 1.0 , CYCLE ) ;
    if TEMP < 0.0 then
      TEMP := CYCLE / 2.0 + TEMP ;
    end if ;
    return TEMP ;
  end ARCCOS ;



-- ARCTAN (natural cycle)
-- Declaration:
--   function ARCTAN (Y : FLOAT_TYPE;
--                    X : FLOAT_TYPE := 1.0) return FLOAT_TYPE;
-- Definition:
--   (a) ARCTAN(Y) ~ arctan(Y)
--   (b) ARCTAN(Y,X) ~ arctan(Y/X) when X >= 0.0
--                     arctan(Y/X)+Pi when X < 0.0 and Y >= 0.0
--                     arctan(Y/X)-Pi when X < 0.0 and Y < 0.0
-- Usage:
--   Z := ARCTAN(Y);      -- Z, in radians, is the angle (in the
--                        -- quadrant containing the point (1.0,Y))
--                        -- whose tangent is Y
--   Z := ARCTAN(Y, X);   -- Z, in radians, is the angle (in the
--                        -- quadrant containing the point (X,Y))
--                        -- whose tangent is Y/X
-- Domain:
--   X /= 0.0 when Y = 0.0
-- Range:
--     (a) |ARCTAN(Y)| <= Pi/2
--     (b) 0.0 < ARCTAN(Y,X) <= Pi when Y >= 0.0
--     (c) -Pi <= ARCTAN(Y,X) <= 0.0 when Y < 0.0
-- Accuracy:
--   (a) Maximum relative error = 4.0 * FLOAT_TYPE'BASE'EPSILON
--   (b) ARCTAN(0.0) = 0.0
--   (c) ARCTAN((0.0,X) = 0.0 when X > 0.0
--                        Pi when X < 0.0
--   (d) ARCTAN(Y,0.0) = Pi/2 when Y > 0.0
--                       -Pi/2 when Y < 0.0
-- Notes:
--   -Pi,-Pi/2,Pi/2 and Pi are not safe numbers of FLOAT_TYPE. Accordingly,
--   an implementation may exceed the range limits, but only slightly;
--   cf.Section 9 for a precise statement of the requirements. Similarly,
--   when accuracy requirement (c) or (d) applies, an implementation may
--   approximate the prescribed result, but only within narrow limits;
--   cf.Section 10 for a precise statement of the requirements.



 function ARCTAN ( Y : FLOAT_TYPE ;
                    X : FLOAT_TYPE := 1.0 ) return FLOAT_TYPE is
    T : FLOAT_TYPE ;
  begin
    if X = 0.0 and then Y = 0.0 then
      raise ARGUMENT_ERROR ;
    elsif Y = 0.0 then
      if X > 0.0 then
        return 0.0 ;
      else -- X < 0.0
        return PI ;
      end if;
    elsif X = 0.0 then
      if Y > 0.0 then
        return HALF_PI ;
      else -- Y < 0.0
        return - HALF_PI ;
      end if;
    else
      return LOCAL_ATAN ( Y , X ) ;
    end if ;
  end ARCTAN ;


-- ARCTAN (arbitrary cycle)
-- Declaration:
--   function ARCTAN (Y     : FLOAT_TYPE;
--                    X     : FLOAT_TYPE := 1.0;
--                    CYCLE : FLOAT_TYPE)        return FLOAT_TYPE;
-- Definition:
--   (a) ARCTAN(Y,CYCLE) ~ arctan(Y)*CYCLE/2Pi
--   (b) ARCTAN(Y,X,CYCLE) ~ arctan(Y/X)*CYCLE/2Pi when X >= 0.0
--                           (arctan(Y/X)+Pi)*CYCLE/2Pi when X < 0.0 and Y >= 0.0
--                           (arctan(Y/X)-Pi)*CYCLE/2Pi when X < 0.0 and Y < 0.0
-- Usage:
--   Z := ARCTAN(Y, CYCLE => 360.0);   -- Z, in degrees, is the
--                                     -- angle (in the quadrant
--                                     -- containing the point
--                                     -- (1.0,Y)) whose tangent is Y
--   Z := ARCTAN(Y, CYCLE => 1.0);     -- Z, in bams, is the
--                                     -- angle (in the quadrant
--                                     -- containing the point
--                                     -- (1.0,Y)) whose tangent is Y
--   Z := ARCTAN(Y, CYCLE => CYCLE);   -- Z, in units such that one
--                                     -- complete cycle of rotation
--                                     -- corresponds to Z = CYCLE,
--                                     -- is the angle (in the
--                                     -- quadrant containing the
--                                     -- point (1.0,Y)) whose
--                                     -- tangent is Y
--   Z := ARCTAN(Y, X, 360.0);         -- Z, in degrees, is the
--                                     -- angle (in the quadrant
--                                     -- containing the point (X,Y))
--                                     -- whose tangent is Y/X
--   Z := ARCTAN(Y, X, 1.0);           -- Z, in bams, is the
--                                     -- angle (in the quadrant
--                                     -- containing the point (X,Y))
--                                     -- whose tangent is Y/X
--   Z := ARCTAN(Y, X, CYCLE);         -- Z, in units such that one
--                                     -- complete cycle of rotation
--                                     -- corresponds to Z = CYCLE,
--                                     -- is the angle (in the
--                                     -- quadrant containing the
--                                     -- point (X,Y)) whose
--                                     -- tangent is Y/X
-- Domain:
--   (a) X /= 0.0 when Y = 0.0
--   (b) CYCLE > 0.0
-- Range:
--   (a) |ARCTAN(Y,CYCLE)| <= CYCLE/4.0
--   (b) 0.0 <= ARCTAN(Y,X,CYCLE) <= CYCLE/2.0 when Y >= 0.0
--   (c) -CYCLE/2.0 <= ARCTAN(Y,X,CYCLE) <= 0.0 when Y < 0.0
-- Accuracy:
--   (a) Maximum relative error = 4.0 * FLOAT_TYPE'BASE'EPSILON
--   (b) ARCTAN(0.0,CYCLE) = 0.0
--   (c) ARCTAN(0.0,X,CYCLE) = 0.0 when X > 0.0
--                             CYCLE/2.0 when X < 0.0
--   (d) ARCTAN(Y,0.0,CYCLE) = CYCLE/4.0 when Y > 0.0
--                             -CYCLE/4.0 when Y < 0.0
-- Notes:
--   -CYCLE/2.0, -CYCLE/4.0, CYCLE/4.0 and CYCLE/2.0
--   might not be safe numbers of FLOAT_TYPE. Accordingly,
--   an implementation may exceed the range limits, but only slightly;
--   cf.Section 9 for a precise statement of the requirements. Similarly,
--   when accuracy requirement (c) or (d) applies, an implementation may
--   approximate the prescribed result, but only within narrow limits;
--   cf.Section 10 for a precise statement of the requirements.

  function ARCTAN ( Y : FLOAT_TYPE ;
                    X : FLOAT_TYPE := 1.0 ;
                    CYCLE : FLOAT_TYPE ) return FLOAT_TYPE is
  begin
    if CYCLE <= 0.0 then
      raise ARGUMENT_ERROR ;
    elsif X = 0.0 and then Y = 0.0 then
      raise ARGUMENT_ERROR ;
    elsif Y = 0.0 then
      if X > 0.0 then
        return 0.0 ;
      else -- X < 0.0
        return CYCLE / 2.0 ;
      end if;
    elsif X = 0.0 then
      if Y > 0.0 then
        return CYCLE / 4.0 ;
      else -- Y < 0.0
        return -CYCLE / 4.0 ;
      end if;
    else
      return LOCAL_ATAN ( Y , X ) * CYCLE / TWO_PI ;
    end if ;
  end ARCTAN ;



-- ARCCOT (natural cycle)
-- Declaration:
--   function ARCCOT (X : FLOAT_TYPE;
--                    Y : FLOAT_TYPE := 1.0) return FLOAT_TYPE;
-- Definition:
--   (a) ARCCOT(X) ~ arccot(X)
--   (b) ARCCOT(X,Y) ~ arccot(X/Y) when Y >= 0.0
--                     arccot(X/Y)-Pi when Y < 0.0
-- Usage:
--   Z := ARCCOT(X);      -- Z, in radians, is the angle (in the
--                        -- quadrant containing the point (X,1.0)
--                        -- whose cotangent is X
--   Z := ARCCOT(X, Y);   -- Z, in radians, is the angle (in the
--                        -- quadrant containing the point (X,Y))
--                        -- whose cotangent is X/Y
-- Domain:
--   Y /= 0.0 when X = 0.0
-- Range:
--   (a) 0.0 <= ARCCOT(X) <= Pi
--   (b) 0.0 <= ARCCOT(X,Y) <= Pi when Y >= 0.0
--   (c) -Pi <= ARCCOT(X,Y) <= 0.0 when Y < 0.0
-- Accuracy:
--   (a) Maximum relative error = 4.0 * FLOAT_TYPE'BASE'EPSILON
--   (b) ARCCOT(0.0) = Pi/2
--   (c) ARCCOT(0.0,Y) = Pi/2 when Y > 0.0
--                       -Pi/2 when Y < 0.0
--   (d) ARCCOT(X,0.0) = 0.0 when X > 0.0
--                       Pi when X < 0.0
-- Notes:
--   -Pi,-Pi/2,Pi/2 and Pi are not safe numbers of FLOAT_TYPE. Accordingly,
--   an implementation may exceed the range limits, but only slightly;
--   cf.Section 9 for a precise statement of the requirements. Similarly,
--   when accuracy requirement (b), (c) or (d) applies, an implementation may
--   approximate the prescribed result, but only within narrow limits;
--   cf.Section 10 for a precise statement of the requirements.

  function ARCCOT ( X : FLOAT_TYPE ;
                    Y : FLOAT_TYPE := 1.0 ) return FLOAT_TYPE is
  begin
    return ARCTAN ( Y , X ) ;  -- just reverse arguments
  end ARCCOT ;



-- ARCCOT (arbitrary cycle)
-- Declaration:
--   function ARCCOT (X     : FLOAT_TYPE;
--                    Y     : FLOAT_TYPE := 1.0;
--                    CYCLE : FLOAT_TYPE)        return FLOAT_TYPE;
-- Definition:
--   (a) ARCCOT(X,CYCLE) ~ arccot(X)*CYCLE/2Pi
--   (b) ARCCOT(X,Y) ~ arccot(X/Y)*CYCLE/2Pi when Y >= 0.0
--                     (arccot(X/Y)-Pi)*CYCLE/2Pi
-- Usage:
--   Z := ARCCOT(X, CYCLE => 360.0);   -- Z, in degrees, is the
--                                     -- angle (in the quadrant
--                                     -- containing the point
--                                     -- (X,1.0)) whose cotangent
--                                     -- is X
--   Z := ARCCOT(X, CYCLE => 1.0);     -- Z, in bams, is the
--                                     -- angle (in the quadrant
--                                     -- containing the point
--                                     -- (X,1.0)) whose cotangent
--                                     -- is X
--   Z := ARCCOT(X, CYCLE => CYCLE);   -- Z, in units such that one
--                                     -- complete cycle of rotation
--                                     -- corresponds to Z = CYCLE,
--                                     -- is the angle (in the
--                                     -- quadrant containing the
--                                     -- point (X,1.0)) whose
--                                     -- cotangent is X
--   Z := ARCCOT(X, Y, 360.0);         -- Z, in degrees, is the
--                                     -- angle (in the quadrant
--                                     -- containing the point (X,Y))
--                                     -- whose cotangent is X/Y
--   Z := ARCCOT(X, Y, 1.0);           -- Z, in bams, is the
--                                     -- angle (in the quadrant
--                                     -- containing the point (X,Y)
--                                     -- whose cotangent is X/Y
--   Z := ARCCOT(X, Y, CYCLE);         -- Z, in units such that one
--                                     -- complete cycle of rotation
--                                     -- corresponds to Z = CYCLE
--                                     -- is the angle (in the
--                                     -- quadrant containing the
--                                     -- point (X,Y)) whose
--                                     -- cotangent is X/Y
-- Domain:
--   (a) Y /= 0.0 when X = 0.0
--   (b) CYCLE > 0.0
-- Range:
--   (a) 0.0 <= ARCCOT((X,CYCLE) <= CYCLE/2.0
--   (b) 0.0 <= ARCCOT(X,Y,CYCLE) <= CYCLE/2.0 when Y >= 0.0
--   (c) -CYCLE/2.0 <= ARCCOT(X,Y,CYCLE) <= 0.0 when Y < 0.0
-- Accuracy:
--   (a) Maximum relative error = 4.0 * FLOAT_TYPE'BASE'EPSILON
--   (b) ARCCOT(0.0,CYCLE) = CYCLE/4.0
--   (c) ARCCOT(0.0,Y,CYCLE) = CYCLE/4.0 when Y > 0.0
--                             -CYCLE/4.0 when Y < 0.0
--   (d) ARCCOT(X,0.0,CYCLE) = 0.0 when X > 0.0
--                             CYCLE/2.0 when X < 0.0
-- Notes:
--   - CYCLE/2.0, - CYCLE/4.0, CYCLE/4.0 and CYCLE/2.0
--   might not be safe numbers of FLOAT_TYPE. Accordingly,
--   an implementation may exceed the range limits, but only slightly;
--   cf.Section 9 for a precise statement of the requirements. Similarly,
--   when accuracy requirement (c) or (d) applies, an implementation may
--   approximate the prescribed result, but only within narrow limits;
--   cf.Section 10 for a precise statement of the requirements.

  function ARCCOT ( X : FLOAT_TYPE ;
                    Y : FLOAT_TYPE := 1.0 ;
                    CYCLE : FLOAT_TYPE ) return FLOAT_TYPE is
  begin
    return ARCTAN ( Y , X , CYCLE ) ;  -- just reverse arguments
  end ARCCOT ;


-- SINH
-- Declaration:
--   function SINH (X : FLOAT_TYPE) return FLOAT_TYPE;
-- Definition:
--   SINH(X) ~ hyperbolic sine of X
-- Usage:
--   Z := SINH(X);
-- Domain:
--   Mathematically unbounded
-- Range:
--   Mathematically unbounded
-- Accuracy:
--   (a) Maximum relative error = 8.0 * FLOAT_TYPE'BASE'EPSILON
--   (b) SINH(0.0) = 0.0
-- Notes:
--   The usable domain of SINH is approximately given by
--   |X| <= ln(FLOAT_TYPE'SAFE_LARGE)+ln(2.0)

  function SINH ( X : FLOAT_TYPE ) return FLOAT_TYPE is
    EXP_X : FLOAT_TYPE ;
  begin
    if X = 0.0 then
      return 0.0 ;
    elsif abs X < SQUARE_ROOT_EPSILON then
      return X ;
    elsif abs X > LOG_LAST + LOG_TWO then
      raise ARGUMENT_ERROR ;
    elsif  X > LOG_INVERSE_EPSILON then
      return EXP ( X - LOG_TWO ) ;
    elsif X < -LOG_INVERSE_EPSILON then
      return - EXP ( (-X) - LOG_TWO ) ;
    end if;
    EXP_X := EXP ( X ) ;
    return 0.5 * ( EXP_X - 1.0 / EXP_X ) ;
  exception
    when others =>
      raise ARGUMENT_ERROR ;
  end SINH ;


-- COSH
-- Declaration:
--   function COSH (X : FLOAT_TYPE) return FLOAT_TYPE;
-- Definition:
--   COSH(X) ~ hyperbolic cosine of X
-- Usage:
--   Z := COSH(X);
-- Domain:
--   Mathematically unbounded
-- Range:
--   COSH(X) >= 1.0
-- Accuracy:
--   (a) Maximum relative error = 8.0 * FLOAT_TYPE'BASE'EPSILON
--   (b) COSH(0.0) = 1.0
-- Notes:
--   The usable domain of COSH is approximately given by
--   |X| <= ln(FLOAT_TYPE'SAFE_LARGE)+ln(2.0)

  function COSH ( X : FLOAT_TYPE ) return FLOAT_TYPE is
    EXP_X : FLOAT_TYPE ;
  begin
    if X = 0.0 then
      return 1.0 ;
    elsif abs X < SQUARE_ROOT_EPSILON then
      return 1.0 ;
    elsif abs X > LOG_LAST + LOG_TWO then
      raise ARGUMENT_ERROR ;
    elsif  abs X > LOG_INVERSE_EPSILON then
      return EXP ( (abs X) - LOG_TWO ) ;
    end if ;
    EXP_X := EXP ( X ) ;
    return 0.5 * ( EXP_X + 1.0 / EXP_X ) ;
  exception
    when others =>
      raise ARGUMENT_ERROR ;
  end COSH ;



-- TANH
-- Declaration:
--   function TANH (X : FLOAT_TYPE) return FLOAT_TYPE;
-- Definition:
--   TANH(X) ~ hyperbolic tangent of X
-- Usage:
--   Z := TANH(X);
-- Domain:
--   Mathematically unbounded
-- Range:
--   |TANH(X)| <= 1.0
-- Accuracy:
--   (a) Maximum relative error = 8.0 * FLOAT_TYPE'BASE'EPSILON
--   (b) TANH(0.0) = 0.0

  function TANH ( X : FLOAT_TYPE ) return FLOAT_TYPE is
    EXP_2X : FLOAT_TYPE ;
  begin
    if X < HALF_LOG_EPSILON then
      return - 1.0 ;
    elsif X > - HALF_LOG_EPSILON then
      return 1.0 ;
    elsif X = 0.0 then
      return 0.0 ;
    elsif abs X < SQUARE_ROOT_EPSILON then
      return X ;
    end if ;
    EXP_2X := EXP ( - 2.0 * X ) ;
    return ( 1.0 - EXP_2X ) / ( 1.0 + EXP_2X ) ;
  end TANH ;



-- COTH
-- Declaration:
--   function COTH (X : FLOAT_TYPE) return FLOAT_TYPE;
-- Definition:
--   COTH(X) ~ hyperbolic cotangent of X
-- Usage:
--   Z := COTH(X);
-- Domain:
--   X /= 0.0
-- Range:
--   |COTH(X)| >= 1.0
-- Accuracy:
--   Maximum relative error = 8.0 * FLOAT_TYPE'BASE'EPSILON

  function COTH ( X : FLOAT_TYPE ) return FLOAT_TYPE is
    EXP_2X : FLOAT_TYPE ;
  begin
    if X < HALF_LOG_EPSILON then
      return - 1.0 ;
    elsif X > - HALF_LOG_EPSILON then
      return 1.0 ;
    elsif X = 0.0 then
      raise ARGUMENT_ERROR ;
    elsif abs X < SQUARE_ROOT_EPSILON then
      return 1.0 / X ;
    end if ;
    EXP_2X := EXP ( - 2.0 * X ) ;
    return ( 1.0 + EXP_2X ) / ( 1.0 - EXP_2X ) ;
  end COTH ;



-- ARCSINH
-- Declaration:
--   function ARCSINH (X : FLOAT_TYPE) return FLOAT_TYPE;
-- Definition:
--   ARCSINH(X) ~ hyperbolic arcsine of X
-- Usage:
--   Z := ARCSINH(X);
-- Domain:
--   Mathematically unbounded
-- Range:
--   Mathematically unbounded
-- Accuracy:
--   (a) Maximum relative error = 8.0 * FLOAT_TYPE'BASE'EPSILON
--   (b) ARCSINH(0.0) = 0.0
-- Notes:
--   The reachable range of ARCSINH is approximately given by
--     |ARCSINH(X)| <= ln(FLOAT_TYPE'SAFE_LARGE)+ln(2.0)
--   Maximum relative error = 8.0 * FLOAT_TYPE'BASE'EPSILON

  function ARCSINH ( X : FLOAT_TYPE ) return FLOAT_TYPE is
  begin
    if abs X < SQUARE_ROOT_EPSILON then
      return X ;
    elsif X > 1.0 / SQUARE_ROOT_EPSILON then
      return LOG ( X ) + LOG_TWO ;
    elsif X < - 1.0 / SQUARE_ROOT_EPSILON then
      return - ( LOG ( - X ) + LOG_TWO ) ;
    elsif X < 0.0 then
      return -LOG ( abs X + SQRT( X * X + 1.0 )) ;
    else
      return LOG ( X + SQRT( X * X + 1.0 )) ;
    end if ;
  end ARCSINH ;


-- ARCCOSH
-- Declaration:
--   function ARCCOSH (X : FLOAT_TYPE) return FLOAT_TYPE;
-- Definition:
--   ARCCOSH(X) ~ hyperbolic arccosine of X
-- Usage:
--   Z := ARCCOSH(X);
-- Domain:
--   X >= 1.0
-- Range:
--   ARCCOSH(X) >= 0.0
-- Accuracy:
--   (a) Maximum relative error = 8.0 * FLOAT_TYPE'BASE'EPSILON
--   (b) ARCCOSH(1.0) = 0.0
-- Notes:
--   The upper bound of the reachable range of ARCCOSH is approximately given
--   by  ARCCOSH(X) <= ln(FLOAT_TYPE'SAFE_LARGE)+ln(2.0)
--   Maximum relative error = 8.0 * FLOAT_TYPE'BASE'EPSILON

  function ARCCOSH ( X : FLOAT_TYPE ) return FLOAT_TYPE is
  begin
  --    return LOG ( X - SQRT( X * X - 1.0 )) ;  double valued,
  --                                             only positive value returned
    if X < 1.0 then
      raise ARGUMENT_ERROR ;
    elsif X < 1.0 + SQUARE_ROOT_EPSILON then
      return X - 1.0 ;
    elsif abs X > 1.0 / SQUARE_ROOT_EPSILON then
      return LOG ( X ) + LOG_TWO ;
    else
      return LOG ( X + SQRT( X * X - 1.0 )) ;
    end if ;
  end ARCCOSH ;



-- ARCTANH
-- Declaration:
--   function ARCTANH (X : FLOAT_TYPE) return FLOAT_TYPE;
-- Definition:
--   ARCTANH(X) ~ hyperbolic arctangent of X
-- Usage:
--   Z := ARCTANH(X);
-- Domain:
--   |X| < 1.0
-- Range:
--   Mathematically unbounded
-- Accuracy:
--   (a) Maximum relative error = 8.0 * FLOAT_TYPE'BASE'EPSILON
--   (b) ARCTANH(0.0) = 0.0

  function ARCTANH ( X : FLOAT_TYPE ) return FLOAT_TYPE is
  begin
    if abs X >= 1.0 then
      raise ARGUMENT_ERROR ;
    elsif abs X < EPSILON then
      return X ;
    else
      return 0.5 * LOG (( 1.0 + X ) / ( 1.0 - X )) ;
    end if ;
  end ARCTANH ;


-- ARCCOTH
-- Declaration:
--   function ARCCOTH (X : FLOAT_TYPE) return FLOAT_TYPE;
-- Definition:
--   ARCCOTH(X) ~ hyperbolic arc cotangent X
-- Usage:
--   Z := ARCCOTH(X);
-- Domain:
--   |X| > 1.0
-- Range:
--   Mathematically unbounded
-- Accuracy:
--   Maximum relative error = 8.0 * FLOAT_TYPE'BASE'EPSILON

  function ARCCOTH ( X : FLOAT_TYPE ) return FLOAT_TYPE is
  begin
    if abs X <= 1.0 then
      raise ARGUMENT_ERROR ;
    elsif abs X > 1.0 / EPSILON then
      return 0.0 ;
    else
      return 0.5 * LOG (( 1.0 + X ) / ( X - 1.0 )) ;
    end if ;
  end ARCCOTH ;

begin -- pseudo constants computed during elaboration
  EPSILON := FLOAT_TYPE ( FLOAT_TYPE'MACHINE_RADIX ) ** ( - FLOAT_TYPE'
    MACHINE_MANTISSA ) ;
  PREVIOUS := EPSILON ;
  while EPSILON + 1.0 /= 1.0 loop
    PREVIOUS := EPSILON ;
    EPSILON := EPSILON / FLOAT_TYPE ( FLOAT_TYPE'MACHINE_RADIX ) ;
  end loop ;
  EPSILON := PREVIOUS ;
  SQUARE_ROOT_EPSILON := SQRT ( EPSILON ) ;
  LOG_INVERSE_EPSILON := LOG ( 1.0 / EPSILON ) ;
  LOG_LAST := LOG ( FLOAT_TYPE'LAST ) ;
  HALF_LOG_EPSILON := 0.5 * LOG ( EPSILON ) ;
end GENERIC_ELEMENTARY_FUNCTIONS ;