Programmer's ROM - The Computer Language Library/programmersrom.iso/ada/math/mlibts.ada

Text File | 1988-05-03 | 49.4 KB | 1,776 lines

--::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
--NUMIO.TXT
--::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

with TEXT_IO; use TEXT_IO;
package NUMERIC_IO is

  procedure GET(FILE : in FILE_TYPE; ITEM : out INTEGER);
  procedure GET(ITEM : out INTEGER);
  procedure GET(FILE : in FILE_TYPE; ITEM : out FLOAT);
  procedure GET(ITEM : out FLOAT);
  procedure PUT(FILE : in FILE_TYPE; ITEM : in INTEGER);
  procedure PUT(ITEM : in INTEGER; WIDTH : in FIELD);
  procedure PUT(ITEM : in INTEGER);
  procedure PUT(FILE : in FILE_TYPE; ITEM : in FLOAT);
  procedure PUT(ITEM : in FLOAT);

end NUMERIC_IO;


with TEXT_IO;
use TEXT_IO;
package body NUMERIC_IO is
-- This ought to be done by instantiating the FLOAT_IO and INTEGER_IO
--  But if you dont yet have the generic TEXT_IO implemented yet
--  then something like this does the job on the DEC-10 IAPC
--  But it is a kludge
--  No effort has been put into making it pretty or portable
use INTEGER_IO; use FLOAT_IO;

  procedure GET(FILE : in FILE_TYPE; ITEM : out INTEGER) is
  begin
INTEGER_IO.GET(FILE, ITEM);
  end GET;

  procedure GET(ITEM : out INTEGER) is
  begin
INTEGER_IO.GET(ITEM);
  end GET;

  procedure GET(FILE : in FILE_TYPE; ITEM : out FLOAT) is
  begin
FLOAT_IO.GET(FILE, ITEM);
  end GET;

  procedure GET(ITEM : out FLOAT) is
  begin
FLOAT_IO.GET(ITEM);
  end GET;

  procedure PUT(FILE : in FILE_TYPE; ITEM : in INTEGER) is
  begin
INTEGER_IO.PUT(FILE, ITEM);
  end PUT;

  procedure PUT(ITEM : in INTEGER; WIDTH : in FIELD) is
    J, K, M : INTEGER := 0;
  begin
    if WIDTH = 1  then
      case ITEM is
        when 0  => PUT('0');
        when 1  => PUT('1');
        when 2  => PUT('2');
        when 3  => PUT('3');
        when 4  => PUT('4');
        when 5  => PUT('5');
        when 6  => PUT('6');
        when 7  => PUT('7');
        when 8  => PUT('8');
        when 9  => PUT('9');
        when others  => PUT('*');
      end case;
    else
      if ITEM < 0  then
        PUT('-');
        J := -ITEM;
      else
        PUT(' ');
        J := ITEM;
      end if;
      for I in 1..WIDTH-1  loop
        M := 10**(WIDTH - 1 - I);
        K := J / M;
        J := J - K*M;
        NUMERIC_IO.PUT(K, 1);
      end loop;
    end if;
  end PUT;

  procedure PUT(ITEM : in INTEGER) is
  begin
INTEGER_IO.PUT(ITEM);
  end PUT;

  procedure PUT(FILE : in FILE_TYPE; ITEM : in FLOAT) is
  begin
FLOAT_IO.PUT(FILE, ITEM);
  end PUT;

  procedure PUT(ITEM : in FLOAT) is
  begin
FLOAT_IO.PUT(ITEM);
  end PUT;

end NUMERIC_IO;



--::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
--FLOATCH.TXT
--::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

package FLOATING_CHARACTERISTICS is
--  This package is a floating mantissa definition of a binary FLOAT
--  It was first used on the DEC-10 and the VAX but should work for any
--  since the parameters are obtained by initializing on the actual hardware
--  Otherwise the parameters could be set in the spec if known
--  This is a preliminary package that defines the properties
--  of the particular floating point type for which we are going to
--  generate the math routines
--  The constants are those required by the routines described in
--  "Software Manual for the Elementary Functions" W. Cody & W. Waite
--  Prentice-Hall 1980
--  Actually most are needed only for the test programs
--  rather than the functions themselves, but might as well be here
--  Most of these could be in the form of attributes if
--  all the floating types to be considered were those built into the
--  compiler, but we also want to be able to support user defined types
--  such as software floating types of greater precision than
--  the hardware affords, or types defined on one machine to
--  simulate another
--  So we use the Cody-Waite names and derive them from an adaptation of the
--  MACHAR routine as given by Cody-Waite in Appendix B

    IBETA : INTEGER;
    --  The radix of the floating-point representation

    IT : INTEGER;
    --  The number of base IBETA digits in the DIS_FLOAT significand

    IRND : INTEGER;
    --  TRUE (1) if floating addition rounds, FALSE (0) if truncates

    NGRD : INTEGER;
    --  Number of guard digits for multiplication

    MACHEP : INTEGER;
    --  The largest negative integer such that
    --    1.0 + FLOAT(IBETA) ** MACHEP /= 1.0
    --  except that MACHEP is bounded below by -(IT + 3)

    NEGEP : INTEGER;
    --  The largest negative integer such that
    --    1.0 -0 FLOAT(IBETA) ** NEGEP /= 1.0
    --  except that NEGEP is bounded below by -(IT + 3)

    IEXP : INTEGER;
    --  The number of bits (decimal places if IBETA = 10)
    --  reserved for the representation of the exponent (including
    --  the bias or sign) of a floating-point number

    MINEXP : INTEGER;
    --  The largest in magnitude negative integer such that
    --  FLOAT(IBETA) ** MINEXP is a positive floating-point number


    MAXEXP : INTEGER;
    --  The largest positive exponent for a finite floating-point number

    EPS : FLOAT;
    --  The smallest positive floating-point number such that
    --                              1.0 + EPS /= 1.0
    --  In particular, if IBETA = 2 or IRND = 0,
    --  EPS = FLOAT(IBETA) ** MACHEP
    --  Otherwise, EPS = (FLOAT(IBETA) ** MACHEP) / 2


    EPSNEG : FLOAT;
    --  A small positive floating-point number such that 1.0-EPSNEG /= 1.0

    XMIN : FLOAT;
    --  The smallest non-vanishing floating-point power of the radix
    --  In particular, XMIN = FLOAT(IBETA) ** MINEXP

    XMAX : FLOAT;
    --  The largest finite floating-point number

--  Here the structure of the floating type is defined
--  I have assumed that the exponent is always some integer form
--  The mantissa can vary
--  Most often it will be a fixed type or the same floating type
--  depending on the most efficient machine implementation
--  Most efficient implementation may require details of the machine hardware
--  In this version the simplest representation is used
--  The mantissa is extracted into a FLOAT and uses the predefined operations
  subtype EXPONENT_TYPE is INTEGER;    --  should be derived  ##########
  subtype MANTISSA_TYPE is FLOAT;     --   range -1.0..1.0;
--  A consequence of the rigorous constraints on MANTISSA_TYPE is that
--  operations must be very carefully examined to make sure that no number
--  greater than one results
--  Actually this limitation is important in constructing algorithms
--  which will also run when MANTISSA_TYPE is a fixed point type

--  If we are not using the STANDARD type, we have to define all the
--  operations at this point
--  We also need PUT for the type if it is not otherwise available

--  Now we do something strange
--  Since we do not know in the following routines whether the mantissa
--  will be carried as a fixed or floating type, we have to make some
--  provision for dividing by two
--  We cannot use the literals, since FIXED/2.0 and FLOAT/2 will fail
--  We define a type-dependent factor that will work
  MANTISSA_DIVISOR_2 : constant FLOAT := 2.0;
  MANTISSA_DIVISOR_3 : constant FLOAT := 3.0;
