--::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
--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

   package INTEGER_IO is
      new TEXT_IO.INTEGER_IO(INTEGER);

      package FLOAT_IO is
         new TEXT_IO.FLOAT_IO(FLOAT); 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 : in out EXPONENT_TYPE;
                           F : in out MANTISSA_TYPE);

         procedure REFLOAT(N : in EXPONENT_TYPE;
                           F : in MANTISSA_TYPE;
                           X : in 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 : in out EXPONENT_TYPE;
                        F : in 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 : in 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;