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

Text File | 1988-05-03 | 75.6 KB | 1,857 lines

with Text_Io;
generic
    type Floating is digits <>;
package Generic_Math_Functions is
    use Text_Io;
    Pi        : constant := 3.14159_26535_89793_23846_26433_83279_50288_41972;
    E         : constant := 2.71828_18284_59045_23536_02874_71352_66249_77572;
    Log_Of_2  : constant := 0.69314_71805_59945_30941_72321_24158_17656_80755;
    Log_Of_10 : constant := 2.30258_50929_94045_68401_77914_54684_36420_76011;
    type Error_Response is
        (Return_Default,
         Return_Default_And_Log,
         Raise_Math_Function_Exception,
         Propagate_Exception);
    What_To_Do_When_There_Is_An_Error :
        Error_Response := Return_Default_And_Log;
    Math_Function_Exception           : exception;
    Error_Log : Text_Io.File_Type;
    function Sign (X, Y : Floating) return Floating;
--  Returns the value of X with the sign of Y
    function Max (X, Y : Floating) return Floating;
--  Returns the algebraicly larger of X and Y
    function Min (X, Y : Floating) return Floating;
--  Returns the algebraicly smaller of X and Y
    function Truncate (X : Floating) return Floating;
--  Returns the floating value of the integer no larger than X
--  Truncates toward zero
    function Round (X : Floating) return Floating;
--  Returns the floating value of the integer nearest X
    procedure Set_Ran_Key (K : in Floating := Floating (0.0));
--  Can reset the random number generator
    function Ran return Floating;
--  A random number between zero and one
    function Sqrt (X : Floating) return Floating;
    function Cbrt (X : Floating) return Floating;
    function Log (X : Floating) return Floating;
    function Log10 (X : Floating) return Floating;
    function Exp (X : Floating) return Floating;
    function "**" (X, Y : Floating) return Floating;
    function Sin (X : Floating) return Floating;
    function Cos (X : Floating) return Floating;
    function Tan (X : Floating) return Floating;
    function Cot (X : Floating) return Floating;
    function Asin (X : Floating) return Floating;
    function Acos (X : Floating) return Floating;
    function Atan (X : Floating) return Floating;
    function Atan2 (V, U : Floating) return Floating;
    function Sinh (X : Floating) return Floating;
    function Cosh (X : Floating) return Floating;
    function Tanh (X : Floating) return Floating;
end Generic_Math_Functions;
with Text_Io;
package body Generic_Math_Functions is
--  The purpose of this package is to provide a portable Ada facility
--  for new installations, considerable optimization could be done
--  for particular machines or if fewer functions were required or if
--  error reporting were eliminated or exceptions suppressed
--  The routines herein are fairly machine-independent, safe, and
--  accurate to about digits 18
--  Each routine is stand-alone from the others - they do not share fragments
--  but LOG10 uses LOG, ATAN2 uses TAN, and SINH, COSH, TANH uses EXP
--  and LOG and EXP are used in constants
--  Test routines are also available
--  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
--      16 JULY 1982       W A WHITAKER  AFATL EGLIN AFB FL 32542
--                         T C EICHOLTZ  USAFA
--      30  APRIL1986      UPDATED FOR MODERN COMPILERS - WHITAKER
    package Integer_Io is new Text_Io.Integer_Io (Integer);
    package Floating_Io is new Text_Io.Float_Io (Floating);
    package Boolean_Io is new Text_Io.Enumeration_Io (Boolean);
    use Text_Io;
--  The following constants are used in the functions
    One           : constant Floating := Floating (1.0);
    Zero          : constant Floating := One - One;
    Half          : constant Floating := One / (One + One);
    One_Over_Pi   : constant Floating := One / Pi;
    Two_Over_Pi   : constant Floating := (One + One) / Pi;
    Pi_Over_Two   : constant Floating := Pi / (One + One);
    Pi_Over_Three : constant Floating := Pi / (One + One + One);
    Pi_Over_Four  : constant Floating := Pi / (One + One + One + One);
    Beta          : Floating;
    Epsilon       : Floating;
    Ymax          : Floating;
    Exp_Large     : Floating;
    Exp_Small     : Floating;
--  These are included so that the package can be stand-alone and still
--  put out diagnostics on errors
    procedure Put (File : in File_Type; Item : in Integer;
                   Width : in Field := Integer_Io.Default_Width;
                   Base : in Number_Base := Integer_Io.Default_Base)
                   renames Integer_Io.Put;
    procedure Put (Item : in Integer;
                   Width : in Field := Integer_Io.Default_Width;
                   Base : in Number_Base := Integer_Io.Default_Base)
                   renames Integer_Io.Put;
    procedure Put (File : in File_Type; Item : in Floating;
                   Fore : in Field := Floating_Io.Default_Fore;
                   Aft : in Field := Floating_Io.Default_Aft;
                   Exp : in Field := Floating_Io.Default_Exp)
                   renames Floating_Io.Put;
    procedure Put (Item : in Floating;
                   Fore : in Field := Floating_Io.Default_Fore;
                   Aft : in Field := Floating_Io.Default_Aft;
                   Exp : in Field := Floating_Io.Default_Exp)
                   renames Floating_Io.Put;
    package Floating_Characteristics is
--  This package is a floating mantissa definition of a binary floating
--  The parameters are obtained by initializing on the actual hardware
--  Otherwise the parameters could be set in the spec if known
--  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 and are reproduced there
--  Many of these could be in the form of attributes for hardware float types
--  but some are not easily available for the machine hardware base
--  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_floating 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 + floating(IBETA) ** MACHEP /= 1.0
--  except that MACHEP is bounded below by -(IT + 3)
        Negep  : Integer;
--  The largest negative integer such that
--    1.0 -0 floating(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
--  floating(IBETA) ** MINEXP is a positive floating-point number
        Maxexp : Integer;
--  The largest positive exponent for a finite floating-point number
        Eps    : Floating;
--  The smallest positive floating-point number such that
--                              1.0 + EPS /= 1.0
--  In particular, if IBETA = 2 or IRND = 0,
--  EPS = floating(IBETA) ** MACHEP
--  Otherwise, EPS = (floating(IBETA) ** MACHEP) / 2
        Epsneg : Floating;
--  A small positive floating-poin number such that 1.0-EPSNEG /= 1.0
        Xmin   : Floating;
--  The smallest non-vanishing floating-point power of the radix
--  In particular, XMIN = floating(IBETA) ** MINEXP
        Xmax   : Floating;
--  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, but is floating here
--  Most efficient implementation may require details of the machine hardware
--  In this version the simplest representation is used
--  The mantissa is extracted into a floating and uses the predefined operations
        type Exponent_Type is new Integer;
        type Mantissa_Type is new Floating;
        --  Now we do something strange
--  Since we may 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 FLOATING/2 will fail
--  We define a type-dependent factor that will work
        Mantissa_Divide_By_2 : constant Mantissa_Type := 1.0 / 2.0;
        Mantissa_Divide_By_3 : constant Mantissa_Type := 1.0 / 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 Floating;
                           N : in out Exponent_Type; F : in out Mantissa_Type);
        procedure Refloat (N : in Exponent_Type; F : in Mantissa_Type;
                           X : in out Floating);
    end Floating_Characteristics;
    package body Floating_Characteristics is