--  This will work for the MANTISSA_TYPE defined above
--  The alternative of defining an operation "/" to take care of it
--  is too sweeping and would allow unAda-like errors

  MANTISSA_HALF : constant MANTISSA_TYPE := 0.5;


  procedure DEFLOAT(X : in FLOAT;
                    N : out EXPONENT_TYPE; F : out MANTISSA_TYPE);
  procedure REFLOAT(N : in EXPONENT_TYPE; F : in MANTISSA_TYPE;
                                                   X : out FLOAT);
--  Since the user may wish to define a floating type by some other name
--  CONVERT_TO_FLOAT is used rather than just FLOAT for explicit coersion
  function CONVERT_TO_FLOAT(K : INTEGER) return FLOAT;
--  function CONVERT_TO_FLOAT(N : EXPONENT_TYPE) return FLOAT;
  function CONVERT_TO_FLOAT(F : MANTISSA_TYPE) return FLOAT;

end FLOATING_CHARACTERISTICS;



with TEXT_IO; use TEXT_IO;
package body FLOATING_CHARACTERISTICS is
--  This package is a floating mantissa definition of a binary FLOAT

    A, B, Y, Z : FLOAT;
    I, K, MX, IZ : INTEGER;
    BETA, BETAM1, BETAIN : FLOAT;
    ONE : FLOAT := 1.0;
    ZERO : FLOAT := 0.0;

  procedure DEFLOAT(X : in FLOAT;
                    N : out EXPONENT_TYPE; F : out MANTISSA_TYPE) is
--  This is admittedly a slow method - but portable - for breaking down
--  a floating point number into its exponent and mantissa
--  Obviously with knowledge of the machine representation
--  it could be replaced with a couple of simple extractions
    EXPONENT_LENGTH : INTEGER := IEXP;
    M : EXPONENT_TYPE;
    W, Y, Z : FLOAT;
  begin
    N := 0;
    F := 0.0;
    Y := ABS(X);
    if Y = 0.0  then
      return;
    elsif Y < 0.5  then
      for J in reverse 0..(EXPONENT_LENGTH - 2)  loop
      --  Dont want to go all the way to 2.0**(EXPONENT_LENGTH - 1)
      --  Since that (or its reciprocal) will overflow if exponent biased
      --  Ought to use talbular values rather than compute each time
        M := EXPONENT_TYPE(2 ** J);
        Z := 1.0 / (2.0**M);
        W := Y / Z;
        if W < 1.0  then
          Y := W;
          N := N - M;
        end if;
      end loop;
    else
      for J in reverse 0..(EXPONENT_LENGTH - 2)  loop
        M := EXPONENT_TYPE(2 ** J);
        Z := 2.0**M;
        W := Y / Z;
        if W >= 0.5  then
          Y := W;
          N := N + M;
        end if;
      end loop;
    --  And just to clear up any loose ends from biased exponents
    end if;
    while Y < 0.5  loop
      Y := Y * 2.0;
      N := N - 1;
    end loop;
    while Y >= 1.0  loop
      Y := Y / 2.0;
      N := N + 1;
    end loop;
    F := MANTISSA_TYPE(Y);
    if X < 0.0  then
      F := -F;
    end if;
    return;
  exception
  when others =>
    N := 0;
    F := 0.0;
    return;
  end DEFLOAT;


  procedure REFLOAT(N : in EXPONENT_TYPE; F : in MANTISSA_TYPE;
                                                   X : out FLOAT) is
--  Again a brute force method - but portable
--  Watch out near MAXEXP
    M : INTEGER;
    Y : FLOAT;
  begin
    if F = 0.0  then
      X := ZERO;
      return;
    end if;
    M := INTEGER(N);
    Y := ABS(FLOAT(F));
    while Y < 0.5  loop
      M := M - 1;
      if M < MINEXP  then
        X := ZERO;
      end if;
      Y := Y + Y;
      exit when M <= MINEXP;
    end loop;
    if M = MAXEXP  then
      M := M - 1;
      X := Y * 2.0**M;
      X := X * 2.0;
    elsif M <= MINEXP + 2  then
      M := M + 3;
      X := Y * 2.0**M;
      X := ((X / 2.0) / 2.0) / 2.0;
    else
      X := Y * 2.0**M;
    end if;
    if F < 0.0  then
      X := -X;
    end if;
    return;
  end REFLOAT;

  function CONVERT_TO_FLOAT(K : INTEGER) return FLOAT is
  begin
    return FLOAT(K);
  end CONVERT_TO_FLOAT;

  function CONVERT_TO_FLOAT(N : EXPONENT_TYPE) return FLOAT is
  begin
    return FLOAT(N);
  end CONVERT_TO_FLOAT;

  function CONVERT_TO_FLOAT(F : MANTISSA_TYPE) return FLOAT is
  begin
    return FLOAT(F);
  end CONVERT_TO_FLOAT;


begin
--  Initialization for the VAX with values derived by MACHAR
--  In place of running MACHAR as the actual initialization

    IBETA :=    2;
    IT :=    24;
    IRND :=    1;
    NEGEP :=    -24;
    EPSNEG :=    5.9604644E-008;
    MACHEP :=    -24;
    EPS :=    5.9604644E-008;
    NGRD :=    0;
    XMIN := 5.9E-39;
    MINEXP :=    -126;
    IEXP :=    8;
    MAXEXP :=    127;
    XMAX :=    8.5E37 * 2.0;


----  This initialization is the MACHAR routine of Cody and Waite Appendix B.
--PUT("INITIALIZATING WITH MACHAR     -     ");
--    A := ONE;
--    while (((A + ONE) - A) - ONE) = ZERO  loop
--      A := A + A;
--    end loop;
--    B := ONE;
--    while ((A + B) - A) = ZERO  loop
--      B := B + B;
--    end loop;
--    IBETA := INTEGER((A + B) - A);
--    BETA := CONVERT_TO_FLOAT(IBETA);
--
--
--    IT := 0;
--    B := ONE;
--    while (((B + ONE) - B) - ONE) = ZERO  loop
--      IT := IT + 1;
--      B := B * BETA;
--    end loop;
--
--
--    IRND := 0;
--    BETAM1 := BETA - ONE;
--    if ((A + BETAM1) - A) /= ZERO  then
--      IRND := 1;
--    end if;
--
--
--    NEGEP := IT + 3;
--    BETAIN := ONE / BETA;
--    A := ONE;
--  --  for I in 1..NEGEP  loop
--  for I in 1..50  loop
--  exit when I > NEGEP;
--      A := A * BETAIN;
--    end loop;
--    B := A;
--    while ((ONE - A) - ONE) = ZERO  loop
--      A := A * BETA;
--      NEGEP := NEGEP - 1;
--    end loop;
--    NEGEP := -NEGEP;
--
--
--    EPSNEG := A;
--    if (IBETA /= 2) and (IRND /= 0)  then
--      A := (A * (ONE + A)) / (ONE + ONE);
--      if ((ONE - A) - ONE) /= ZERO  then
--        EPSNEG := A;
--      end if;
--    end if;
--
--
--    MACHEP := -IT - 3;
--    A := B;
--    while ((ONE + A) - ONE) = ZERO  loop
--      A := A * BETA;
--      MACHEP := MACHEP + 1;
--    end loop;
--
--
--    EPS := A;
--    if (IBETA /= 2) and (IRND /= 0)  then
--      A := (A * (ONE + A)) / (ONE + ONE);
--      if ((ONE + A) - ONE) /= ZERO  then
--        EPS := A;
--      end if;
--    end if;
--
--
--    NGRD := 0;
--    if ((IRND = 0) and ((ONE + EPS) * ONE - ONE) /= ZERO)  then
--      NGRD := 1;
--    end if;
--
--
--    I := 0;
--    K := 1;
--    Z := BETAIN;
--    loop
--      Y := Z;
--      Z := Y * Y;
--      A := Z * ONE;
--      exit when ((A + A) = ZERO) or (ABS(Z) >= Y);
--      I := I + 1;
--      K := K + K;
--    end loop;
--    if (IBETA /= 10)  then
--      IEXP := I + 1;
--      MX := K + K;
--    else
--      IEXP := 2;
--      IZ := IBETA;
--      while (K >= IZ)  loop
--        IZ := IZ * IBETA;
--        IEXP := IEXP + 1;
--      end loop;
--      MX := IZ + IZ - 1;
--    end if;
--
--    loop
--      XMIN := Y;
--      Y := Y * BETAIN;
--      A := Y * ONE;
--      exit when ((A + A) = ZERO) or (ABS(Y) >= XMIN);
--      K := K + 1;
--    end loop;
--
--
--    MINEXP := -K;
--
--
--    if ((MX <= (K + K - 3)) and (IBETA /= 10))  then
--      MX := MX + MX;
--      IEXP := IEXP + 1;
--    end if;
--
--
--    MAXEXP := MX + MINEXP;
--    I := MAXEXP + MINEXP;
--    if ((IBETA = 2) and (I = 0))  then
--      MAXEXP := MAXEXP - 1;
--    end if;
--    if (I > 20)  then
--      MAXEXP := MAXEXP - 1;
--    end if;
--    if (A /= Y)  then
--      MAXEXP := MAXEXP - 2;
--    end if;
--
--
--    XMAX := ONE - EPSNEG;
--    if ((XMAX * ONE) /= XMAX)  then
--      XMAX := ONE - BETA * EPSNEG;
--    end if;
--    XMAX := XMAX / (BETA * BETA * BETA * XMIN);
--    I := MAXEXP + MINEXP + 3;
--    if I > 0  then
--      for J in 1..50  loop
--  exit when J > I;
--        if IBETA = 2  then
--          XMAX := XMAX + XMAX;
--        else
--          XMAX := XMAX * BETA;
--        end if;
--      end loop;
--    end if;
--
--PUT("INITIALIZED"); NEW_LINE;

end FLOATING_CHARACTERISTICS;




--::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
--NUMPM.TXT
--::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

with FLOATING_CHARACTERISTICS; use FLOATING_CHARACTERISTICS;
package NUMERIC_PRIMITIVES is

  --  This may seem a little much but is put in this form to allow the
  --  same form to be used for a generic package
  --  If that is not needed, simple litterals could be substituted
  ZERO  : FLOAT := CONVERT_TO_FLOAT(INTEGER(0));
  ONE   : FLOAT := CONVERT_TO_FLOAT(INTEGER(1));
  TWO   : FLOAT := ONE + ONE;
  THREE : FLOAT := ONE + ONE + ONE;
  HALF  : FLOAT := ONE / TWO;

  --  The following "constants" are effectively deferred to
  --  the initialization part of the package body
  --  This is in order to make it possible to generalize the floating type
  --  If that capability is not desired, constants may be included here
  PI            : FLOAT;
  ONE_OVER_PI   : FLOAT;
  TWO_OVER_PI   : FLOAT;
  PI_OVER_TWO   : FLOAT;
  PI_OVER_THREE : FLOAT;
  PI_OVER_FOUR  : FLOAT;
  PI_OVER_SIX   : FLOAT;


  function SIGN(X, Y : FLOAT) return FLOAT;
    --  Returns the value of X with the sign of Y
  function MAX(X, Y :  FLOAT) return FLOAT;
    --  Returns the algebraicly larger of X and Y
  function MIN(X, Y :  FLOAT) return FLOAT;
    --  Returns the algebraicly smaller of X and Y
  function TRUNCATE(X : FLOAT) return FLOAT;
    --  Returns the floating value of the integer no larger than X
    --  AINT(X)
  function ROUND(X : FLOAT) return FLOAT;
    --  Returns the floating value nearest X
    --  AINTRND(X)
  function RAN return FLOAT;
    --  This uses a portable algorithm and is included at this point
    --  Algorithms that presume unique machine hardware information
    --  should be initiated in FLOATING_CHARACTERISTICS

end NUMERIC_PRIMITIVES;



with FLOATING_CHARACTERISTICS; use FLOATING_CHARACTERISTICS;
package body NUMERIC_PRIMITIVES is


  function SIGN(X, Y : FLOAT) return FLOAT is
    --  Returns the value of X with the sign of Y
  begin
    if Y >= 0.0  then
      return X;
    else
      return -X;
    end if;
  end SIGN;

  function MAX(X, Y : FLOAT) return FLOAT is
  begin
    if X >= Y  then
      return X;
    else
      return Y;
    end if;
  end MAX;

  function MIN(X, Y : FLOAT) return FLOAT is
  begin
    if X <= Y  then
      return X;
    else
      return Y;
    end if;
  end MIN;

  function TRUNCATE(X : FLOAT) return FLOAT is
  --  Optimum code depends on how the system rounds at exact halves
  begin
    if FLOAT(INTEGER(X)) = X  then
      return X;
    end if;
    if X > ZERO  then
      return FLOAT(INTEGER(X - HALF));
    elsif X = ZERO  then
      return ZERO;
    else
      return FLOAT(INTEGER(X + HALF));
    end if;
  end TRUNCATE;

  function ROUND(X : FLOAT) return FLOAT is
  begin
    return FLOAT(INTEGER(X));
  end ROUND;


  package KEY is
    X : INTEGER := 10_001;
    Y : INTEGER := 20_001;
    Z : INTEGER := 30_001;
  end KEY;

  function RAN return FLOAT is
  --  This rectangular random number routine is adapted from a report
  --  "A Pseudo-Random Number Generator" by B. A. Wichmann and I. D. Hill
  --  NPL Report DNACS XX (to be published)
  --  In this stripped version, it is suitable for machines supporting
  --  INTEGER at only 16 bits and is portable in Ada
    W : FLOAT;
  begin

    KEY.X := 171 * (KEY.X mod 177 - 177) -  2 * (KEY.X / 177);
    if KEY.X < 0  then
      KEY.X := KEY.X + 30269;
    end if;

    KEY.Y := 172 * (KEY.Y mod 176 - 176) - 35 * (KEY.Y / 176);
    if KEY.Y < 0  then
      KEY.Y := KEY.Y + 30307;
    end if;

    KEY.Z := 170 * (KEY.Z mod 178 - 178) - 63 * (KEY.Z / 178);
    if KEY.Z < 0  then
      KEY.Z := KEY.Z + 30323;
    end if;

    --  CONVERT_TO_FLOAT is used instead of FLOAT since the floating
    --  type may be software defined

    W :=     CONVERT_TO_FLOAT(KEY.X)/30269.0
           + CONVERT_TO_FLOAT(KEY.Y)/30307.0
           + CONVERT_TO_FLOAT(KEY.Z)/30323.0;

    return  W - CONVERT_TO_FLOAT(INTEGER(W - 0.5));

  end RAN;

begin
  PI            := CONVERT_TO_FLOAT(INTEGER(3)) +
                   CONVERT_TO_FLOAT(MANTISSA_TYPE(0.14159_26535_89793_23846));
  ONE_OVER_PI   := CONVERT_TO_FLOAT(MANTISSA_TYPE(0.31830_98861_83790_67154));
  TWO_OVER_PI   := CONVERT_TO_FLOAT(MANTISSA_TYPE(0.63661_97723_67581_34308));
  PI_OVER_TWO   := CONVERT_TO_FLOAT(INTEGER(1)) +
                   CONVERT_TO_FLOAT(MANTISSA_TYPE(0.57079_63267_94896_61923));
  PI_OVER_THREE := CONVERT_TO_FLOAT(INTEGER(1)) +
                   CONVERT_TO_FLOAT(MANTISSA_TYPE(0.04719_75511_96597_74615));
  PI_OVER_FOUR  := CONVERT_TO_FLOAT(MANTISSA_TYPE(0.78539_81633_97448_30962));
  PI_OVER_SIX   := CONVERT_TO_FLOAT(MANTISSA_TYPE(0.52359_87755_98298_87308));

end NUMERIC_PRIMITIVES;




--::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
--COREF.TXT
--::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