--  This package is a floating mantissa definition of a binary floating
        package Exponent_Io is new Text_Io.Integer_Io (Exponent_Type);
        package Mantissa_Io is new Text_Io.Float_Io (Mantissa_Type);
--  Again we define PUT only so we can do diagnostics in DEFLOAT and REFLOAT
        procedure Put (File : in File_Type; Item : in Exponent_Type;
                       Width : in Field := Integer_Io.Default_Width;
                       Base : in Number_Base := Integer_Io.Default_Base)
                       renames Exponent_Io.Put;
        procedure Put (Item : in Exponent_Type;
                       Width : in Field := Integer_Io.Default_Width;
                       Base : in Number_Base := Integer_Io.Default_Base)
                       renames Exponent_Io.Put;
        procedure Put (File : in File_Type; Item : in Mantissa_Type;
                       Fore : in Field := Floating_Io.Default_Fore;
                       Aft : in Field := Floating_Io.Default_Aft;
                       Exp : in Field := Floating_Io.Default_Exp)
                       renames Mantissa_Io.Put;
        procedure Put (Item : in Mantissa_Type;
                       Fore : in Field := Floating_Io.Default_Fore;
                       Aft : in Field := Floating_Io.Default_Aft;
                       Exp : in Field := Floating_Io.Default_Exp)
                       renames Mantissa_Io.Put;
        procedure Defloat (X : in Floating;
                           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         : Floating;
        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 / Floating (2.0) ** Integer (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 := Floating (2.0) ** Integer (M);
                    W := Y / Z;
                    if W >= 0.5 then
                        Y := W;
                        N := N + M;
                    end if;
                end loop;
            end if;
        --  And just to clear up any loose ends from biased exponents
            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 =>
                if What_To_Do_When_There_Is_An_Error =
                    Return_Default_And_Log then
                    Put (Error_Log,
                         "Exception raised in FLOATING_CHARACTERISTICS.DEFLOAT  X = ");
                    Put (Error_Log, X);
                    New_Line (Error_Log);
                elsif What_To_Do_When_There_Is_An_Error = Propagate_Exception then
                    raise Math_Function_Exception;
                end if;
                N := 0;
                F := 0.0;
                return;
        end Defloat;
        procedure Refloat (N : in Exponent_Type; F : in Mantissa_Type;
                           X : in out Floating) is
--  Again a brute force method - but portable
--  Watch out near MAXEXP
            M : Integer;
            Y : Floating;
        begin
            if F = 0.0 then
                X := Zero;
                return;
            end if;
            M := Integer (N);
            Y := abs (Floating (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;
        exception
            when others =>
                if What_To_Do_When_There_Is_An_Error = Return_Default_And_Log then
                    Put (Error_Log,
                         "Exception raised in FLOATING_CHARACTERISTICS.REFLOAT  N = ");
                    Put (Error_Log, N);
                    Put (Error_Log, "  F = ");
                    Put (Error_Log, F);
                    New_Line (Error_Log);
                elsif What_To_Do_When_There_Is_An_Error = Propagate_Exception then
                    raise Math_Function_Exception;
                end if;
                X := Zero;
                return;
        end Refloat;
        procedure Initialize is
--  This initialization is the MACAR routine of Cody and Waite Appendix B.
            A, B, Y, Z           : Floating;
            I, K, Mx, Iz         : Integer;
            Beta, Betam1, Betain : Floating;
        begin
            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  := Floating (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
                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;
            Find_Iexp:
                declare
                    Y : Floating := 1.0;
                    A : Floating := Betain;
                    I : Integer := 0;
                begin
                    loop
                        Y := A * A;
                        exit when Y = 0.0;
                        I := I + 1;
                        A := Y;
                    end loop;
                    Iexp := I;
                end Find_Iexp;
            Find_Smallest:
                declare
                    A, Y : Floating := 1.0;
                    I    : Integer := 0;
                begin
                    loop
                        Y := A / Beta;
                        exit when Y = 0.0;
                        I := I - 1;
                        A := Y;
                    end loop;
                    Minexp := I;
                    Xmin   := A;
                end Find_Smallest;
            Find_Largest:
                declare
                    A, Y : Floating := 1.0;
                    I    : Integer := 1;
                begin
                    loop
                        Y := A * Beta;
                        I := I + 1;
                        A := Y;
                    end loop;
                exception
                    when others =>
                        Maxexp := I;
                        Xmax   := A * ((1.0 - Epsneg) * Beta);
                end Find_Largest;
        end Initialize;
    begin
        Initialize;
    end Floating_Characteristics;
    use Floating_Characteristics;
    function Sign (X, Y : Floating) return Floating 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 : Floating) return Floating is
    begin
        if X >= Y then
            return X;
        else
            return Y;
        end if;
    end Max;
    function Min (X, Y : Floating) return Floating is
    begin
        if X <= Y then
            return X;
        else
            return Y;
        end if;
    end Min;
    function Truncate (X : Floating) return Floating is
    --  Optimum code depends on how the system rounds at exact halves
        Z : Floating := Floating (Integer (X));
    begin
        if X = Z then
            return Z;
        else
            if X > Zero then
        --  For instance, you can get into trouble when 1.0E-20 is overwhelmed
        --  by HALF and X - HALF is -HALF and that is "truncated" to -1
        --  so a simple floating (Integer (X - Half)) does not work in all cases
                if X < Z then
                    return Z - One;
                else
                    return Z;
                end if;
            elsif X = Zero then
                return Zero;
            else
                if X < Z then
                    return Z;
                else
                    return Z + One;
                end if;
            end if;
        end if;
    end Truncate;
    function Round (X : Floating) return Floating is
    begin
        return Floating (Integer (X));
    end Round;
    package Key is
        X : Integer := 10_001;
        Y : Integer := 20_001;
        Z : Integer := 30_001;
    end Key;
    procedure Set_Ran_Key (K : in Floating := Floating (0.0)) is
        Dummy : Floating := 0.0;
    begin
        if K = Zero then
            Key.X := 10_001;
            Key.Y := 20_001;
            Key.Z := 30_001;
        else
            for I in 1..Integer (K - Floating (Integer (abs (K) / 31_247.0) *
                  31_247.0)) mod
                  31_247 loop
                Dummy := Ran;
            end loop;
        end if;
    end Set_Ran_Key;
    function Ran return Floating is
--  This uses a portable algorithm and is included at this point so that
--  an implementation could take advantage of unique machine characteristics
--  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 : Floating;
    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;
        W := Floating (Key.X) / 30269.0
             + Floating (Key.Y) / 30307.0
             + Floating (Key.Z) / 30323.0;
        return W - Floating (Integer (W - 0.5));
    end Ran;
    function Sqrt (X : Floating) return Floating is
        M, N    : Exponent_Type;
        F, Y    : Mantissa_Type;
        Result  : Floating;
        Sqrt_C1 : constant Mantissa_Type := 0.41731;
        Sqrt_C2 : constant Mantissa_Type := 0.59016;
        Sqrt_C3 : constant Mantissa_Type := 0.70710_67811_86547_52440;
    begin
        if X = Zero then
            Result := Zero;
            return Result;
        elsif X = One then
            Result := One;     --  To get exact SQRT(1.0)
            return Result;
        elsif X < Zero then
            case What_To_Do_When_There_Is_An_Error is
                when Return_Default =>
                    null;
                when Return_Default_And_Log =>
                    New_Line (Error_Log);
                    Put (Error_Log, "CALLED SQRT FOR NEGATIVE ARGUMENT   ");
                    Put (Error_Log, X);
                    Put (Error_Log, "   USED ABSOLUTE VALUE");
                when Raise_Math_Function_Exception =>
                    raise Math_Function_Exception;
                when Propagate_Exception =>
                    null;
            end case;
            Result := Sqrt (abs (X));
            return Result;
        else
            Defloat (X, N, F);
            Y := Sqrt_C1 + Mantissa_Type (Sqrt_C2 * F);
            if N mod 2 = 1 then
                Y := Y * Sqrt_C3;
                F := F * Mantissa_Divide_By_2;
                N := N + 1;
            end if;
            M := N / 2;
--      Since always looping to 3, unroll the loop
            Y := Mantissa_Type (Floating (Y) + Floating (F) / Floating (Y)) *
                Mantissa_Divide_By_2;
            Y := Mantissa_Type (Floating (Y) + Floating (F) / Floating (Y)) *
                Mantissa_Divide_By_2;
            Y := Mantissa_Type (Floating (Y) + Floating (F) / Floating (Y)) *
                Mantissa_Divide_By_2;
            Refloat (M, Y, Result);
            return Result;
        end if;
    exception
        when others =>
            case What_To_Do_When_There_Is_An_Error is
                when Return_Default =>
                    null;
                when Return_Default_And_Log =>
                    New_Line (Error_Log);
                    Put (Error_Log, " EXCEPTION IN SQRT, X = ");
                    Put (Error_Log, X);
                    Put_Line (Error_Log, "  RETURNED 1.0");
                when Raise_Math_Function_Exception =>
                    raise Math_Function_Exception;
                when Propagate_Exception =>
                    raise;
            end case;
            return One;
    end Sqrt;
    function Cbrt (X : Floating) return Floating is
        M, N    : Exponent_Type;
        F, Y    : Mantissa_Type;
        Result  : Floating;
        Cbrt_C1 : constant Mantissa_Type := 0.5874009;
        Cbrt_C2 : constant Mantissa_Type := 0.4125990;
        Cbrt_C3 : constant Mantissa_Type := 0.62996_05249;
        Cbrt_C4 : constant 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);
            case (N mod 3) is
                when 0 =>
                    null;
                when 1 =>
                    Y := Mantissa_Type (Cbrt_C3 * Y);
                    F := F / 4.0;
                    N := N + 2;
                when 2 =>
                    Y := Mantissa_Type (Cbrt_C4 * Y);
                    F := F / 2.0;
                    N := N + 1;
                when others =>
                    null;
            end case;
            M := N / 3;
            Y := Y
                 - (Y * Mantissa_Divide_By_3
                 - Mantissa_Type
                 ((F * Mantissa_Divide_By_3) / Mantissa_Type (Y *
                Y)));
            Y := Y
                 - (Y * Mantissa_Divide_By_3
                 - Mantissa_Type
                 ((F * Mantissa_Divide_By_3) / Mantissa_Type (Y *
                Y)));
            Y := Y
                 - (Y * Mantissa_Divide_By_3
                 - Mantissa_Type
                 ((F * Mantissa_Divide_By_3) / Mantissa_Type (Y *
                Y)));
            Y := Y
                 - (Y * Mantissa_Divide_By_3
                 - Mantissa_Type
                 ((F * Mantissa_Divide_By_3) / Mantissa_Type (Y *
                Y)));
            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 (Error_Log);
            Put (Error_Log, "EXCEPTION IN CBRT, X = ");
            Put (Error_Log, X);
            Put (Error_Log, "  RETURNED  ");
            Put (Error_Log, Result);
            New_Line (Error_Log);
            return Result;
    end Cbrt;
    function Log (X : Floating) return Floating is
        Result        : Floating;
        N             : Exponent_Type;
        Xn            : Floating;
        Y             : Floating;
        F             : Mantissa_Type;
        Z, Zden, Znum : Mantissa_Type;
        C0            : constant Mantissa_Type := 0.70710_67811_86547_52440;
        C1            : constant Floating := 8#0.543#;
        C2            : constant Floating := - 2.12194_44005_46905_82767_9E-4;
        function R (Z : Mantissa_Type) return Mantissa_Type is
            W : Mantissa_Type := Z * Z;
        begin
            if Floating'Base'Digits < 10 then
                declare
                    A0 : constant Mantissa_Type := - 0.46490_62303_464e+0;
                    A1 : constant Mantissa_Type := 0.13600_95468_621e-1;
                    B0 : constant Mantissa_Type := - 0.55788_73750_242e+1;
                    B1 : constant Mantissa_Type := 0.10000_00000_000e+1;
                begin
                    return
                        Z + Z * W * (A0 + A1 * W) / ((B0 + W));
                end;
            else
                declare
                    A0 : constant Mantissa_Type :=
                        - 0.64124_94342_37455_81147e+2;
                    A1 : constant Mantissa_Type :=
                        0.16383_94356_30215_34222e+2;
                    A2 : constant Mantissa_Type :=
                        - 0.78956_11288_74912_57267e+0;
                    B0 : constant Mantissa_Type :=
                        - 0.76949_93210_84948_79777e+3;
                    B1 : constant Mantissa_Type :=
                        0.31203_22209_19245_32844e+3;
                    B2 : constant Mantissa_Type :=
                        - 0.35667_97773_90346_46171e+2;
                    B3 : constant Mantissa_Type :=
                        0.10000_00000_00000_00000e+1;
                begin
                    return
                        Z + Z * W * (A0 + (A1 + A2 * W) * W) /
                        (B0 + (B1 + (B2 + W) * W) * W);
                end;
            end if;
        end R;
    begin
        if X < Zero then
            New_Line (Error_Log);
            Put (Error_Log, "CALLED LOG FOR NEGATIVE ");
            Put (Error_Log, X);
            Put (Error_Log, "   USE ABS => ");
            Result := Log (abs (X));
            Put (Error_Log, Result);
            New_Line (Error_Log);
        elsif X = Zero then
            New_Line (Error_Log);
            Put (Error_Log, "CALLED LOG FOR ZERO ARGUMENT, RETURNED ");
            Result := - Xmax;      --  SUPPOSED TO BE -LARGE
            Put (Error_Log, Result);
            New_Line (Error_Log);
        else
            Defloat (X, N, F);
            Znum := F - Mantissa_Half;
            if F > C0 then
                Znum := Znum - Mantissa_Half;
                Zden := F * Mantissa_Divide_By_2 + Mantissa_Half;
            else
                N    := N - 1;
                Zden := Znum * Mantissa_Divide_By_2 + Mantissa_Half;
            end if;
            Z := Mantissa_Type (Znum / Zden);
            if N = 0 then
                Result := Floating (R (Z));
            else
                Xn     := Floating (N);
                Result := (Xn * C2 + Floating (R (Z))) + Xn * C1;
            end if;
        end if;
        return Result;
    exception
        when others =>
           New_Line (Error_Log);
            Put (Error_Log, " EXCEPTION IN LOG, X = ");
            Put (Error_Log, X);
            Put (Error_Log, "  RETURNED 0.0");
            New_Line (Error_Log);
            return Zero;
    end Log;
    function Log10 (X : Floating) return Floating is
        Log_Base_10_Of_E : constant Floating := 0.43429_44819_03251_82765;
    begin
        return Log_Base_10_Of_E * Log (X);
    end Log10;
    function Exp (X : Floating) return Floating is
        Result         : Floating;
        N              : Exponent_Type;
        Xg, Xn, X1, X2 : Floating;
        F, G           : Mantissa_Type;
        One_Over_Log_2 : constant Floating := 1.4426_95040_88896_34074;
        C1             : constant Floating := 0.69335_9375;
        C2             : constant Floating := - 2.1219_44400_54690_58277E-4;
        function R (G : Mantissa_Type) return Mantissa_Type is
            Z, Gp, Q : Mantissa_Type;
        begin
            Z := Mantissa_Type (G * G);
            if Floating'Base'Digits < 9 then
     --  Good for IT <30
                declare
                    P0 : constant Mantissa_Type := 0.24999_99995_0e+0;
                    P1 : constant Mantissa_Type := 0.41602_88626_8e-2;
                    Q0 : constant Mantissa_Type := 0.5;
                    Q1 : constant Mantissa_Type := 0.49987_17877_8e-1;
                begin
                    Gp := Mantissa_Type ((Mantissa_Type (P1 * Z) + P0) * G);
                    Q  := Mantissa_Type (Q1 * Z) + Q0;
                end;
            else
                              --  Good for IT < 57
                declare
                    P0 : constant Mantissa_Type := 0.24999_99999_99999_993e+0;
                    P1 : constant Mantissa_Type := 0.69436_00015_11792_852e-2;
                    P2 : constant Mantissa_Type := 0.16520_33002_68279_130e-4;
                    Q0 : constant Mantissa_Type := 0.5;
                    Q1 : constant Mantissa_Type := 0.55553_86669_69001_188e-1;
                    Q2 : constant Mantissa_Type := 0.49586_28849_05441_294e-3;
                begin
                    Gp := Mantissa_Type ((Mantissa_Type ((P2 * Z + P1) * Z) +
                        P0) * G);
                    Q  := Mantissa_Type ((Mantissa_Type (Q2 * Z) + Q1) * Z) +
                        Q0;
                end;
            end if;
            return Mantissa_Half + Mantissa_Type (Gp / (Q - Gp));
        end R;
    begin
        if X > Exp_Large then
            New_Line (Error_Log);
            Put (Error_Log,
                 "  EXP CALLED WITH TOO BIG A POSITIVE ARGUMENT, ");
            Put (Error_Log, X);
            New_Line (Error_Log);
            Put (Error_Log, "   RETURNED XMAX");
            New_Line (Error_Log);
            Result := Xmax;
        elsif X < Exp_Small then
            New_Line (Error_Log);
            Put (Error_Log,
                 "  EXP CALLED WITH TOO BIG A NEGATIVE ARGUMENT,  ");
            Put (Error_Log, X);
            Put (Error_Log, "    RETURNED ZERO");
            New_Line (Error_Log);
            Result := Zero;
        elsif abs (X) < Eps then
            Result := One;
        else
            N  := Exponent_Type (X * One_Over_Log_2);
            Xn := Floating (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 (Error_Log);
            Put (Error_Log, " EXCEPTION IN EXP, X = ");
            Put (Error_Log, X);
            Put (Error_Log, "  RETURNED 1.0");
            New_Line (Error_Log);
            return One;
    end Exp;
    function "**" (X, Y : Floating) return Floating is
        M, N : Exponent_Type;
        G : Mantissa_Type;
        P, Temp, Iw1, I : Integer;
        Z, V, R, U1, U2, W1, W2, W3, Y1, Y2 : Mantissa_Type;
        A1p1 : Mantissa_Type;
        W, Result : Floating;
        K : constant := 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));
     --  2**(-i/16) from which A1 and A2 are constructed
        B1 : constant := 8#1.00000_00000_00000_00000_00000_00000#;
        B2 : constant := 8#0.75222_57505_22220_66302_61734_72062#;
        B3 : constant := 8#0.72540_30671_75644_41622_73201_32727#;
        B4 : constant := 8#0.70146_33673_02522_47021_04062_61124#;
        B5 : constant := 8#0.65642_37462_55323_53257_20715_15057#;
        B6 : constant := 8#0.63422_21405_21760_44016_17421_53016#;
        B7 : constant := 8#0.61263_45204_25240_66655_64761_25503#;
        B8 : constant := 8#0.57204_24347_65401_41553_72504_02177#;
        B9 : constant := 8#0.55202_36314_77473_63110_21313_73047#;
        B10 : constant := 8#0.53254_07672_44124_12254_31114_01243#;
        B11 : constant := 8#0.51377_32652_33052_11616_50345_72776#;
        B12 : constant := 8#0.47572_46230_11064_10432_66404_42174#;
        B13 : constant := 8#0.46033_76024_30667_05226_75065_32214#;
        B14 : constant := 8#0.44341_72334_72542_16063_30176_55544#;
        B15 : constant := 8#0.42712_70170_76521_36556_33737_10612#;
        B16 : constant := 8#0.41325_30331_74611_03661_23056_22556#;
        B17 : constant := 8#0.40000_00000_00000_00000_00000_00000#;
        function Reduce (V : Floating) return Mantissa_Type is
        begin
            return Mantissa_Type (Integer (16.0 * V)) * 0.0625;
        end Reduce;
    begin
        if X <= Zero then
            if X < Zero then
                Result := (abs (X)) ** Y;
                New_Line (Error_Log);
                Put (Error_Log, "X**Y CALLED WITH X = ");
                Put (Error_Log, X);
                New_Line (Error_Log);
                Put (Error_Log, "USED ABS, RETURNED ");
                Put (Error_Log, Result);
                New_Line (Error_Log);
            else
                if Y <= Zero then
                    if Y = Zero then
                        Result := Zero;
                    else
                        Result := Xmax;
                    end if;
                    New_Line (Error_Log);
                    Put (Error_Log, "X**Y CALLED WITH X = 0, Y = ");
                    Put (Error_Log, Y);
                    New_Line (Error_Log);
                    Put (Error_Log, "RETURNED ");
                    Put (Error_Log, Result);
                    New_Line (Error_Log);
                else
                    Result := Zero;
                end if;
            end if;
        else
            Defloat (X, M, G);
            if Floating'Base'Digits < 9 then
                declare
                    A1 : array (1..17) of Mantissa_Type :=
                        (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 Mantissa_Type :=
                        (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#);
                begin
                    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));
                end;
            else
                declare
                    A1 : array (1..17) of Mantissa_Type :=
                        (8#1.00000_00000_00000_000#,
                         8#0.75222_57505_22220_662#,
                         8#0.72540_30671_75644_416#,
                         8#0.70146_33673_02522_470#,
                         8#0.65642_37462_55323_532#,
                         8#0.63422_21405_21760_440#,
                         8#0.61263_45204_25240_666#,
                         8#0.57204_24347_65401_414#,
                         8#0.55202_36314_77473_630#,
                         8#0.53254_07672_44124_122#,
                         8#0.51377_32652_33052_116#,
                         8#0.47572_46230_11064_104#,
                         8#0.46033_76024_30667_052#,
                         8#0.44341_72334_72542_160#,
                         8#0.42712_70170_76521_364#,
                         8#0.41325_30331_74611_036#,
                         8#0.40000_00000_00000_000#);
                    A2 : array (1..8) of Mantissa_Type :=
                        (8#0.00000_00000_00000_00102_61734_72062#,
                         8#0.00000_00000_00000_00021_04062_61124#,
                         8#0.00000_00000_00000_00016_17421_53016#,
                         8#0.00000_00000_00000_00153_72504_02177#,
                         8#0.00000_00000_00000_00054_31114_01243#,
                         8#0.00000_00000_00000_00032_66404_42174#,
                         8#0.00000_00000_00000_00063_30176_55544#,
                         8#0.00000_00000_00000_00061_23056_22556#);
                begin
                    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));
                end;
            end if;
            Z := Z + Z;
            V := Z * Z;
            if Floating'Base'Digits < 11 then
                declare
                    P1 : constant := 0.83333_32862_45E-1;
                    P2 : constant := 0.12506_48500_52E-1;
                begin
                    R := (P2 * V + P1) * V * Z;
                end;
            else
                declare
                    P1 : constant := 0.83333_33333_33332_11405e-1;
                    P2 : constant := 0.12500_00000_05037_99174e-1;
                    P3 : constant := 0.22321_42128_59242_58967e-2;
                    P4 : constant := 0.43445_77567_21631_19635e-3;
                begin
                    R := (((P4 * V + P3) * V + P2) * V + P1) * V * Z;
                end;
            end if;
            R   := R + K * R;
            U2  := (R + Z * K) + Z;
            U1  := Mantissa_Type (Integer (M) * 16 - P) * 0.0625;
            Y1  := Reduce (Y);
            Y2  := Mantissa_Type (Y - Floating (Y1));
            W   := Floating (U2) * Y + Floating (U1 * Y2);
            W1  := Reduce (W);
            W2  := Mantissa_Type (W - Floating (W1));
            W   := Floating (W1 + U1 * Y1);
            W1  := Reduce (W);
            W2  := W2 + Mantissa_Type ((W - Floating (W1)));
            W3  := Reduce (Floating (W2));
            Iw1 := Integer (Truncate (16.0 * Floating (W1 + W3)));
            W2  := W2 - W3;
            if W > Floating (Ibigx) then
                Result := Xmax;
                New_Line (Error_Log);
                Put (Error_Log, "X**Y CALLED  X =");
                Put (Error_Log, X);
                Put (Error_Log, "   Y =");
                Put (Error_Log, Y);
                Put (Error_Log, "   TOO LARGE  RETURNED ");
                Put (Error_Log, Result);
                New_Line (Error_Log);
            elsif W < Floating (Ismallx) then
                Result := Zero;
                New_Line (Error_Log);
                Put (Error_Log, "X**Y CALLED  X =");
                Put (Error_Log, X);
                Put (Error_Log, "   Y =");
                Put (Error_Log, Y);
                Put (Error_Log, "   TOO SMALL  RETURNED ");
                Put (Error_Log, Result);
                New_Line (Error_Log);
            else
                if W2 > Mantissa_Type (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;
                if Floating'Base'Digits < 13 then
                    declare
                        Q1 : constant := 0.69314_71805_56341;
                        Q2 : constant := 0.24022_65061_44710;
                        Q3 : constant := 0.55504_04881_30765E-1;
                        Q4 : constant := 0.96162_06595_83789E-2;
                        Q5 : constant := 0.13052_55159_42810E-2;
                    begin
                        Z := ((((Q5 * W2 + Q4) * W2 + Q3) * W2 + Q2)
                             * W2 + Q1) * W2;
                    end;
                else
                    declare
                        Q1 : constant := 0.69314_71805_59945_29629e+0;
                        Q2 : constant := 0.24022_65069_59095_37056e+0;
                        Q3 : constant := 0.55504_10866_40855_95326e-1;
                        Q4 : constant := 0.96181_29059_51724_16964e-2;
                        Q5 : constant := 0.13333_54131_35857_84703e-2;
                        Q6 : constant := 0.15400_29044_09897_64601e-3;
                        Q7 : constant := 0.14928_85268_05956_08186e-4;
                    begin
                        Z := ((((((Q7 * W2 + Q6) * W2 + Q5) * W2 + Q4)
                             * W2 + Q3) * W2 + Q2) * W2 + Q1) * W2;
                    end;
                end if;
                if Floating'Base'Digits < 9 then
                    declare
                        A1 : array (1..17) of Mantissa_Type :=
                            (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#);
                    begin
                        Z := A1 (P + 1) + (A1 (P + 1) * Z);
                    end;
                else
                    declare
                        A1 : array (1..17) of Mantissa_Type :=
                            (8#1.00000_00000_00000_000#,
                             8#0.75222_57505_22220_662#,
                             8#0.72540_30671_75644_416#,
                             8#0.70146_33673_02522_470#,
                             8#0.65642_37462_55323_532#,
                             8#0.63422_21405_21760_440#,
                             8#0.61263_45204_25240_666#,
                             8#0.57204_24347_65401_414#,
                             8#0.55202_36314_77473_630#,
                             8#0.53254_07672_44124_122#,
                             8#0.51377_32652_33052_116#,
                             8#0.47572_46230_11064_104#,
                             8#0.46033_76024_30667_052#,
                             8#0.44341_72334_72542_160#,
                             8#0.42712_70170_76521_364#,
                             8#0.41325_30331_74611_036#,
                             8#0.40000_00000_00000_000#);
                    begin
                        Z := A1 (P + 1) + (A1 (P + 1) * Z);
                    end;
                end if;
                Refloat (M, Z, Result);
            end if;
        end if;
        return Result;
    end "**";
    function Sin (X : Floating) return Floating is
        Sgn, Y       : Floating;
        N            : Integer;
        Xn           : Floating;
        F, G, X1, X2 : Floating;
        Result       : Floating;
        C1           : constant Floating := 8#3.1104#;
        C2           : constant Floating := - 0.89089_10206_76153_73566_17e-5;
        function R (G : Floating) return Floating is
        begin
            if Floating'Base'Digits < 10 then
                declare
                    R1 : constant := - 0.16666_66660_883e-0;
                    R2 : constant := 0.83333_30720_556E-2;
                    R3 : constant := - 0.19840_83282_313E-3;
                    R4 : constant := 0.27523_97106_775E-5;
                    R5 : constant := - 0.23868_34640_601E-7;
                begin
                    return ((((R5 * G + R4) * G + R3) * G + R2) * G + R1) * G;
                end;
            else
                declare
                    R1 : constant := - 0.16666_66666_66666_65052e+0;
                    R2 : constant := 0.83333_33333_33316_50314e-2;
                    R3 : constant := - 0.19841_26984_12018_40457e-3;
                    R4 : constant := 0.27557_31921_01527_56119e-5;
                    R5 : constant := - 0.25052_10679_82745_84544e-7;
                    R6 : constant := 0.16058_93649_03715_89114e-9;
                    R7 : constant := - 0.76429_17806_89104_67734e-12;
                    R8 : constant := 0.27204_79095_78888_46175e-14;
                begin
                    return (((((((R8 * G + R7) * G + R6) * G + R5) * G + R4) *
                        G + R3) * G + R2) * G + R1) * G;
                end;
            end if;
        end R;
    begin
        if X < Zero then
            Sgn := - One;
            Y   := - X;
        else
            Sgn := One;
            Y   := X;
        end if;
        if Y > Ymax then
            New_Line (Error_Log);
            Put (Error_Log,
                 " SIN CALLED WITH ARGUMENT TOO LARGE FOR ACCURACY ");
            Put (Error_Log, X);
            New_Line (Error_Log);
            begin
                N := Integer (Y);
            exception
                when others =>
                    Put_Line (Error_Log, " MUCH TOO LARGE!  RETURNED ZERO");
                    return Zero;
            end;
        end if;
        N  := Integer (Y * One_Over_Pi);
        Xn := Floating (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
            Result := F + F * R (F * F);
        end if;
        return (Sgn * Result);
    end Sin;
    function Cos (X : Floating) return Floating is
        Sgn, Y       : Floating;
        N            : Integer;
        Xn           : Floating;
        F, G, X1, X2 : Floating;
        Result       : Floating;
        C1           : constant Floating := 8#3.1104#;
        C2           : constant Floating := - 0.89089_10206_76153_73566_17e-5;
        function R (G : Floating) return Floating is
        begin
            if Floating'Base'Digits < 10 then
                declare
                    R1 : constant := - 0.16666_66660_883e-0;
                    R2 : constant := 0.83333_30720_556E-2;
                    R3 : constant := - 0.19840_83282_313E-3;
                    R4 : constant := 0.27523_97106_775E-5;
                    R5 : constant := - 0.23868_34640_601E-7;
                begin
                    return ((((R5 * G + R4) * G + R3) * G + R2) * G + R1) * G;
                end;
            else
                declare
                    R1 : constant := - 0.16666_66666_66666_65052e+0;
                    R2 : constant := 0.83333_33333_33316_50314e-2;
                    R3 : constant := - 0.19841_26984_12018_40457e-3;
                    R4 : constant := 0.27557_31921_01527_56119e-5;
                    R5 : constant := - 0.25052_10679_82745_84544e-7;
                    R6 : constant := 0.16058_93649_03715_89114e-9;
                    R7 : constant := - 0.76429_17806_89104_67734e-12;
                    R8 : constant := 0.27204_79095_78888_46175e-14;
                begin
                    return (((((((R8 * G + R7) * G + R6) * G + R5) * G + R4) *
                        G +
                        R3) * G + R2) * G + R1) * G;
                end;
            end if;
        end R;
    begin
        Sgn := 1.0;
        Y   := abs (X) + Pi_Over_Two;
        if Y > Ymax then
            New_Line (Error_Log);
            Put (Error_Log,
                 " COS CALLED WITH ARGUMENT TOO LARGE FOR ACCURACY ");
            Put (Error_Log, X);
            New_Line (Error_Log);
            begin
                N := Integer (Y);
            exception
                when others =>
                    Put_Line (" MUCH TOO LARGE!  RETURNED ONE");
                    return One;
            end;
        end if;
        N  := Integer (Y * One_Over_Pi);
        Xn := Floating (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 : Floating) return Floating is
        Sgn, Y       : Floating;
        N            : Integer;
        Xn           : Floating;
        F, G, X1, X2 : Floating;
        Result       : Floating;
        C1           : constant Floating := 8#1.4442#;
        C2           : constant Floating := - 0.44544_55103_38076_86783_08e-5;
        function R (G : Floating) return Floating is
        begin
            if Floating'Base'Digits < 11 then
                declare
                    P0 : constant Floating := 1.0;
                    P1 : constant Floating := - 0.11136_14403_566;
                    P2 : constant Floating := 0.10751_54738_488E-2;
                    Q0 : constant Floating := 1.0;
                    Q1 : constant Floating := - 0.44469_47720_281;
                    Q2 : constant Floating := 0.15973_39213_300E-1;
                begin
                    return ((P2 * G + P1) * G * F + F) /
                        (((Q2 * G + Q1) * G + 0.5) + 0.5);
                end;
            else
                declare
                    P0 : constant := 1.0;
                    P1 : constant := - 0.13338_35000_64219_60681e+0;
                    P2 : constant := 0.34248_87823_58905_89960e-2;
                    P3 : constant := - 0.17861_70734_22544_26711e-4;
                    Q0 : constant := 1.0;
                    Q1 : constant := - 0.46671_68333_97552_94240e+0;
                    Q2 : constant := 0.25663_83228_94401_12864e-1;
                    Q3 : constant := - 0.31181_53190_70100_27307e-3;
                    Q4 : constant := 0.49819_43399_37865_12270e-6;
                begin
                    return (((P3 * G + P2) * G + P1) * G * F + F) /
                        (((((Q4 * G + Q3) * G + Q2) * G + Q1) * G + 0.5) + 0.5);
                end;
            end if;
        end R;
    begin
        Y := abs (X);
        if Y > Ymax / 2.0 then
            New_Line (Error_Log);
            Put (Error_Log,
                 " TAN CALLED WITH ARGUMENT TOO LARGE FOR ACCURACY ");
            Put (Error_Log, X);
            New_Line (Error_Log);
            begin
                N := Integer (Y);
            exception
                when others =>
                    Put_Line (" MUCH TOO LARGE!  RETURNED ZERO");
                    return Zero;
            end;
        end if;
        N  := Integer (X * Two_Over_Pi);
        Xn := Floating (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 : Floating) return Floating is
        Sgn, Y       : Floating;
        N            : Integer;
        Xn           : Floating;
        F, G, X1, X2 : Floating;
        Result       : Floating;
        Epsilon1     : constant Floating := 1.0 / Xmax;
        C1           : constant Floating := 8#1.4442#;
        C2           : constant Floating := - 0.44544_55103_38076_86783_08e-5;
        function R (G : Floating) return Floating is
        begin
            if Floating'Base'Digits < 11 then
                declare
                    P0 : constant Floating := 1.0;
                    P1 : constant Floating := - 0.11136_14403_566;
                    P2 : constant Floating := 0.10751_54738_488E-2;
                    Q0 : constant Floating := 1.0;
                    Q1 : constant Floating := - 0.44469_47720_281;
                    Q2 : constant Floating := 0.15973_39213_300E-1;
                begin
                    return ((P2 * G + P1) * G * F + F) /
                        (((Q2 * G + Q1) * G + 0.5) + 0.5);
                end;
            else
                declare
                    P0 : constant := 1.0;
                    P1 : constant := - 0.13338_35000_64219_60681e+0;
                    P2 : constant := 0.34248_87823_58905_89960e-2;
                    P3 : constant := - 0.17861_70734_22544_26711e-4;
                    Q0 : constant := 1.0;
                    Q1 : constant := - 0.46671_68333_97552_94240e+0;
                    Q2 : constant := 0.25663_83228_94401_12864e-1;
                    Q3 : constant := - 0.31181_53190_70100_27307e-3;
                    Q4 : constant := 0.49819_43399_37865_12270e-6;
                begin
                    return (((P3 * G + P2) * G + P1) * G * F + F) /
                        (((((Q4 * G + Q3) * G + Q2) * G + Q1) * G + 0.5) + 0.5);
                end;
            end if;
        end R;
    begin
        Y := abs (X);
        if Y < Epsilon1 then
            New_Line (Error_Log);
            Put (Error_Log, " COT CALLED WITH ARGUMENT TOO NEAR ZERO ");
            Put (Error_Log, X);
            New_Line (Error_Log);
            if X < 0.0 then
                return - Xmax;
            else
                return Xmax;
            end if;
        end if;
        if Y > Ymax / 2.0 then
            New_Line (Error_Log);
            Put (Error_Log,
                 " COT CALLED WITH ARGUMENT TOO LARGE FOR ACCURACY ");
            Put (Error_Log, X);
            New_Line (Error_Log);
        end if;
        N  := Integer (X * Two_Over_Pi);
        Xn := Floating (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 : Floating) return Floating is
        G, Y   : Floating;
        Result : Floating;
        function R (G : Floating) return Floating is
        begin
            if Floating'Base'Digits < 11 then
                declare
                    P1 : constant := - 0.27516_55529_0596E1;
                    P2 : constant := 0.29058_76237_4859E1;
                    P3 : constant := - 0.59450_14419_3246;
                    Q0 : constant := - 0.16509_93320_2424E2;
                    Q1 : constant := 0.24864_72896_9164E2;
                    Q2 : constant := - 0.10333_86707_2113E2;
                    Q3 : constant := 1.0;
                begin
                    return (((P3 * G + P2) * G + P1) * G) /
                        (((G + Q2) * G + Q1) * G + Q0);
                end;
            else
                declare
                    P1 : constant := - 0.27368_49452_41642_55994e+2;
                    P2 : constant := 0.57208_22787_78917_31407e+2;
                    P3 : constant := - 0.39688_86299_75048_77339e+2;
                    P4 : constant := 0.10152_52223_38064_63645e+2;
                    P5 : constant := - 0.69674_57344_73506_46411e+0;
                    Q0 : constant := - 0.16421_09671_44985_60795e+3;
                    Q1 : constant := 0.41714_43024_82604_12556e+3;
                    Q2 : constant := - 0.38186_30336_17501_49284e+3;
                    Q3 : constant := 0.15095_27084_10306_04719e+3;
                    Q4 : constant := - 0.23823_85915_36702_38830e+2;
                    Q5 : constant := 1.0;
                begin
                    return (((((P5 * G + P4) * G + P3) * G + P2) * G + P1) * G) /
                        (((((G + Q4) * G + Q3) * G + Q2) * G + Q1) * G + Q0);
                end;
            end if;
        end R;
    begin
        Y := abs (X);
        if Y > Half then
            if Y > 1.0 then
                New_Line (Error_Log);
                Put (Error_Log, " ASIN CALLED FOR ");
                Put (Error_Log, X);
                Put (Error_Log, " (> 1)  TRUNCATED TO 1, CONTINUED");
                New_Line (Error_Log);
                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 : Floating) return Floating is
        G, Y   : Floating;
        Result : Floating;
        function R (G : Floating) return Floating is
        begin
            if Floating'Base'Digits < 11 then
                declare
                    P1 : constant := - 0.27516_55529_0596E1;
                    P2 : constant := 0.29058_76237_4859E1;
                    P3 : constant := - 0.59450_14419_3246;
                    Q0 : constant := - 0.16509_93320_2424E2;
                    Q1 : constant := 0.24864_72896_9164E2;
                    Q2 : constant := - 0.10333_86707_2113E2;
                    Q3 : constant := 1.0;
                begin
                    return (((P3 * G + P2) * G + P1) * G) /
                        (((G + Q2) * G + Q1) * G + Q0);
                end;
            else
                declare
                    P1 : constant := - 0.27368_49452_41642_55994e+2;
                    P2 : constant := 0.57208_22787_78917_31407e+2;
                    P3 : constant := - 0.39688_86299_75048_77339e+2;
                    P4 : constant := 0.10152_52223_38064_63645e+2;
                    P5 : constant := - 0.69674_57344_73506_46411e+0;
                    Q0 : constant := - 0.16421_09671_44985_60795e+3;
                    Q1 : constant := 0.41714_43024_82604_12556e+3;
                    Q2 : constant := - 0.38186_30336_17501_49284e+3;
                    Q3 : constant := 0.15095_27084_10306_04719e+3;
                    Q4 : constant := - 0.23823_85915_36702_38830e+2;
                    Q5 : constant := 1.0;
                begin
                    return (((((P5 * G + P4) * G + P3) * G + P2) * G + P1) * G) /
                        (((((G + Q4) * G + Q3) * G + Q2) * G + Q1) * G + Q0);
                end;
            end if;
        end R;
    begin
        Y := abs (X);
        if Y > Half then
            if Y > 1.0 then
                New_Line (Error_Log);
                Put (Error_Log, " ACOS CALLED FOR ");
                Put (Error_Log, X);
                Put (Error_Log, " (> 1)  TRUNCATED TO 1, CONTINUED");
                New_Line (Error_Log);
                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 : Floating) return Floating is
        F, G : Floating;
        type Region is range 0..3;
        N                : Region;
        Result           : Floating;
        Pi_Over_Six      : constant Floating := Pi_Over_Three / (One + One);
        Sqrt_3           : constant Floating := 1.73205_08075_68877_29353;
        Sqrt_3_Minus_1   : constant Floating := 0.73205_08075_68877_29353;
        Two_Minus_Sqrt_3 : constant Floating := 0.26794_91924_31122_70647;
        function R (G : Floating) return Floating is
        begin
            if Floating'Base'Digits < 10 then
                declare
                    P0 : constant Floating := - 0.14400_83448_74E1;
                    P1 : constant Floating := - 0.72002_68488_98;
                    Q0 : constant Floating := 0.43202_50389_19E1;
                    Q1 : constant Floating := 0.47522_25845_99E1;
                    Q2 : constant Floating := 1.0;
                begin
                    return ((P1 * G + P0) * G) / ((G + Q1) * G + Q0);
                end;
            else
                declare
                    P0 : constant := - 0.13688_76889_41919_26929e+2;
                    P1 : constant := - 0.20505_85519_58616_51981e+2;
                    P2 : constant := - 0.84946_24035_13206_83534e+1;
                    P3 : constant := - 0.83758_29936_81500_59274e+0;
                    Q0 : constant := 0.41066_30668_25757_81263e+2;
                    Q1 : constant := 0.86157_34959_71302_42515e+2;
                    Q2 : constant := 0.59578_43614_25973_44465e+2;
                    Q3 : constant := 0.15024_00116_00285_76121e+2;
                    Q4 : constant := 1.0;
                begin
                    return ((((P3 * G + P2) * G + P1) * G + P0) * G) /
                        ((((G + Q3) * G + Q2) * G + Q1) * G + Q0);
                end;
            end if;
        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 : Floating) return Floating is
        X, Result : Floating;
    begin
        if U = 0.0 then
            if V = 0.0 then
                Result := 0.0;
                New_Line (Error_Log);
                Put (Error_Log, " ATAN2 CALLED WITH 0/0   RETURNED ");
                Put (Error_Log, Result);
                New_Line (Error_Log);
            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 : Floating) return Floating is
        G, W, Y, Z       : Floating;
        Result           : Floating;
        Ybar             : constant Floating := Exp_Large;
        Ln_V             : constant Floating := 8#0.542714#;
        V_Over_2_Minus_1 : constant Floating := 0.13830_27787_96019_02638E-4;
        Wmax             : constant Floating := Ybar - Ln_V + 0.69;
        function R (G : Floating) return Floating is
        begin
            if Floating'Base'Digits < 13 then
                declare
                    P0 : constant Floating := 0.10622_28883_7151E4;
                    P1 : constant Floating := 0.31359_75645_6058E2;
                    P2 : constant Floating := 0.34364_14035_8506;
                    Q0 : constant Floating := 0.63733_73302_1822E4;
                    Q1 : constant Floating := - 0.13051_01250_9199E3;
                    Q2 : constant Floating := 1.0;
                begin
                    return (((P2 * G + P1) * G + P0) * G) / ((G + Q1) * G + Q0);
                end;
            else
                declare
                    P0 : constant := - 0.35181_28343_01771_17881e+6;
                    P1 : constant := - 0.11563_52119_68517_68270e+5;
                    P2 : constant := - 0.16375_79820_26307_51372e+3;
                    P3 : constant := - 0.78966_12741_73570_99479e+0;
                    Q0 : constant := - 0.21108_77005_81062_71242e+7;
                    Q1 : constant := 0.36162_72310_94218_36460e+5;
                    Q2 : constant := - 0.27773_52311_96507_01667e+3;
                    Q3 : constant := 1.0;
                begin
                    return ((((P3 * G + P2) * G + P1) * G + P0) * G) /
                        (((G + Q2) * G + Q1) * G + Q0);
                end;
            end if;
        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 (Error_Log);
                    Put (Error_Log,
                         " SINH CALLED WITH TOO LARGE ARGUMENT  ");
                    Put (Error_Log, X);
                    Put (Error_Log, " RETURN BIG");
                    New_Line (Error_Log);
                    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 : Floating) return Floating is
        G, W, Y, Z       : Floating;
        Result           : Floating;
        Ybar             : constant Floating := Exp_Large;
        Ln_V             : constant := 8#0.542714#;
        V_Over_2_Minus_1 : constant := 0.13830_27787_96019_02638E-4;
        Wmax             : constant Floating := Ybar - Ln_V + 0.69;
    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 (Error_Log);
                Put (Error_Log, " COSH CALLED WITH TOO LARGE ARGUMENT  ");
                Put (Error_Log, X);
                Put (Error_Log, " RETURN BIG");
                New_Line (Error_Log);
                W := Wmax;
            end if;
            Z      := Exp (W);
            Result := Z + V_Over_2_Minus_1 * Z;
        end if;
        return Result;
    end Cosh;
    function Tanh (X : Floating) return Floating is
        G, W, Y, Z  : Floating;
        Result      : Floating;
        Xbig        : constant Floating :=
                    (Log (2.0) + Floating (It + 1) * Log (Beta)) / 2.0;
        Ln_3_Over_2 : constant Floating := 0.54930_61443_34054_84570;
        function R (G : Floating) return Floating is
        begin
            if Floating'Base'Digits < 11 then
                declare
                    P0 : constant Floating := - 0.21063_95800_0245E2;
                    P1 : constant Floating := - 0.93363_47565_2401;
                    Q0 : constant Floating := 0.63191_87401_5582E2;
                    Q1 : constant Floating := 0.28077_65347_0471E2;
                    Q2 : constant Floating := 1.0;
                begin
                    return ((P1 * G + P0) * G) / ((G + Q1) * G + Q0);
                end;
            else
                declare
                    P0 : constant := - 0.16134_11902_39962_28053e+4;
                    P1 : constant := - 0.99225_92967_22360_83313e+2;
                    P2 : constant := - 0.96437_49277_72254_69787e+0;
                    Q0 : constant := 0.48402_35707_19886_88686e+4;
                    Q1 : constant := 0.22337_72071_89623_12926e+4;
                    Q2 : constant := 0.11274_47438_05349_49335e+3;
                    Q3 : constant := 1.0;
                begin
                    return (((P2 * G + P1) * G + P0) * G) /
                        (((G + Q2) * G + Q1) * G + Q0);
                end;
            end if;
        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
    Beta      := Floating (Ibeta);
    Epsilon   := Beta ** (- It / 2);
    Ymax      := Floating (Integer (Pi * (One + One) ** (It / 2)));
    Exp_Large := Log (Xmax) * (One - Epsneg);
    Exp_Small := Log (Xmin) * (One - Epsneg);
    Set_Ran_Key;         --  Set to the default anyway
    if What_To_Do_When_There_Is_An_Error = Return_Default_And_Log then
        Create (Error_Log, Out_File, Name (Standard_Output));  --  Set default
    end if;
end Generic_Math_Functions;