with FLOATING_CHARACTERISTICS; use FLOATING_CHARACTERISTICS;
package CORE_FUNCTIONS is

  EXP_LARGE : FLOAT;
  EXP_SMALL : FLOAT;

  function SQRT(X : FLOAT) return FLOAT;
  function CBRT(X : FLOAT) return FLOAT;
  function LOG(X : FLOAT) return FLOAT;
  function LOG10(X : FLOAT) return FLOAT;
  function EXP(X : FLOAT) return FLOAT;
  function "**"(X, Y : FLOAT) return FLOAT;

end CORE_FUNCTIONS;




with TEXT_IO; use TEXT_IO;
with FLOATING_CHARACTERISTICS; use FLOATING_CHARACTERISTICS;
with NUMERIC_IO; use NUMERIC_IO;
with NUMERIC_PRIMITIVES; use NUMERIC_PRIMITIVES;
package body CORE_FUNCTIONS is

--  The following routines are coded directly from the algorithms and
--  coeficients given in "Software Manual for the Elementry Functions"
--  by William J. Cody, Jr. and William Waite, Prentice_Hall, 1980
--  CBRT by analogy
--  A more general formulation uses MANTISSA_TYPE, etc.
--  The coeficients are appropriate for 25 to 32 bits floating significance
--  They will work for less but slightly shorter versions are possible
--  The routines are coded to stand alone so they need not be compiled together

--  These routines have been coded to accept a general MANTISSA_TYPE
--  That is, they are designed to work with a manitssa either fixed of float
--  There are some explicit conversions which are required but these will
--  not cause any extra code to be generated

--      16 JULY 1982       W A WHITAKER  AFATL EGLIN AFB FL 32542
--                         T C EICHOLTZ  USAFA


  function SQRT(X : FLOAT) return FLOAT is
    M, N : EXPONENT_TYPE;
    F, Y : MANTISSA_TYPE;
    RESULT : FLOAT;

    subtype INDEX is INTEGER range 0..100;    --  #########################
    SQRT_L1 : INDEX := 3;
    --  Could get away with SQRT_L1 := 2 for 28 bits
    --  Using the better Cody-Waite coeficients overflows MANTISSA_TYPE
    SQRT_C1 : MANTISSA_TYPE := 8#0.3317777777#;
    SQRT_C2 : MANTISSA_TYPE := 8#0.4460000000#;
    SQRT_C3 : MANTISSA_TYPE := 8#0.55202_36314_77747_36311_0#;

  begin
    if X = ZERO  then
      RESULT := ZERO;
      return RESULT;
    elsif X = ONE  then            --  To get exact SQRT(1.0)
      RESULT := ONE;
      return RESULT;
    elsif X < ZERO  then
      NEW_LINE;
      PUT("CALLED SQRT FOR NEGATIVE ARGUMENT   ");
      PUT(X);
      PUT("   USED ABSOLUTE VALUE");
      NEW_LINE;
      RESULT := SQRT(ABS(X));
      return RESULT;
    else
      DEFLOAT(X, N, F);
      Y := SQRT_C1 + MANTISSA_TYPE(SQRT_C2 * F);
      for J in 1..SQRT_L1  loop
        Y := Y/MANTISSA_DIVISOR_2 + MANTISSA_TYPE((F/MANTISSA_DIVISOR_2)/Y);
      end loop;
      if (N mod 2) /= 0  then
        Y := MANTISSA_TYPE(SQRT_C3 * Y);
        N := N + 1;
      end if;
      M := N/2;
      REFLOAT(M,Y,RESULT);
      return RESULT;
    end if;
  exception
  when others =>
    NEW_LINE; PUT(" EXCEPTION IN SQRT, X = "); PUT(X);
    PUT("  RETURNED 1.0"); NEW_LINE;
    return ONE;
  end SQRT;


  function CBRT(X : FLOAT) return FLOAT is
    M, N : EXPONENT_TYPE;
    F, Y : MANTISSA_TYPE;
    RESULT : FLOAT;

    subtype INDEX is INTEGER range 0..100;    --  #########################
    CBRT_L1 : INDEX := 3;
    CBRT_C1 : MANTISSA_TYPE := 0.5874009;
    CBRT_C2 : MANTISSA_TYPE := 0.4125990;
    CBRT_C3 : MANTISSA_TYPE := 0.62996_05249;
    CBRT_C4 : MANTISSA_TYPE := 0.79370_05260;

  begin
    if X = ZERO then
      RESULT := ZERO;
      return RESULT;
    else
      DEFLOAT(X, N, F);
      F := ABS(F);
      Y := CBRT_C1 + MANTISSA_TYPE(CBRT_C2 * F);
      for J in 1 .. CBRT_L1 loop
        Y :=     Y
            - (  Y/MANTISSA_DIVISOR_3
               - MANTISSA_TYPE((F/MANTISSA_DIVISOR_3) / MANTISSA_TYPE(Y*Y)) );
      end loop;
      case (N mod 3) is
        when 0 =>
          null;
        when 1 =>
          Y := MANTISSA_TYPE(CBRT_C3 * Y);
          N := N + 2;
        when 2 =>
          Y := MANTISSA_TYPE(CBRT_C4 * Y);
          N := N + 1;
        when others =>
          null;
      end case;
      M := N/3;
      if X < ZERO  then
        Y := -Y;
      end if;
      REFLOAT(M, Y, RESULT);
      return RESULT;
    end if;
  exception
    when others =>
      RESULT := ONE;
      if X < ZERO then
      RESULT := - ONE;
      end if;
      NEW_LINE; PUT("EXCEPTION IN CBRT, X = "); PUT(X);
      PUT("  RETURNED  "); PUT(RESULT); NEW_LINE;
      return RESULT;
  end CBRT;

    function LOG(X : FLOAT) return FLOAT is
  --  Uses fixed formulation for generality

    RESULT : FLOAT;
    N : EXPONENT_TYPE;
    XN : FLOAT;
    Y : FLOAT;
    F : MANTISSA_TYPE;
    Z, ZDEN, ZNUM : MANTISSA_TYPE;

    C0 : constant MANTISSA_TYPE := 0.20710_67811_86547_52440;
                                               --  SQRT(0.5) - 0.5
    C1 : constant FLOAT := 8#0.543#;
    C2 : constant FLOAT :=-2.12194_44005_46905_82767_9E-4;

    function R(Z : MANTISSA_TYPE) return MANTISSA_TYPE is
    --  Use fixed formulation here because the float coeficents are > 1.0
    --  and would exceed the limits on a MANTISSA_TYPE
      A0 : constant MANTISSA_TYPE := 0.04862_85276_587;
      B0 : constant MANTISSA_TYPE := 0.69735_92187_803;
      B1 : constant MANTISSA_TYPE :=-0.125;
      C  : constant MANTISSA_TYPE := 0.01360_09546_862;
    begin
      return Z + MANTISSA_TYPE(Z *
          MANTISSA_TYPE(MANTISSA_TYPE(Z * Z) * (C +
          MANTISSA_TYPE(A0/(B0 + MANTISSA_TYPE(B1 * MANTISSA_TYPE(Z * Z)))))));
    end R;

  begin

    if X < ZERO      then
      NEW_LINE;
      PUT("CALLED LOG FOR NEGATIVE ");
      PUT(X);
      PUT("   USE ABS => ");
      RESULT := LOG(ABS(X));
      PUT(RESULT);
      NEW_LINE;
    elsif X = ZERO  then
      NEW_LINE;
      PUT("CALLED LOG FOR ZERO ARGUMENT, RETURNED ");
      RESULT := -XMAX;      --  SUPPOSED TO BE -LARGE
      PUT(RESULT);
      NEW_LINE;
    else
      DEFLOAT(X,N,F);
      ZNUM := F - MANTISSA_HALF;
      Y := CONVERT_TO_FLOAT(ZNUM);
      ZDEN := ZNUM / MANTISSA_DIVISOR_2 + MANTISSA_HALF;
      if ZNUM > C0  then
        Y := Y - MANTISSA_HALF;
        ZNUM := ZNUM - MANTISSA_HALF;
        ZDEN := ZDEN + MANTISSA_HALF/MANTISSA_DIVISOR_2;
      else
        N := N -1;
      end if;
      Z    := MANTISSA_TYPE(ZNUM / ZDEN);
      RESULT := CONVERT_TO_FLOAT(R(Z));
      if N /= 0  then
        XN := CONVERT_TO_FLOAT(N);
        RESULT := (XN * C2 + RESULT) + XN * C1;
      end if;
    end if;
    return RESULT;

  exception
  when others =>
    NEW_LINE; PUT(" EXCEPTION IN LOG, X = "); PUT(X);
    PUT("  RETURNED 0.0"); NEW_LINE;
    return ZERO;
  end LOG;


  function LOG10(X : FLOAT) return FLOAT is
    LOG_10_OF_2 : constant FLOAT :=
             CONVERT_TO_FLOAT(MANTISSA_TYPE(8#0.33626_75425_11562_41615#));
  begin
    return LOG(X) * LOG_10_OF_2;
  end LOG10;

  function EXP(X : FLOAT) return FLOAT is

    RESULT : FLOAT;
    N : EXPONENT_TYPE;
    XG, XN, X1, X2 : FLOAT;
    F, G : MANTISSA_TYPE;

    BIGX : FLOAT := EXP_LARGE;
    SMALLX : FLOAT := EXP_SMALL;

    ONE_OVER_LOG_2 : constant FLOAT :=  1.4426_95040_88896_34074;
    C1 : constant FLOAT :=  0.69335_9375;
    C2 : constant FLOAT := -2.1219_44400_54690_58277E-4;

    function R(G : MANTISSA_TYPE) return MANTISSA_TYPE is
      Z , GP, Q : MANTISSA_TYPE;

      P0 : constant MANTISSA_TYPE :=  0.24999_99999_9992;
      P1 : constant MANTISSA_TYPE :=  0.00595_04254_9776;
      Q0 : constant MANTISSA_TYPE :=  0.5;
      Q1 : constant MANTISSA_TYPE :=  0.05356_75176_4522;
      Q2 : constant MANTISSA_TYPE :=  0.00029_72936_3682;
    begin
      Z  := MANTISSA_TYPE(G * G);
      GP := MANTISSA_TYPE( (MANTISSA_TYPE(P1 * Z) + P0) * G );
      Q  := MANTISSA_TYPE( (MANTISSA_TYPE(Q2 * Z) + Q1) * Z ) + Q0;
      return MANTISSA_HALF + MANTISSA_TYPE( GP /(Q - GP) );
    end R;


  begin

    if X > BIGX  then
      NEW_LINE;
      PUT("  EXP CALLED WITH TOO BIG A POSITIVE ARGUMENT, ");
          PUT(X); PUT("   RETURNED XMAX");
      NEW_LINE;
      RESULT := XMAX;
    elsif X < SMALLX  then
      NEW_LINE;
      PUT("  EXP CALLED WITH TOO BIG A NEGATIVE ARGUMENT,  ");
          PUT(X); PUT("    RETURNED ZERO");
      NEW_LINE;
      RESULT := ZERO;
    elsif ABS(X) < EPS  then
      RESULT := ONE;
    else
      N  := EXPONENT_TYPE(X * ONE_OVER_LOG_2);
      XN := CONVERT_TO_FLOAT(N);
      X1 := ROUND(X);
      X2 := X - X1;
      XG := ( (X1 - XN * C1) + X2 ) - XN * C2;
      G  := MANTISSA_TYPE(XG);
      N  := N + 1;
      F := R(G);
      REFLOAT(N, F, RESULT);
    end if;
    return RESULT;

  exception
  when others =>
    NEW_LINE; PUT(" EXCEPTION IN EXP, X = "); PUT(X);
    PUT("  RETURNED 1.0"); NEW_LINE;
    return ONE;
  end EXP;

function "**" (X, Y : FLOAT) return FLOAT is
--  This is the last function to be coded since it appeared that it really
--  was un-Ada-like and ought not be in the regular package
--  Nevertheless it was included in this version
--  It is specific for FLOAT and does not have the MANTISSA_TYPE generality
  M, N : EXPONENT_TYPE;
  G : MANTISSA_TYPE;
  P, TEMP, IW1, I : INTEGER;
  RESULT, Z, V, R, U1, U2, W, W1, W2, W3, Y1, Y2 : FLOAT;

  K : constant FLOAT := 0.44269_50408_88963_40736;
  IBIGX : constant INTEGER := INTEGER(TRUNCATE(16.0 * LOG(XMAX) - 1.0));
  ISMALLX : constant INTEGER := INTEGER(TRUNCATE(16.0 * LOG(XMIN) + 1.0));

  P1 : constant FLOAT := 0.83333_32862_45E-1;
  P2 : constant FLOAT := 0.12506_48500_52E-1;

  Q1 : constant FLOAT := 0.69314_71805_56341;
  Q2 : constant FLOAT := 0.24022_65061_44710;
  Q3 : constant FLOAT := 0.55504_04881_30765E-1;
  Q4 : constant FLOAT := 0.96162_06595_83789E-2;
  Q5 : constant FLOAT := 0.13052_55159_42810E-2;

  A1 : array (1 .. 17) of FLOAT:=
     (  8#1.00000_0000#,
        8#0.75222_5750#,
        8#0.72540_3067#,
        8#0.70146_3367#,
        8#0.65642_3746#,
        8#0.63422_2140#,
        8#0.61263_4520#,
        8#0.57204_2434#,
        8#0.55202_3631#,
        8#0.53254_0767#,
        8#0.51377_3265#,
        8#0.47572_4623#,
        8#0.46033_7602#,
        8#0.44341_7233#,
        8#0.42712_7017#,
        8#0.41325_3033#,
        8#0.40000_0000#  );

  A2 : array (1 .. 8) of FLOAT :=
     (  8#0.00000_00005_22220_66302_61734_72062#,
        8#0.00000_00003_02522_47021_04062_61124#,
        8#0.00000_00005_21760_44016_17421_53016#,
        8#0.00000_00007_65401_41553_72504_02177#,
        8#0.00000_00002_44124_12254_31114_01243#,
        8#0.00000_00000_11064_10432_66404_42174#,
        8#0.00000_00004_72542_16063_30176_55544#,
        8#0.00000_00001_74611_03661_23056_22556#  );


  function REDUCE (V : FLOAT) return FLOAT is
  begin
    return FLOAT(INTEGER(16.0 * V)) * 0.0625;
  end REDUCE;

  begin
    if X <= ZERO then
      if X < ZERO then
        RESULT := (ABS(X))**Y;
        NEW_LINE;
        PUT("X**Y CALLED WITH X = "); PUT(X); NEW_LINE;
        PUT("USED ABS, RETURNED "); PUT(RESULT); NEW_LINE;
      else
        if Y <= ZERO then
          if Y = ZERO then
            RESULT := ZERO;
          else
            RESULT := XMAX;
          end if;
          NEW_LINE;
          PUT("X**Y CALLED WITH X = 0, Y = "); PUT(Y); NEW_LINE;
          PUT("RETURNED "); PUT(RESULT); NEW_LINE;
        else
          RESULT := ZERO;
        end if;
      end if;
    else
      DEFLOAT(X, M, G);
      P := 1;
      if G <= A1(9) then
        P := 9;
      end if;
      if G <= A1(P+4) then
        P := P + 4;
      end if;
      if G <= A1(P+2) then
        P := P + 2;
      end if;
      Z := ((G - A1(P+1)) - A2((P+1)/2))/(G + A1(P+1));
      Z := Z + Z;
      V := Z * Z;
      R := (P2 * V + P1) * V * Z;
      R := R + K * R;
      U2 := (R + Z * K) + Z;
      U1 := FLOAT(INTEGER(M) * 16 - P) * 0.0625;
      Y1 := REDUCE(Y);
      Y2 := Y - Y1;
      W := U2 * Y + U1 * Y2;
      W1 := REDUCE(W);
      W2 := W - W1;
      W := W1 + U1 * Y1;
      W1 := REDUCE(W);
      W2 := W2 + (W - W1);
      W3 := REDUCE(W2);
      IW1 := INTEGER(TRUNCATE(16.0 * (W1 + W3)));
      W2 := W2 - W3;
      if W > FLOAT(IBIGX) then
        RESULT := XMAX;
        PUT("X**Y CALLED  X ="); PUT(X); PUT("   Y ="); PUT(Y);
        PUT("   TOO LARGE  RETURNED "); PUT(RESULT); NEW_LINE;
      elsif W < FLOAT(ISMALLX) then
        RESULT := ZERO;
        PUT("X**Y CALLED  X ="); PUT(X); PUT("   Y ="); PUT(Y);
        PUT("   TOO SMALL  RETURNED "); PUT(RESULT); NEW_LINE;
      else
        if W2 > ZERO then
          W2 := W2 - 0.0625;
          IW1 := IW1 + 1;
        end if;
        if IW1 < INTEGER(ZERO) then
          I := 0;
        else
          I := 1;
        end if;
        M := EXPONENT_TYPE(I + IW1/16);
        P := 16 * INTEGER(M) - IW1;
        Z := ((((Q5 * W2 + Q4) * W2 + Q3) * W2 + Q2) * W2 + Q1) * W2;
        Z := A1(P+1) + (A1(P+1) * Z);

        REFLOAT(M, Z, RESULT);
      end if;
    end if;
    return RESULT;
  end "**";

begin
  EXP_LARGE := LOG(XMAX) * (ONE - EPS);
  EXP_SMALL := LOG(XMIN) * (ONE - EPS);
end CORE_FUNCTIONS;



--::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
--TRIGF.TXT
--::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

package TRIG_LIB is
  function SIN(X : FLOAT) return FLOAT;
  function COS(X : FLOAT) return FLOAT;
  function TAN(X : FLOAT) return FLOAT;
  function COT(X : FLOAT) return FLOAT;
  function ASIN(X : FLOAT) return FLOAT;
  function ACOS(X : FLOAT) return FLOAT;
  function ATAN(X : FLOAT) return FLOAT;
  function ATAN2(V, U : FLOAT) return FLOAT;
  function SINH(X : FLOAT) return FLOAT;
  function COSH(X : FLOAT) return FLOAT;
  function TANH(X : FLOAT) return FLOAT;
end TRIG_LIB;

with TEXT_IO; use TEXT_IO;
with FLOATING_CHARACTERISTICS; use FLOATING_CHARACTERISTICS;
with NUMERIC_IO; use NUMERIC_IO;
with NUMERIC_PRIMITIVES; use NUMERIC_PRIMITIVES;
with CORE_FUNCTIONS; use CORE_FUNCTIONS;
package body TRIG_LIB is

--  PRELIMINARY VERSION *********************************

--  The following routines are coded directly from the algorithms and
--  coeficients given in "Software Manual for the Elementry Functions"
--  by William J. Cody, Jr. and William Waite, Prentice_Hall, 1980
--  This particular version is stripped to work with FLOAT and INTEGER
--  and uses a mantissa represented as a FLOAT
--  A more general formulation uses MANTISSA_TYPE, etc.
--  The coeficients are appropriate for 25 to 32 bits floating significance
--  They will work for less but slightly shorter versions are possible
--  The routines are coded to stand alone so they need not be compiled together

--      16 JULY 1982       W A WHITAKER  AFATL EGLIN AFB FL 32542
--                         T C EICHOLTZ  USAFA


  function SIN(X : FLOAT) return FLOAT is
    SGN, Y : FLOAT;
    N : INTEGER;
    XN : FLOAT;
    F, G, X1, X2 : FLOAT;
    RESULT : FLOAT;

    YMAX : FLOAT := FLOAT(INTEGER(PI * TWO**(IT/2)));
    BETA : FLOAT := CONVERT_TO_FLOAT(IBETA);
    EPSILON : FLOAT := BETA ** (-IT/2);

    C1 : constant FLOAT :=  3.140625;
    C2 : constant FLOAT :=  9.6765_35897_93E-4;

    function R(G : FLOAT) return FLOAT is
      R1 : constant FLOAT := -0.16666_66660_883;
      R2 : constant FLOAT :=  0.83333_30720_556E-2;
      R3 : constant FLOAT := -0.19840_83282_313E-3;
      R4 : constant FLOAT :=  0.27523_97106_775E-5;
      R5 : constant FLOAT := -0.23868_34640_601E-7;
    begin
      return ((((R5*G + R4)*G + R3)*G + R2)*G + R1)*G;
    end R;

  begin
    if X < ZERO  then
      SGN := -ONE;
      Y := -X;
    else
      SGN := ONE;
      Y := X;
    end if;

    if Y > YMAX  then
      NEW_LINE;
      PUT(" SIN CALLED WITH ARGUMENT TOO LARGE FOR ACCURACY ");
      PUT(X); NEW_LINE;
    end if;

    N := INTEGER(Y * ONE_OVER_PI);
    XN := CONVERT_TO_FLOAT(N);
    if N mod 2 /= 0  then
      SGN := -SGN;
    end if;
    X1 := TRUNCATE(ABS(X));
    X2 := ABS(X) - X1;
    F := ((X1 - XN*C1) + X2) - XN*C2;
    if ABS(F) < EPSILON  then
      RESULT := F;
    else
      G := F * F;
      RESULT := F + F*R(G);
    end if;
    return (SGN * RESULT);
  end SIN;


  function COS(X : FLOAT) return FLOAT is
    SGN, Y : FLOAT;
    N : INTEGER;
    XN : FLOAT;
    F, G, X1, X2 : FLOAT;
    RESULT : FLOAT;

    YMAX : FLOAT := FLOAT(INTEGER(PI * TWO**(IT/2)));
    BETA : FLOAT := CONVERT_TO_FLOAT(IBETA);
    EPSILON : FLOAT := BETA ** (-IT/2);

    C1 : constant FLOAT :=  3.140625;
    C2 : constant FLOAT :=  9.6765_35897_93E-4;

    function R(G : FLOAT) return FLOAT is
      R1 : constant FLOAT := -0.16666_66660_883;
      R2 : constant FLOAT :=  0.83333_30720_556E-2;
      R3 : constant FLOAT := -0.19840_83282_313E-3;
      R4 : constant FLOAT :=  0.27523_97106_775E-5;
      R5 : constant FLOAT := -0.23868_34640_601E-7;
    begin
      return ((((R5*G + R4)*G + R3)*G + R2)*G + R1)*G;
    end R;

  begin
    SGN := 1.0;
    Y := ABS(X) + PI_OVER_TWO;

    if Y > YMAX  then
      NEW_LINE;
      PUT(" COS CALLED WITH ARGUMENT TOO LARGE FOR ACCURACY ");
      PUT(X); NEW_LINE;
    end if;

    N := INTEGER(Y * ONE_OVER_PI);
    XN := CONVERT_TO_FLOAT(N);
    if N mod 2 /= 0  then
      SGN := -SGN;
    end if;
    XN := XN - 0.5;          -- TO FORM COS INSTEAD OF SIN
    X1 := TRUNCATE(ABS(X));
    X2 := ABS(X) - X1;
    F := ((X1 - XN*C1) + X2) - XN*C2;
    if ABS(F) < EPSILON  then
      RESULT := F;
    else
      G := F * F;
      RESULT := F + F*R(G);
    end if;
    return (SGN * RESULT);
  end COS;


  function TAN(X : FLOAT) return FLOAT is
    SGN, Y : FLOAT;
    N : INTEGER;
    XN : FLOAT;
    F, G, X1, X2 : FLOAT;
    RESULT : FLOAT;

    YMAX : FLOAT := FLOAT(INTEGER(PI * TWO**(IT/2))) /2.0;
    BETA : FLOAT := CONVERT_TO_FLOAT(IBETA);
    EPSILON : FLOAT := BETA ** (-IT/2);

    C1 : constant FLOAT :=  8#1.444#;
    C2 : constant FLOAT :=  4.8382_67948_97E-4;

    function R(G : FLOAT) return FLOAT is
      P0 : constant FLOAT :=  1.0;
      P1 : constant FLOAT := -0.11136_14403_566;
      P2 : constant FLOAT :=  0.10751_54738_488E-2;
      Q0 : constant FLOAT :=  1.0;
      Q1 : constant FLOAT := -0.44469_47720_281;
      Q2 : constant FLOAT :=  0.15973_39213_300E-1;
    begin
      return ((P2*G + P1)*G*F + F) / (((Q2*G + Q1)*G +0.5) + 0.5);
    end R;

  begin
    Y := ABS(X);
    if Y > YMAX  then
      NEW_LINE;
      PUT(" TAN CALLED WITH ARGUMENT TOO LARGE FOR ACCURACY ");
      PUT(X); NEW_LINE;
    end if;

    N := INTEGER(X * TWO_OVER_PI);
    XN := CONVERT_TO_FLOAT(N);
    X1 := TRUNCATE(X);
    X2 := X - X1;
    F := ((X1 - XN*C1) + X2) - XN*C2;
    if ABS(F) < EPSILON  then
      RESULT := F;
    else
      G := F * F;
      RESULT := R(G);
    end if;
    if N mod 2 = 0  then
      return RESULT;
    else
      return -1.0/RESULT;
    end if;
  end TAN;

  function COT(X : FLOAT) return FLOAT is
    SGN, Y : FLOAT;
    N : INTEGER;
    XN : FLOAT;
    F, G, X1, X2 : FLOAT;
    RESULT : FLOAT;


    YMAX : FLOAT := FLOAT(INTEGER(PI * TWO**(IT/2))) /2.0;
    BETA : FLOAT := CONVERT_TO_FLOAT(IBETA);
    EPSILON : FLOAT := BETA ** (-IT/2);
    EPSILON1 : FLOAT :=  1.0/XMAX;

    C1 : constant FLOAT :=  8#1.444#;
    C2 : constant FLOAT :=  4.8382_67948_97E-4;

    function R(G : FLOAT) return FLOAT is
      P0 : constant FLOAT :=  1.0;
      P1 : constant FLOAT := -0.11136_14403_566;
      P2 : constant FLOAT :=  0.10751_54738_488E-2;
      Q0 : constant FLOAT :=  1.0;
      Q1 : constant FLOAT := -0.44469_47720_281;
      Q2 : constant FLOAT :=  0.15973_39213_300E-1;
    begin
      return ((P2*G + P1)*G*F + F) / (((Q2*G + Q1)*G +0.5) + 0.5);
    end R;

  begin
    Y := ABS(X);
    if Y < EPSILON1  then
      NEW_LINE;
      PUT(" COT CALLED WITH ARGUMENT TOO NEAR ZERO ");
      PUT(X); NEW_LINE;
      if X < 0.0  then
        return -XMAX;
      else
        return XMAX;
      end if;
    end if;
    if Y > YMAX  then
      NEW_LINE;
      PUT(" COT CALLED WITH ARGUMENT TOO LARGE FOR ACCURACY ");
      PUT(X); NEW_LINE;
    end if;

    N := INTEGER(X * TWO_OVER_PI);
    XN := CONVERT_TO_FLOAT(N);
    X1 := TRUNCATE(X);
    X2 := X - X1;
    F := ((X1 - XN*C1) + X2) - XN*C2;
    if ABS(F) < EPSILON  then
      RESULT := F;
    else
      G := F * F;
      RESULT := R(G);
    end if;
    if N mod 2 /= 0  then
      return -RESULT;
    else
      return 1.0/RESULT;
    end if;
  end COT;


  function ASIN(X : FLOAT) return FLOAT is
    G, Y : FLOAT;
    RESULT : FLOAT;
    BETA : FLOAT := CONVERT_TO_FLOAT(IBETA);
    EPSILON : FLOAT := BETA ** (-IT/2);

    function R(G : FLOAT) return FLOAT is
    P1 : constant FLOAT := -0.27516_55529_0596E1;
    P2 : constant FLOAT :=  0.29058_76237_4859E1;
    P3 : constant FLOAT := -0.59450_14419_3246;
    Q0 : constant FLOAT := -0.16509_93320_2424E2;
    Q1 : constant FLOAT :=  0.24864_72896_9164E2;
    Q2 : constant FLOAT := -0.10333_86707_2113E2;
    Q3 : constant FLOAT :=  1.0;
    begin
      return (((P3*G + P2)*G + P1)*G) / (((G + Q2)*G + Q1)*G + Q0);
    end R;

  begin
    Y := ABS(X);

    if Y > HALF  then
      if Y > 1.0  then
        NEW_LINE; PUT(" ASIN CALLED FOR "); PUT(X);
            PUT(" (> 1)  TRUNCATED TO 1, CONTINUED"); NEW_LINE;
        Y := 1.0;
      end if;
      G := ((0.5 - Y) + 0.5) / 2.0;
      Y := -2.0 * SQRT(G);
      RESULT := Y + Y * R(G);
      RESULT := (PI_OVER_FOUR + RESULT) + PI_OVER_FOUR;
    else
      if Y < EPSILON  then
        RESULT := Y;
      else
        G := Y * Y;
        RESULT := Y + Y * R(G);
      end if;
    end if;
    if X < 0.0  then
      RESULT := -RESULT;
    end if;

    return RESULT;
  end ASIN;

  function ACOS(X : FLOAT) return FLOAT is
    G, Y : FLOAT;
    RESULT : FLOAT;
    BETA : FLOAT := CONVERT_TO_FLOAT(IBETA);
    EPSILON : FLOAT := BETA ** (-IT/2);

    function R(G : FLOAT) return FLOAT is
    P1 : constant FLOAT := -0.27516_55529_0596E1;
    P2 : constant FLOAT :=  0.29058_76237_4859E1;
    P3 : constant FLOAT := -0.59450_14419_3246;
    Q0 : constant FLOAT := -0.16509_93320_2424E2;
    Q1 : constant FLOAT :=  0.24864_72896_9164E2;
    Q2 : constant FLOAT := -0.10333_86707_2113E2;
    Q3 : constant FLOAT :=  1.0;
    begin
      return (((P3*G + P2)*G + P1)*G) / (((G + Q2)*G + Q1)*G + Q0);
    end R;

  begin
    Y := ABS(X);

    if Y > HALF  then
      if Y > 1.0  then
        NEW_LINE; PUT(" ACOS CALLED FOR "); PUT(X);
            PUT(" (> 1)  TRUNCATED TO 1, CONTINUED"); NEW_LINE;
        Y := 1.0;
      end if;
      G := ((0.5 - Y) + 0.5) / 2.0;
      Y := -2.0 * SQRT(G);
      RESULT := Y + Y * R(G);
      if X < 0.0  then
        RESULT := (PI_OVER_TWO + RESULT) + PI_OVER_TWO;
      else
        RESULT := -RESULT;
      end if;

    else
      if Y < EPSILON  then
        RESULT := Y;
      else
        G := Y * Y;
        RESULT := Y + Y * R(G);
      end if;
      if X < 0.0  then
        RESULT := (PI_OVER_FOUR + RESULT) + PI_OVER_FOUR;
      else
        RESULT := (PI_OVER_FOUR - RESULT) + PI_OVER_FOUR;
      end if;
    end if;

    return RESULT;
  end ACOS;


  function ATAN(X : FLOAT) return FLOAT is
    F, G : FLOAT;
    subtype REGION is INTEGER range 0..3;    --  ##########
    N : REGION;
    RESULT : FLOAT;

    BETA : FLOAT := CONVERT_TO_FLOAT(IBETA);

    EPSILON : FLOAT := BETA ** (-IT/2);

    SQRT_3           : constant FLOAT :=  1.73205_08075_68877_29353;
    SQRT_3_MINUS_1   : constant FLOAT :=  0.73205_08075_68877_29353;
    TWO_MINUS_SQRT_3 : constant FLOAT :=  0.26794_91924_31122_70647;

    function R(G : FLOAT) return FLOAT is
      P0 : constant FLOAT := -0.14400_83448_74E1;
      P1 : constant FLOAT := -0.72002_68488_98;
      Q0 : constant FLOAT :=  0.43202_50389_19E1;
      Q1 : constant FLOAT :=  0.47522_25845_99E1;
      Q2 : constant FLOAT :=  1.0;
    begin
      return ((P1*G + P0)*G) / ((G + Q1)*G + Q0);
    end R;

  begin
    F := ABS(X);

    if F > 1.0  then
      F := 1.0 / F;
      N := 2;
    else
      N := 0;
    end if;

    if F > TWO_MINUS_SQRT_3  then
      F := (((SQRT_3_MINUS_1 * F - 0.5) - 0.5) + F) / (SQRT_3 + F);
      N := N + 1;
    end if;

    if ABS(F) < EPSILON  then
      RESULT := F;
    else
      G := F * F;
      RESULT := F + F * R(G);
    end if;

    if N > 1  then
      RESULT := - RESULT;
    end if;

    case N is
    when 0  =>
      RESULT := RESULT;
    when 1  =>
      RESULT := PI_OVER_SIX + RESULT;
    when 2  =>
      RESULT := PI_OVER_TWO + RESULT;
    when 3  =>
      RESULT := PI_OVER_THREE + RESULT;
    end case;

    if X < 0.0  then
      RESULT := - RESULT;
    end if;

    return RESULT;

  end ATAN;



  function ATAN2(V, U : FLOAT) return FLOAT is
    X, RESULT : FLOAT;

  begin

    if U = 0.0  then
      if V = 0.0  then
        RESULT := 0.0;
        NEW_LINE;
        PUT(" ATAN2 CALLED WITH 0/0   RETURNED "); PUT(RESULT);
        NEW_LINE;
      elsif V > 0.0  then
        RESULT := PI_OVER_TWO;
      else
        RESULT := - PI_OVER_TWO;
      end if;

    else
      X := ABS(V/U);
      --  If underflow or overflow is detected, go to the exception
      RESULT := ATAN(X);
      if U < 0.0  then
        RESULT := PI - RESULT;
      end if;
      if V < 0.0  then
        RESULT := - RESULT;
      end if;
    end if;
    return RESULT;
  exception
  when NUMERIC_ERROR  =>
    if ABS(V) > ABS(U)  then
      RESULT := PI_OVER_TWO;
      if V < 0.0  then
        RESULT := - RESULT;
      end if;
    else
      RESULT := 0.0;
      if U < 0.0  then
        RESULT := PI - RESULT;
      end if;
    end if;
    return RESULT;
  end ATAN2;


  function SINH(X : FLOAT) return FLOAT is
    G, W, Y, Z : FLOAT;
    RESULT : FLOAT;
    BETA : FLOAT := CONVERT_TO_FLOAT(IBETA);
    EPSILON : FLOAT := BETA ** (-IT/2);

    YBAR : FLOAT := EXP_LARGE;
    LN_V : FLOAT := 8#0.542714#;
    V_OVER_2_MINUS_1 : FLOAT :=  0.13830_27787_96019_02638E-4;
    WMAX : FLOAT := YBAR - LN_V + 0.69;

    function R(G : FLOAT) return FLOAT is
    P0 : constant FLOAT :=  0.10622_28883_7151E4;
    P1 : constant FLOAT :=  0.31359_75645_6058E2;
    P2 : constant FLOAT :=  0.34364_14035_8506;
    Q0 : constant FLOAT :=  0.63733_73302_1822E4;
    Q1 : constant FLOAT := -0.13051_01250_9199E3;
    Q2 : constant FLOAT :=  1.0;
    begin
      return (((P2*G + P1)*G + P0)*G) / ((G + Q1)*G + Q0);
    end R;

  begin
    Y := ABS(X);

    if Y <= 1.0  then
      if Y < EPSILON  then
        RESULT := X;
      else
        G := X * X;
        RESULT := X + X * R(G);
      end if;

    else
      if Y <= YBAR  then
        Z := EXP(Y);
        RESULT := (Z - 1.0/Z) / 2.0;
      else
        W := Y - LN_V;
        if W > WMAX  then
          NEW_LINE;
          PUT(" SINH CALLED WITH TOO LARGE ARGUMENT  "); PUT(X);
          PUT(" RETURN BIG"); NEW_LINE;
          W := WMAX;
        end if;
        Z := EXP(W);
        RESULT := Z + V_OVER_2_MINUS_1 * Z;
      end if;
      if X < 0.0  then
        RESULT := -RESULT;
      end if;

    end if;
    return RESULT;
  end SINH;


  function COSH(X : FLOAT) return FLOAT is
    G, W, Y, Z : FLOAT;
    RESULT : FLOAT;
    BETA : FLOAT := CONVERT_TO_FLOAT(IBETA);
    EPSILON : FLOAT := BETA ** (-IT/2);

    YBAR : FLOAT := EXP_LARGE;
    LN_V : FLOAT := 8#0.542714#;
    V_OVER_2_MINUS_1 : FLOAT :=  0.13830_27787_96019_02638E-4;
    WMAX : FLOAT := YBAR - LN_V + 0.69;

    function R(G : FLOAT) return FLOAT is
    P0 : constant FLOAT :=  0.10622_28883_7151E4;
    P1 : constant FLOAT :=  0.31359_75645_6058E2;
    P2 : constant FLOAT :=  0.34364_14035_8506;
    Q0 : constant FLOAT :=  0.63733_73302_1822E4;
    Q1 : constant FLOAT := -0.13051_01250_9199E3;
    Q2 : constant FLOAT :=  1.0;
    begin
      return (((P2*G + P1)*G + P0)*G) / ((G + Q1)*G + Q0);
    end R;

  begin
    Y := ABS(X);

    if Y <= YBAR  then
      Z := EXP(Y);
      RESULT := (Z + 1.0/Z) / 2.0;
    else
      W := Y - LN_V;
      if W > WMAX  then
        NEW_LINE;
        PUT(" COSH CALLED WITH TOO LARGE ARGUMENT  "); PUT(X);
        PUT(" RETURN BIG"); NEW_LINE;
        W := WMAX;
      end if;
      Z := EXP(W);
      RESULT := Z + V_OVER_2_MINUS_1 * Z;
    end if;

    return RESULT;
  end COSH;


  function TANH(X : FLOAT) return FLOAT is
    G, W, Y, Z : FLOAT;
    RESULT : FLOAT;
    BETA : FLOAT := CONVERT_TO_FLOAT(IBETA);
    EPSILON : FLOAT := BETA ** (-IT/2);

    XBIG : FLOAT := (LOG(2.0) + CONVERT_TO_FLOAT(IT + 1) * LOG(BETA))/2.0;
    LN_3_OVER_2 : FLOAT :=  0.54930_61443_34054_84570;

    function R(G : FLOAT) return FLOAT is
    P0 : constant FLOAT := -0.21063_95800_0245E2;
    P1 : constant FLOAT := -0.93363_47565_2401;
    Q0 : constant FLOAT :=  0.63191_87401_5582E2;
    Q1 : constant FLOAT :=  0.28077_65347_0471E2;
    Q2 : constant FLOAT :=  1.0;
    begin
      return ((P1*G + P0)*G) / ((G + Q1)*G + Q0);
    end R;

  begin
    Y := ABS(X);

    if Y > XBIG  then
      RESULT := 1.0;
    else
      if Y > LN_3_OVER_2  then
        RESULT := 0.5 - 1.0 / (EXP(Y + Y) + 1.0);
        RESULT := RESULT + RESULT;
      else
        if Y < EPSILON  then
          RESULT := Y;
        else
          G := Y * Y;
          RESULT := Y + Y * R(G);
        end if;
      end if;
    end if;
    if X < 0.0  then
      RESULT := - RESULT;
    end if;

    return RESULT;
  end TANH;


begin
  null;
end TRIG_LIB;