DataMuseum.dk

Presents historical artifacts from the history of:

Rational R1000/400 Tapes

This is an automatic "excavation" of a thematic subset of
artifacts from Datamuseum.dk's BitArchive.

See our Wiki for more about Rational R1000/400 Tapes

Excavated with: AutoArchaeologist - Free & Open Source Software.


top - download
Index: ┃ B T

⟦dd3287cd9⟧ TextFile

    Length: 110870 (0x1b116)
    Types: TextFile
    Names: »B«

Derivation

└─⟦85b835f43⟧ Bits:30000549 8mm tape, Rational 1000, Xlib rev 6.00
    └─ ⟦0c20f784e⟧ »DATA« 
        └─⟦1abbe589f⟧ 
            └─⟦306851c02⟧ 
                └─⟦this⟧ 

TextFile

--/ if R1000 then
with Debug_Tools;
--/ elsif TeleGen2 and then Unix then
--// with Float_Text_Io;
--/ end if;

with Text_Io;

with Xlbt_Arithmetic;  
use Xlbt_Arithmetic;

package body Trig is
------------------------------------------------------------------------------
-- Date      - /Name/ Comment
--
-- 01-FEB-89 - /GEB/ Created based upon the book
--           -  Software Manual for the Elementary Functions by
--           -  William J. Cody, Jr. and William Waite
--           -  Prentice-Hall Series in Computational Mathematics, 1980
-- 05-JUL-90 - /GEB/ Add structured comments, port to TeleGen2/TeleSoft.
------------------------------------------------------------------------------

--/ if not (TeleGen2 and then Unix) then
    package Float_Text_Io is new Text_Io.Float_Io (Float);
--/ end if;

    type Float_Array is array (Natural32 range <>) of Float;

    Degrees_To_Radians : constant Float := Pi / 180.0;

    Sqrt_One_Tenth : constant Float :=  
       0.316227766016837933199889354443271853371955513727833;  
    Sqrt_One_Half  : constant Float :=  
       0.7071067811865475244008443621048490392848359375608;

    Ln_Float_Large : Float;  
    Ln_Float_Small : Float;  
    Ln_2           : constant Float :=  
       0.693147180559945309417232121458176568075500133107812;  
    Log10_E        : constant Float :=  
       0.434294481903251827651128918916605082294397005627727;  
    Ln_10          : constant Float :=  
       2.30258509299404568401799145468436420760110148794495;

    Expx_A1 : Float_Array (1 .. 17);  
    Expx_A2 : Float_Array (1 .. 8);

    Tanh_Xbig : Float;  
    Cot_Eps1  : Float;

--\f

    function Truncate (X : Float) return Float is
------------------------------------------------------------------------------
-- Chop off any non-integral portion of X.  We do this by converting it to
-- an "integer" form and then back again without rounding.
------------------------------------------------------------------------------
        Result : Float;  
        Fr     : Float := abs X;  
        Ex     : Integer32;  
        Tmp    : Integer32;  
        Delt   : Integer32;  
    begin

----Anything less than 1.0 in magnitude is 0.0 for our result.

        if Fr < 1.0 then  
            return 0.0;  
        end if;

----Anything < Integer32'Last can be done using Integer32s.

        if Fr < Float (Integer32'Last) then  
            Tmp := Integer32 (Fr);  
            if Float (Tmp) > Fr then  
                Tmp := Tmp - 1;  
            end if;  
            if X >= 0.0 then  
                return Float (Tmp);  
            else  
                return -Float (Tmp);  
            end if;  
        end if;

----Separate the fraction and the exponent in the number.  Each time we take
--  some bit out of the fraction (in order to put them into our result) we
--  will take out out Delt Float'Machine_Radix "digits" worth.  This guarantees
--  that we don't overflow our Integer32 temporary.

        Fraction_Exponent (abs X, Fr, Ex);  
        case Float'Machine_Radix is  
            when 2 =>  
                Delt := (Integer32'Size - 1);  
            when 4 =>  
                Delt := (Integer32'Size - 1) / 2;  
            when 8 =>  
                Delt := (Integer32'Size - 1) / 3;  
            when 16 =>  
                Delt := (Integer32'Size - 1) / 4;  
            when 10 =>  
                Delt := (Integer32'Size - 1) / 4;  
            when others =>  
                Delt := 1;  
        end case;

----Loop extracting "large" chunks of bits from the fraction.  Keep placing
--  them into the result.  We will be extracting only integral bits.

        Result := 0.0;  
        while Ex > Delt loop  
            Ex  := Ex - Delt;  
            Fr  := Fr * Float (Float'Machine_Radix) ** Integer (Delt);  
            Tmp := Integer32 (Fr);  
            if Float (Tmp) > Fr then        -- Reverse any rounding affects.
                Tmp := Tmp - 1;  
            end if;  
            Fr     := Fr - Float (Tmp);  
            Result := Result * Float (Float'Machine_Radix) ** Integer (Delt) +  
                         Float (Tmp);  
        end loop;

----Remove the last few Ex Float'Machine_Radix "digits" from the fraction.
--  This removes the last of the integral portion of the original number.

        Fr  := Fr * Float (Float'Machine_Radix) ** Integer (Ex);  
        Tmp := Integer32 (Fr);  
        if Float (Tmp) > Fr then  
            Tmp := Tmp - 1;                 -- Reverse any rounding affects.
        end if;  
        Result := Result * Float (Float'Machine_Radix) ** Integer (Ex) +  
                     Float (Tmp);

----Return a result with the correct sign.

        if X < 0.0 then  
            return -Result;  
        else  
            return Result;  
        end if;

    end Truncate;

--\f

    function Truncate (X : Float) return Integer32 is
------------------------------------------------------------------------------
-- Convert the Float to an Integer32; truncate the fraction; do not round.
------------------------------------------------------------------------------
        I : Integer32 := Integer32 (X);  
    begin

        if I >= 0 then  
            if Float (I) > X then  
                return I - 1;  
            else  
                return I;  
            end if;  
        else  
            if Float (I) < X then  
                return I + 1;  
            else  
                return I;  
            end if;  
        end if;

    end Truncate;

--\f

    function Round (X : Float) return Float is
------------------------------------------------------------------------------
-- Chop off any non-integral portion of X.  We do this by converting it to
-- an "integer" form and then back again with rounding.
------------------------------------------------------------------------------
        Result : Float;  
        Fr     : Float := abs X;  
        Ex     : Integer32;  
        Tmp    : Integer32;  
        Delt   : Integer32;  
    begin

----Anything < Integer32'Last can be done using Integer32s.

        if Fr < Float (Integer32'Last) then  
            Tmp := Integer32 (Fr);  
            if X >= 0.0 then  
                return Float (Tmp);  
            else  
                return -Float (Tmp);  
            end if;  
        end if;

----Separate the fraction and the exponent in the number.  Each time we take
--  some bit out of the fraction (in order to put them into our result) we
--  will take out out Delt Float'Machine_Radix "digits" worth.  This guarantees
--  that we don't overflow our Integer32 temporary.

        Fraction_Exponent (abs X, Fr, Ex);  
        case Float'Machine_Radix is  
            when 2 =>  
                Delt := (Integer32'Size - 1);  
            when 4 =>  
                Delt := (Integer32'Size - 1) / 2;  
            when 8 =>  
                Delt := (Integer32'Size - 1) / 3;  
            when 16 =>  
                Delt := (Integer32'Size - 1) / 4;  
            when 10 =>  
                Delt := (Integer32'Size - 1) / 4;  
            when others =>  
                Delt := 1;  
        end case;

----Loop extracting "large" chunks of bits from the fraction.  Keep placing
--  them into the result.  We will be extracting only integral bits.

        Result := 0.0;  
        while Ex > Delt loop  
            Ex  := Ex - Delt;  
            Fr  := Fr * Float (Float'Machine_Radix) ** Integer (Delt);  
            Tmp := Integer32 (Fr);  
            if Float (Tmp) > Fr then        -- Reverse any rounding affects.
                Tmp := Tmp - 1;  
            end if;  
            Fr     := Fr - Float (Tmp);  
            Result := Result * Float (Float'Machine_Radix) ** Integer (Delt) +  
                         Float (Tmp);  
        end loop;

----Remove the last few Ex Float'Machine_Radix "digits" from the fraction.
--  This removes the last of the integral portion of the original number.

        Fr     := Fr * Float (Float'Machine_Radix) ** Integer (Ex);  
        Tmp    := Integer32 (Fr);                -- Allow rounding to occur.
        Result := Result * Float (Float'Machine_Radix) ** Integer (Ex) +  
                     Float (Tmp);

----Return a result with the correct sign.

        if X < 0.0 then  
            return -Result;  
        else  
            return Result;  
        end if;

    end Round;

--\f

    function Round (X : Float) return Integer32 is  
    begin  
        return Integer32 (X);  
    end Round;

--\f

    procedure Fraction_Exponent (X        :     Float;  
                                 Fraction : out Float;  
                                 Exponent : out Integer32) is
------------------------------------------------------------------------------
-- Separate the number into a fraction (0.5 <= abs F < 1.0) and an exponent
-- such that X == Fraction * Float'Machine_Radix ** Exponent
------------------------------------------------------------------------------
        F : Float     := abs X;  
        E : Integer32 := 0;  
    begin

----A zero is 0.0 with 0 exponent.

        if X = 0.0 then  
            Fraction := 0.0;  
            Exponent := 0;  
            return;  
        end if;

----Reduce values > 1.0 to values < 1.0.  We scale by Machine_Radix so there
--  is no roundoff error whatsoever.

        if F >= 1.0 then  
            while F >= Float (Float'Machine_Radix) ** 50 loop  
                E := E + 50;  
                F := F / Float (Float'Machine_Radix) ** 50;  
            end loop;  
            while F >= Float (Float'Machine_Radix) ** 10 loop  
                E := E + 10;  
                F := F / Float (Float'Machine_Radix) ** 10;  
            end loop;  
            while F >= 1.0 loop  
                E := E + 1;  
                F := F / Float (Float'Machine_Radix);  
            end loop;

----Increase values < 1.0/Machine_Radix to values >= 1.0/Machine_Radix.  We
--  scale by Machine_Radix so there is no roundoff error whatsoever.

        elsif F < 1.0 / Float (Float'Machine_Radix) then  
            while F < 1.0 / Float (Float'Machine_Radix) ** 50 loop  
                E := E - 50;  
                F := F * Float (Float'Machine_Radix) ** 50;  
            end loop;  
            while F < 1.0 / Float (Float'Machine_Radix) ** 10 loop  
                E := E - 10;  
                F := F * Float (Float'Machine_Radix) ** 10;  
            end loop;  
            while F < 1.0 / Float (Float'Machine_Radix) loop  
                E := E - 1;  
                F := F * Float (Float'Machine_Radix);  
            end loop;  
        end if;

----Return our results; preserve signs.

        if X < 0.0 then  
            Fraction := -F;  
        else  
            Fraction := F;  
        end if;  
        Exponent := E;

    end Fraction_Exponent;

--\f

    function Odd (X : Float) return Boolean is
------------------------------------------------------------------------------
-- Returns TRUE if Aint(X) is an ODD value.  Any value < 1.0 is even
-- because it is less than 1;
-- and any value > Float'Machine_Radix ** Float'Machine_Mantissa is
-- even because the even/odd bit has been lost due to insufficient
-- precision in the Float representation.
------------------------------------------------------------------------------
        Y     : Float     := abs X;  
        Z     : Float;  
        Shift : Integer32 := 0;  
        Int   : Integer32;  
    begin

----If the value is too small to have an Integer32 value or if it has lost the
--  1-bit that determines odd-ness then we are de-facto even.

        if Y < 1.0 or else  
           Y >= Float (Float'Machine_Radix) ** Float'Machine_Mantissa then  
            return False;  
        end if;

----If necessary we will shift the number around until we can Integer32() it.

        if Y >= Float (Integer32'Last) then

----Determine our shift amount.  Guarantee that Y mod Z is < Integer32'Last.

            case Float'Machine_Radix is  
                when 2 =>  
                    Z := 2.0 ** (Integer32'Size - 2);  
                when 4 =>  
                    Z := 4.0 ** ((Integer32'Size - 2) / 2);  
                when 8 =>  
                    Z := 8.0 ** ((Integer32'Size - 2) / 3);  
                when 16 =>  
                    Z := 16.0 ** ((Integer32'Size - 1) / 4);  
                when 10 =>  
                    Z := 10.0 ** 9;  
                when others =>  
                    raise Program_Error;  
            end case;

----Keep shifting down until we are small enough to fit in an Integer32 when
--  we are Integer32()'ed.  Keep a count of the number of shifts we perform.

            while Y >= Float (Integer32'Last) loop  
                Shift := Shift + 1;  
                Y     := Y / Z;  
            end loop;

----Now recover those shifts.  Remove successive groups of bits until we are
--  left with just the bottom Integer32'Size (more or less) bits from the
--  original Float.

            while Shift > 0 loop  
                Int := Integer32 (Y);  
                if Float (Int) > Y then  
                    Int := Int - 1;  
                end if;  
                Y := Y - Float (Int);  
                Y := Y * Z;  
            end loop;  
        end if;

----Check the last few bits and see if they constitute odd'ness.

        Int := Integer32 (Y);  
        if Float (Int) > Y then  
            Int := Int - 1;  
        end if;  
        if Int rem 2 /= 0 then  
            return True;  
        else  
            return False;  
        end if;

    end Odd;

--\f

    function Sqrt (A : Float) return Float is
------------------------------------------------------------------------------
-- Square-Root  - returns the positive root of positive numbers
--              - raises Constraint_Error for negative number
-- Based upon elements from:
-- Computer Approximations by John F. Hart et. al.; Krieger Publishing 1968
-- pg 89-96
-- Software Manual for the Elementary Functions by W. J. Cody, Jr. and
-- W. Waite; Prentice-Hall 1980; pg 17-34
------------------------------------------------------------------------------
        Adj  : Float;  
        F    : Float;  
        N    : Integer32;  
        Y    : Float;  
        Bits : Natural32;  
    begin

----We don't do negative numbers and sqrt(0.0) = 0.0.

        if A <= 0.0 then  
            if A < 0.0 then  
                raise Constraint_Error;  
            end if;  
            return 0.0;  
        end if;

----Break A into a fraction and an exponent.

        Fraction_Exponent (A, F, N);

----Compute an approximation of the answer for sqrt(f).

        case Float'Machine_Radix is  
            when 2 =>  
                Y := 0.41731 + 0.59016 * F;  
            when others =>  
                Y := 0.580661 + F * 0.5 - 0.086462 / (F + 0.175241);  
        end case;

----Now we put back the exponential information.

        if N rem 2 /= 0 then  
            case Float'Machine_Radix is  
                when 2 =>  
                    Y := Y * 8#0.55202_36314_77473_63110#;  
                when 4 =>  
                    Y := Y * 4#0.2#;  
                when 8 =>  
                    Y := Y * 8#0.26501_17146_37635_71444#;  
                when 10 =>  
                    Y := Y * 0.31622_77660_16837_93320;  
                when 16 =>  
                    Y := Y * 16#0.4#;  
                when others =>  
                    raise Program_Error;  
            end case;  
            N := N + 1;  
        end if;  
        N := N / 2;  
        Y := Y * Float (Float'Machine_Radix) ** Integer (N);

----Now use Newton's approximation to compute the result of sqrt(A).

        Bits := 8;                          -- Currently valid bits (minimum)
        while Bits < Float'Size loop  
            Y    := 0.5 * (Y + A / Y);  
            Bits := Bits + Bits;  
        end loop;

----We now have an answer that is correct to within the bottom bit (or two).
--  See if we can improve on that.  The Newton's approximation cannot
--  do any better because the hardware does (or does not do) rounding and this
--  introduces a "granularity" to the continuous Newton approximation.
--  So, we will sit and diddle with the bottom bit of the answer.  We want a
--  value that when squared is strictly less than our original argument.

        F := Y * Y;  
        if F /= A then  
            Adj := Float (Float'Machine_Radix) **  
                      Integer (N - Float'Machine_Mantissa);  
            if not (Float'Machine_Rounds or else Y = Y + Adj / 2.0) or else  
               not (not Float'Machine_Rounds or else  
                    Y = Y + Adj / 2.00001) or else  
               not (Y /= Y + Adj) or else not (Y = (Y + Adj) - Adj) then  
                raise Program_Error;  
            end if;
            ----Adj is so small that it only diddles the bottom bit of the
            --  mantissa; we so assert.
            begin  
                while F < A loop  
                    Y := Y + Adj;  
                    F := Y * Y;  
                end loop;  
            exception  
                when Numeric_Error =>
                    ----We got an agument that was at/close-to Float'Large
                    --  Reduce Y back by Adj and return it.
                    Y := Y - Adj;  
                    return Y;  
            end;  
            while F > A loop  
                Y := Y - Adj;  
                F := Y * Y;  
            end loop;  
        end if;

----Return the correct result.

        return Y;

    end Sqrt;

--\f

    function Cbrt (A : Float) return Float is
------------------------------------------------------------------------------
-- Cube-Root    - returns the real root
--
-- Based upon elements from:
-- Computer Approximations by John F. Hart et. al.; Krieger Publishing 1968
-- pg 89-96
-- Software Manual for the Elementary Functions by W. J. Cody, Jr. and
-- W. Waite; Prentice-Hall 1980; pg 17-34
------------------------------------------------------------------------------
        Aa   : Float := abs A;  
        Adj  : Float;  
        F    : Float;  
        N    : Integer32;  
        Y    : Float;  
        K    : Integer32;  
        Bits : Natural32;  
    begin

----cbrt(0.0) = 0.0.

        if A = 0.0 then  
            return 0.0;  
        end if;

----Break A into a fraction and an exponent.

        Fraction_Exponent (Aa, F, N);

----Compute an approximation of the answer.

        Y := (((-0.82075_488E-01 * F) + 0.42038_11298) * F + 0.78999_824) * F +  
                0.43261_46E-01;  
        Y := Y / (F + 0.17178_3523);  
        K := N rem 3;  
        if K > 0 then
            ----1/cbrt(3) = .69336.....
            K := 3 - K;  
            Y := Y * 7.937005259840997373758528196362E-001 ** Integer (K);  
            N := N + K;  
        elsif K < 0 then  
            Y := Y * 7.937005259840997373758528196362E-001 ** Integer (abs K);  
            N := N - K;  
        end if;  
        N := N / 3;  
        Y := Y * Float (Float'Machine_Radix) ** Integer (N);

----Each time we do the equation we double our correct digits.  We need
--  Float_Digits worth correct.  Loop until this is true.

        Bits := 8;  
        while Bits < Float'Size loop  
            Y    := (Y + Y + Aa / (Y * Y)) / 3.0;  
            Bits := Bits + Bits;  
        end loop;

----We now have an answer that is correct to within the bottom bit (or two).
--  See if we can improve on that.  The Newton's approximation cannot
--  do any better because the hardware does (or does not do) rounding and this
--  introduces a "granularity" to the continuous Newton approximation.
--  So, we will sit and diddle with the bottom bit of the answer.  We want a
--  value that when squared is strictly less than our original argument.

        F := Y * Y * Y;  
        if F /= Aa then  
            Adj := Float (Float'Machine_Radix) **  
                      Integer (N - Float'Machine_Mantissa);  
            if not (Float'Machine_Rounds or else Y = Y + Adj / 2.0) or else  
               not (not Float'Machine_Rounds or else  
                    Y = Y + Adj / 2.00001) or else  
               not (Y /= Y + Adj) or else  
               not (Y = (Y + Adj) - Adj) then  
                raise Program_Error;  
            end if;
            ----Adj is so small that it only diddles the bottom bit of the
            --  mantissa; we so assert.
            begin  
                while F < Aa loop  
                    Y := Y + Adj;  
                    F := Y * Y * Y;  
                end loop;  
            exception  
                when Numeric_Error =>
                    ----We got an agument that was at/close-to Float'Large
                    --  Reduce Y back by Adj and return it.
                    Y := Y - Adj;  
                    return Y;  
            end;  
            while F > Aa loop  
                Y := Y - Adj;  
                F := Y * Y * Y;  
            end loop;  
        end if;

----Return the correct result.

        if A < 0.0 then  
            return -Y;  
        else  
            return Y;  
        end if;

    end Cbrt;

--\f

    function Sin (A : Float) return Float is  
        R    : Float;  
        F    : Float;  
        G    : Float;  
        Y    : Float;  
        Xn   : Float;  
        Tmp  : Float;  
        Sign : Boolean;  
    begin

----Get Y := abs A where Rad is also in radians.

        if In_Radians then  
            Y := A;  
        else  
            Y := A * Degrees_To_Radians;  
        end if;  
        if Y < 0.0 then  
            Y    := -Y;  
            Sign := True;  
        else  
            Sign := False;  
        end if;

----Check for ridiculously large values of Y.

        if Y >= Pi * Float (Float'Machine_Radix) **  
                        (Float'Machine_Mantissa / 2) then  
            return 0.0;  
        end if;

----Form XN and check for even/odd-ness.

        Xn := Round (Y * 0.31830_98861_83790_67154);  
        if Odd (Xn) then  
            Sign := not Sign;  
        end if;

----Compute the f portion of the argument.

        Tmp := Truncate (Y);  
        if (Float'Machine_Radix = 2 and then  
            Float'Machine_Mantissa >= 33) or else  
           (Float'Machine_Radix = 4 and then  
            Float'Machine_Mantissa * 2 >= 33) or else  
           (Float'Machine_Radix = 8 and then  
            Float'Machine_Mantissa * 3 >= 33) or else  
           (Float'Machine_Radix = 16 and then  
            Float'Machine_Mantissa * 4 >= 33) then  
            F := ((Tmp - Xn * 8#3.1104#) + (Y - Tmp)) +  
                    Xn * 8.9089_10206_76153_73566_17E-06;  
        else  
            F := ((Tmp - Xn * 8#3.11#) + (Y - Tmp)) - Xn * 9.6765_35897_93E-04;  
        end if;

----See if there is no point to computing sin(x).

        if abs F < Float (Float'Machine_Radix) **  
                      (-Float'Machine_Mantissa / 2) then  
            if Sign then  
                return -F;  
            else  
                return F;  
            end if;  
        end if;

----Compute an answer.

        G := F * F;  
        if (Float'Machine_Radix = 2 and then  
            Float'Machine_Mantissa >= 51) or else  
           (Float'Machine_Radix = 4 and then  
            Float'Machine_Mantissa * 2 >= 51) or else  
           (Float'Machine_Radix = 8 and then  
            Float'Machine_Mantissa * 3 >= 51) or else  
           (Float'Machine_Radix = 16 and then  
            Float'Machine_Mantissa * 4 >= 51) or else  
           (Float'Machine_Radix = 10 and then Float'Machine_Mantissa >= 16) then  
            R := (((((((0.27204_79095_78888_46175E-014 * G -  
                        0.76429_17806_89104_67734E-012) * G +  
                       0.16058_93649_03715_89114E-09) * G -  
                      0.25052_10679_82745_84544E-07) * G +  
                     0.27557_31921_01527_56119E-05) * G -  
                    0.19841_26984_12018_40457E-03) * G +  
                   0.83333_33333_33316_50314E-02) * G -  
                  0.16666_66666_66666_65052E+00) * G;  
        elsif (Float'Machine_Radix = 2 and then  
               Float'Machine_Mantissa >= 33) or else  
              (Float'Machine_Radix = 4 and then  
               Float'Machine_Mantissa * 2 >= 33) or else  
              (Float'Machine_Radix = 8 and then  
               Float'Machine_Mantissa * 3 >= 33) or else  
              (Float'Machine_Radix = 16 and then  
               Float'Machine_Mantissa * 4 >= 33) or else  
              (Float'Machine_Radix = 10 and then  
               Float'Machine_Mantissa >= 11) then  
            R := ((((((-0.73706_62775_07114_174E-012 * G +  
                       0.16047_84463_23816_900E-09) * G -  
                      0.25051_87088_34705_760E-07) * G +  
                     0.27557_31642_12926_457E-05) * G -  
                    0.19841_26982_32225_068E-03) * G +  
                   0.8333_33333_27592_139E-02) * G -  
                  0.16666_66666_66659_653E+00) * G;  
        elsif (Float'Machine_Radix = 2 and then  
               Float'Machine_Mantissa >= 25) or else  
              (Float'Machine_Radix = 4 and then  
               Float'Machine_Mantissa * 2 >= 25) or else  
              (Float'Machine_Radix = 8 and then  
               Float'Machine_Mantissa * 3 >= 25) or else  
              (Float'Machine_Radix = 16 and then  
               Float'Machine_Mantissa * 4 >= 25) or else  
              (Float'Machine_Radix = 10 and then  
               Float'Machine_Mantissa >= 9) then  
            R := ((((-0.23868_34640_601E-07 * G +  
                     0.27523_97106_775E-05) * G -  
                    0.19840_83282_313E-03) * G +  
                   0.83333_30720_556E-02) * G -  
                  0.16666_66660_883E+00) * G;  
        else  
            R := (((0.26019_03036E-05 * G -  
                    0.19807_41872E-03) * G +  
                   0.83330_25139E-02) * G -  
                  0.16666_65668E+00) * G;  
        end if;

----Now finish off the result.

        R := F + F * R;  
        if Sign then  
            return -R;  
        else  
            return R;  
        end if;

    end Sin;

--\f

    function Cos (A : Float) return Float is  
        R    : Float;  
        F    : Float;  
        G    : Float;  
        Y    : Float;  
        Xn   : Float;  
        Tmp  : Float;  
        Sign : Boolean;  
    begin

----Get Y := abs A where Rad is also in radians.

        if In_Radians then  
            Y := abs A;  
        else  
            Y := abs A * Degrees_To_Radians;  
        end if;  
        Y    := Y + 1.57079_63267_94896_61923;  
        Sign := False;

----Check for ridiculously large values of Y.

        if Y >= Pi * Float (Float'Machine_Radix) **  
                        (Float'Machine_Mantissa / 2) then  
            return 0.0;  
        end if;

----Form XN and check for even/odd-ness.

        Xn := Round (Y * 0.31830_98861_83790_67154);  
        if Odd (Xn) then  
            Sign := not Sign;  
        end if;  
        Xn := Xn - 0.5;

----Compute the f portion of the argument.

        Tmp := Truncate (abs A);  
        if (Float'Machine_Radix = 2 and then  
            Float'Machine_Mantissa >= 33) or else  
           (Float'Machine_Radix = 4 and then  
            Float'Machine_Mantissa * 2 >= 33) or else  
           (Float'Machine_Radix = 8 and then  
            Float'Machine_Mantissa * 3 >= 33) or else  
           (Float'Machine_Radix = 16 and then  
            Float'Machine_Mantissa * 4 >= 33) then  
            F := ((Tmp - Xn * 8#3.1104#) + (abs A - Tmp)) +  
                    Xn * 8.9089_10206_76153_73566_17E-06;  
        else  
            F := ((Tmp - Xn * 8#3.11#) + (abs A - Tmp)) -  
                    Xn * 9.6765_35897_93E-04;  
        end if;

----See if there is no point to computing cos(x).

        if abs F < Float (Float'Machine_Radix) **  
                      (-Float'Machine_Mantissa / 2) then  
            if Sign then  
                return -F;  
            else  
                return F;  
            end if;  
        end if;

----Compute an answer.

        G := F * F;  
        if (Float'Machine_Radix = 2 and then  
            Float'Machine_Mantissa >= 51) or else  
           (Float'Machine_Radix = 4 and then  
            Float'Machine_Mantissa * 2 >= 51) or else  
           (Float'Machine_Radix = 8 and then  
            Float'Machine_Mantissa * 3 >= 51) or else  
           (Float'Machine_Radix = 16 and then  
            Float'Machine_Mantissa * 4 >= 51) or else  
           (Float'Machine_Radix = 10 and then Float'Machine_Mantissa >= 16) then  
            R := (((((((0.27204_79095_78888_46175E-014 * G -  
                        0.76429_17806_89104_67734E-012) * G +  
                       0.16058_93649_03715_89114E-09) * G -  
                      0.25052_10679_82745_84544E-07) * G +  
                     0.27557_31921_01527_56119E-05) * G -  
                    0.19841_26984_12018_40457E-03) * G +  
                   0.83333_33333_33316_50314E-02) * G -  
                  0.16666_66666_66666_65052E+00) * G;  
        elsif (Float'Machine_Radix = 2 and then  
               Float'Machine_Mantissa >= 33) or else  
              (Float'Machine_Radix = 4 and then  
               Float'Machine_Mantissa * 2 >= 33) or else  
              (Float'Machine_Radix = 8 and then  
               Float'Machine_Mantissa * 3 >= 33) or else  
              (Float'Machine_Radix = 16 and then  
               Float'Machine_Mantissa * 4 >= 33) or else  
              (Float'Machine_Radix = 10 and then  
               Float'Machine_Mantissa >= 11) then  
            R := ((((((-0.73706_62775_07114_174E-012 * G +  
                       0.16047_84463_23816_900E-09) * G -  
                      0.25051_87088_34705_760E-07) * G +  
                     0.27557_31642_12926_457E-05) * G -  
                    0.19841_26982_32225_068E-03) * G +  
                   0.8333_33333_27592_139E-02) * G -  
                  0.16666_66666_66659_653E+00) * G;  
        elsif (Float'Machine_Radix = 2 and then  
               Float'Machine_Mantissa >= 25) or else  
              (Float'Machine_Radix = 4 and then  
               Float'Machine_Mantissa * 2 >= 25) or else  
              (Float'Machine_Radix = 8 and then  
               Float'Machine_Mantissa * 3 >= 25) or else  
              (Float'Machine_Radix = 16 and then  
               Float'Machine_Mantissa * 4 >= 25) or else  
              (Float'Machine_Radix = 10 and then  
               Float'Machine_Mantissa >= 9) then  
            R := ((((-0.23868_34640_601E-07 * G +  
                     0.27523_97106_775E-05) * G -  
                    0.19840_83282_313E-03) * G +  
                   0.83333_30720_556E-02) * G -  
                  0.16666_66660_883E+00) * G;  
        else  
            R := (((0.26019_03036E-05 * G -  
                    0.19807_41872E-03) * G +  
                   0.83330_25139E-02) * G -  
                  0.16666_65668E+00) * G;  
        end if;

----Now finish off the result.

        R := F + F * R;  
        if Sign then  
            return -R;  
        else  
            return R;  
        end if;

    end Cos;

--\f

    function Tan (A : Float) return Float is  
        Aa   : Float;  
        F    : Float;  
        X1   : Float;  
        Xnum : Float;  
        Xden : Float;  
        Even : Boolean;  
        Xn   : Float;  
    begin

----Make sure we have the correct units.

        if In_Radians then  
            Aa := A;  
        else  
            Aa := A * Degrees_To_Radians;  
        end if;  
        F := abs Aa;

----See if our argument is much too large to make any sense.

        if F > Truncate (Float (Float'Machine_Radix) **  
                         (Float'Machine_Mantissa / 2) * Pi_Over_2) then  
            return 0.0;  
        end if;

----Form XN and check for even/odd-ness.

        Xn   := Round (Aa * 0.63661_97723_67581_34308);  
        Even := not Odd (Xn);

----Compute f.

        X1 := Truncate (Aa);  
        if Float'Machine_Radix /= 10 then  
            F := ((X1 - Xn * 8#1.4442#) + (Aa - X1)) +  
                    Xn * 4.4544_55103_38706_86783_08E-06;  
        else  
            F := ((X1 - Xn * 1.57079) + (Aa - X1)) +  
                    Xn * 6.32679489661923132169163975144E-06;  
        end if;  
        if not (abs F <= Pi_Over_4) then  
            raise Program_Error;  
        end if;

----See if f is very very small.

        if abs F < Float (Float'Machine_Radix) **  
                      (-Float'Machine_Mantissa / 2) then  
            Xnum := F;  
            Xden := 1.0;  
        else  
            X1 := F * F;  
            if (Float'Machine_Radix = 2 and then  
                Float'Machine_Mantissa >= 53) or else  
               (Float'Machine_Radix = 4 and then  
                Float'Machine_Mantissa * 2 >= 53) or else  
               (Float'Machine_Radix = 8 and then  
                Float'Machine_Mantissa * 3 >= 53) or else  
               (Float'Machine_Radix = 16 and then  
                Float'Machine_Mantissa * 4 >= 53) or else  
               (Float'Machine_Radix = 10 and then  
                Float'Machine_Mantissa >= 17) then  
                Xnum := (((-0.17861_70734_22544_26711E-04 * X1 +  
                           0.34248_87823_58905_89960E-02) * X1 -  
                          0.13338_35000_64219_60681E+00) * X1) * F + F;  
                Xden := (((0.49819_43399_37865_12270E-06 * X1 -  
                           0.31181_53190_70100_27307E-03) * X1 +  
                          0.25663_83228_94401_12864E-01) * X1 -  
                         0.46671_68333_97552_94240E+00) * X1 +  
                        0.5;  
            elsif (Float'Machine_Radix = 2 and then  
                   Float'Machine_Mantissa >= 33) or else  
                  (Float'Machine_Radix = 4 and then  
                   Float'Machine_Mantissa * 2 >= 33) or else  
                  (Float'Machine_Radix = 8 and then  
                   Float'Machine_Mantissa * 3 >= 33) or else  
                  (Float'Machine_Radix = 16 and then  
                   Float'Machine_Mantissa * 4 >= 33) or else  
                  (Float'Machine_Radix = 10 and then  
                   Float'Machine_Mantissa >= 11) then  
                Xnum := (((-0.74836_34966_61206_5149E-05 * X1 +  
                           0.28059_18241_16998_8906E-02) * X1 -  
                          0.12828_34704_09574_3847E+00) * X1) * F + F;  
                Xden := ((-0.20844_80442_20387_0948E-03 * X1 +  
                          0.23344_85282_20687_2802E-01) * X1 -  
                         0.46161_68037_42904_8840E+00) * X1 +  
                        0.5;  
            elsif (Float'Machine_Radix = 2 and then  
                   Float'Machine_Mantissa >= 25) or else  
                  (Float'Machine_Radix = 4 and then  
                   Float'Machine_Mantissa * 2 >= 25) or else  
                  (Float'Machine_Radix = 8 and then  
                   Float'Machine_Mantissa * 3 >= 25) or else  
                  (Float'Machine_Radix = 16 and then  
                   Float'Machine_Mantissa * 4 >= 25) or else  
                  (Float'Machine_Radix = 10 and then  
                   Float'Machine_Mantissa >= 8) then  
                Xnum := ((0.10751_54738_488E-02 * X1 -  
                          0.11136_14403_566E+00) * X1) * F + F;  
                Xden := (0.15973_39213_300E-01 * X1 -  
                         0.44469_47720_281E+00) * X1 +  
                        0.5;  
            else  
                Xnum := (-0.95801_7723E-01 * X1) * F + F;  
                Xden := (0.97168_5835E-02 * X1 -  
                         0.42913_5777E+00) * X1 +  
                        0.5;  
            end if;  
            Xden := Xden + 0.5;  
        end if;

----If XN was odd then change sign of result.

        if not Even then  
            Xnum := Xden / (-Xnum);  
        else  
            Xnum := Xnum / Xden;  
        end if;  
        return Xnum;

    end Tan;

--\f

    function Sec (A : Float) return Float is  
    begin  
        return 1.0 / Cos (A);  
    end Sec;

--\f

    function Csc (A : Float) return Float is  
    begin  
        return 1.0 / Sin (A);  
    end Csc;

--\f

    function Cot (A : Float) return Float is  
        Aa   : Float;  
        F    : Float;  
        X1   : Float;  
        Xnum : Float;  
        Xden : Float;  
        Even : Boolean;  
        Xn   : Float;  
    begin

----Make sure we have the correct units.

        if In_Radians then  
            Aa := A;  
        else  
            Aa := A * Degrees_To_Radians;  
        end if;  
        F := abs Aa;

----See if our argument is much too small to handle.

        if F < Cot_Eps1 then  
            if Aa = 0.0 then  
                return Float'Last;  
            elsif Aa < 0.0 then  
                return Float'First;  
            else  
                return Float'Last;  
            end if;  
        end if;

----See if our argument is much too large to make any sense.

        if F > Truncate (Float (Float'Machine_Radix) **  
                         (Float'Machine_Mantissa / 2) * Pi_Over_2) then  
            return 0.0;  
        end if;

----Form XN and check for even/odd-ness.

        Xn   := Round (Aa * 0.63661_97723_67581_34308);  
        Even := not Odd (Xn);

----Compute f.

        X1 := Truncate (Aa);  
        if Float'Machine_Radix /= 10 then  
            F := ((X1 - Xn * 8#1.4442#) + (Aa - X1)) +  
                    Xn * 4.4544_55103_38706_86783_08E-06;  
        else  
            F := ((X1 - Xn * 1.57079) + (Aa - X1)) +  
                    Xn * 6.32679489661923132169163975144E-06;  
        end if;

----See if f is very very small.

        if abs F < Float (Float'Machine_Radix) **  
                      (-Float'Machine_Mantissa / 2) then  
            Xnum := F;  
            Xden := 1.0;  
        else  
            X1 := F * F;  
            if (Float'Machine_Radix = 2 and then  
                Float'Machine_Mantissa >= 53) or else  
               (Float'Machine_Radix = 4 and then  
                Float'Machine_Mantissa * 2 >= 53) or else  
               (Float'Machine_Radix = 8 and then  
                Float'Machine_Mantissa * 3 >= 53) or else  
               (Float'Machine_Radix = 16 and then  
                Float'Machine_Mantissa * 4 >= 53) or else  
               (Float'Machine_Radix = 10 and then  
                Float'Machine_Mantissa >= 17) then  
                Xnum := (((-0.17861_70734_22544_26711E-04 * X1 +  
                           0.34248_87823_58905_89960E-02) * X1 -  
                          0.13338_35000_64219_60681E+00) * X1) * F + F;  
                Xden := (((0.49819_43399_37865_12270E-06 * X1 -  
                           0.31181_53190_70100_27307E-03) * X1 +  
                          0.25663_83228_94401_12864E-01) * X1 -  
                         0.46671_68333_97552_94240E+00) * X1 +  
                        0.5;  
            elsif (Float'Machine_Radix = 2 and then  
                   Float'Machine_Mantissa >= 33) or else  
                  (Float'Machine_Radix = 4 and then  
                   Float'Machine_Mantissa * 2 >= 33) or else  
                  (Float'Machine_Radix = 8 and then  
                   Float'Machine_Mantissa * 3 >= 33) or else  
                  (Float'Machine_Radix = 16 and then  
                   Float'Machine_Mantissa * 4 >= 33) or else  
                  (Float'Machine_Radix = 10 and then  
                   Float'Machine_Mantissa >= 11) then  
                Xnum := (((-0.74836_34966_61206_5149E-05 * X1 +  
                           0.28059_18241_16998_8906E-02) * X1 -  
                          0.12828_34704_09574_3847E+00) * X1) * F + F;  
                Xden := ((-0.20844_80442_20387_0948E-03 * X1 +  
                          0.23344_85282_20687_2802E-01) * X1 -  
                         0.46161_68037_42904_8840E+00) * X1 +  
                        0.5;  
            elsif (Float'Machine_Radix = 2 and then  
                   Float'Machine_Mantissa >= 25) or else  
                  (Float'Machine_Radix = 4 and then  
                   Float'Machine_Mantissa * 2 >= 25) or else  
                  (Float'Machine_Radix = 8 and then  
                   Float'Machine_Mantissa * 3 >= 25) or else  
                  (Float'Machine_Radix = 16 and then  
                   Float'Machine_Mantissa * 4 >= 25) or else  
                  (Float'Machine_Radix = 10 and then  
                   Float'Machine_Mantissa >= 8) then  
                Xnum := ((0.10751_54738_488E-02 * X1 -  
                          0.11136_14403_566E+00) * X1) * F + F;  
                Xden := (0.15973_39213_300E-01 * X1 -  
                         0.44469_47720_281E+00) * X1 +  
                        0.5;  
            else  
                Xnum := (-0.95801_7723E-01 * X1) * F + F;  
                Xden := (0.97168_5835E-02 * X1 -  
                         0.42913_5777E+00) * X1 +  
                        0.5;  
            end if;  
            Xden := Xden + 0.5;  
        end if;

----If XN was odd then change sign of result.

        if not Even then  
            Xnum := -Xnum / Xden;  
        else  
            Xnum := Xden / Xnum;  
        end if;  
        return Xnum;

    end Cot;

--\f

    function Asin (A : Float) return Float is  
        Y : Float := abs A;  
        F : Float;  
        G : Float;  
        I : Integer32;  
    begin

----Very small arguments are equal to the answer.

        if Y < Float (Float'Machine_Radix) ** (-Float'Machine_Mantissa / 2) then  
            I := 0;  
            G := Y;  
        else

----Compute the intermediate G value differently depending upon the
--  magnitude of the argument.

            if Y > 0.5 then  
                I := 1;  
                if Y > 1.0 then  
                    raise Constraint_Error;  
                end if;  
                F := ((0.5 - Y) + 0.5) / 2.0;  
                Y := Sqrt (F);  
                Y := -(Y + Y);  
                G := F;  
            else  
                I := 0;  
                G := Y * Y;  
            end if;

----Approximate our result.

            if (Float'Machine_Radix = 2 and then  
                Float'Machine_Mantissa >= 49) or else  
               (Float'Machine_Radix = 4 and then  
                Float'Machine_Mantissa * 2 >= 49) or else  
               (Float'Machine_Radix = 8 and then  
                Float'Machine_Mantissa * 3 >= 49) or else  
               (Float'Machine_Radix = 16 and then  
                Float'Machine_Mantissa * 4 >= 49) or else  
               (Float'Machine_Radix = 10 and then  
                Float'Machine_Mantissa >= 15) then  
                F := ((((-0.69674_57344_73506_46411E+00 * G +  
                         0.10152_52223_38064_63645E+02) * G -  
                        0.39688_86299_75048_77339E+02) * G +  
                       0.57208_22787_78917_31407E+02) * G -  
                      0.27368_49452_41642_55994E+02) * G  
                      / (((((G - 0.23823_85915_36702_38830E+02) * G +  
                           0.15095_27084_10306_04719E+03) * G -  
                          0.38186_30336_17501_49284E+03) * G +  
                         0.41714_43024_82604_12556E+03) * G -  
                        0.16421_09671_44985_60795E+03);  
            elsif (Float'Machine_Radix = 2 and then  
                   Float'Machine_Mantissa >= 37) or else  
                  (Float'Machine_Radix = 4 and then  
                   Float'Machine_Mantissa * 2 >= 37) or else  
                  (Float'Machine_Radix = 8 and then  
                   Float'Machine_Mantissa * 3 >= 37) or else  
                  (Float'Machine_Radix = 16 and then  
                   Float'Machine_Mantissa * 4 >= 37) or else  
                  (Float'Machine_Radix = 10 and then  
                   Float'Machine_Mantissa >= 12) then  
                F := (((-0.65404_06899_93350_09E+00 * G +  
                        0.59683_15761_77515_34E+01) * G -  
                       0.13428_70791_34253_12E+02) * G +  
                      0.85372_16436_67719_50E+01) * G  
                      / ((((G - 0.16429_55755_74951_70E+02) * G +  
                          0.68729_59765_38088_06E+02) * G -  
                         0.10362_27318_64014_80E+03) * G +  
                        0.51223_29862_01096_91E+02);  
            elsif (Float'Machine_Radix = 2 and then  
                   Float'Machine_Mantissa >= 25) or else  
                  (Float'Machine_Radix = 4 and then  
                   Float'Machine_Mantissa * 2 >= 25) or else  
                  (Float'Machine_Radix = 8 and then  
                   Float'Machine_Mantissa * 3 >= 25) or else  
                  (Float'Machine_Radix = 16 and then  
                   Float'Machine_Mantissa * 4 >= 25) or else  
                  (Float'Machine_Radix = 10 and then  
                   Float'Machine_Mantissa >= 8) then  
                F := ((-0.59450_14419_3246E+00 * G +  
                       0.29058_76237_4859E+01) * G -  
                      0.27516_55529_0596E+01) * G  
                      / (((G - 0.10333_86707_2113E+02) * G +  
                         0.24864_72896_9164E+02) * G -  
                        0.16509_93320_2424E+02);  
            else  
                F := (-0.50440_0557E+00 * G +  
                      0.93393_5835E+00) * G  
                      / ((G - 0.55484_6723E+01) * G +  
                        0.56036_3004E+01);  
            end if;  
            G := Y + Y * F;

----Adjust our answer to a different quadrant if necessary.

            if I = 1 then  
                G := 8#0.62207_73250_42055_06043# +  
                        (8#0.62207_73250_42055_06043# + G);  
            end if;  
        end if;

----Return the result with the correct sign.

        if not In_Radians then  
            G := G / Degrees_To_Radians;  
        end if;  
        if A < 0.0 then  
            return -G;  
        else  
            return G;  
        end if;

    end Asin;

--\f

    function Acos (A : Float) return Float is  
        Y : Float := abs A;  
        F : Float;  
        G : Float;  
        I : Integer32;  
    begin

----Very small arguments are equal to the answer.

        if Y < Float (Float'Machine_Radix) ** (-Float'Machine_Mantissa / 2) then  
            I := 1;  
            G := Y;  
        else

----Compute the intermediate G value differently depending upon the
--  magnitude of the argument.

            if Y > 0.5 then  
                I := 0;  
                if Y > 1.0 then  
                    raise Constraint_Error;  
                end if;  
                F := ((0.5 - Y) + 0.5) / 2.0;  
                Y := Sqrt (F);  
                Y := -(Y + Y);  
                G := F;  
            else  
                I := 1;  
                G := Y * Y;  
            end if;

----Approximate our result.

            if (Float'Machine_Radix = 2 and then  
                Float'Machine_Mantissa >= 49) or else  
               (Float'Machine_Radix = 4 and then  
                Float'Machine_Mantissa * 2 >= 49) or else  
               (Float'Machine_Radix = 8 and then  
                Float'Machine_Mantissa * 3 >= 49) or else  
               (Float'Machine_Radix = 16 and then  
                Float'Machine_Mantissa * 4 >= 49) or else  
               (Float'Machine_Radix = 10 and then  
                Float'Machine_Mantissa >= 15) then  
                F := ((((-0.69674_57344_73506_46411E+00 * G +  
                         0.10152_52223_38064_63645E+02) * G -  
                        0.39688_86299_75048_77339E+02) * G +  
                       0.57208_22787_78917_31407E+02) * G -  
                      0.27368_49452_41642_55994E+02) * G  
                      / (((((G - 0.23823_85915_36702_38830E+02) * G +  
                           0.15095_27084_10306_04719E+03) * G -  
                          0.38186_30336_17501_49284E+03) * G +  
                         0.41714_43024_82604_12556E+03) * G -  
                        0.16421_09671_44985_60795E+03);  
            elsif (Float'Machine_Radix = 2 and then  
                   Float'Machine_Mantissa >= 37) or else  
                  (Float'Machine_Radix = 4 and then  
                   Float'Machine_Mantissa * 2 >= 37) or else  
                  (Float'Machine_Radix = 8 and then  
                   Float'Machine_Mantissa * 3 >= 37) or else  
                  (Float'Machine_Radix = 16 and then  
                   Float'Machine_Mantissa * 4 >= 37) or else  
                  (Float'Machine_Radix = 10 and then  
                   Float'Machine_Mantissa >= 12) then  
                F := (((-0.65404_06899_93350_09E+00 * G +  
                        0.59683_15761_77515_34E+01) * G -  
                       0.13428_70791_34253_12E+02) * G +  
                      0.85372_16436_67719_50E+01) * G  
                      / ((((G - 0.16429_55755_74951_70E+02) * G +  
                          0.68729_59765_38088_06E+02) * G -  
                         0.10362_27318_64014_80E+03) * G +  
                        0.51223_29862_01096_91E+02);  
            elsif (Float'Machine_Radix = 2 and then  
                   Float'Machine_Mantissa >= 25) or else  
                  (Float'Machine_Radix = 4 and then  
                   Float'Machine_Mantissa * 2 >= 25) or else  
                  (Float'Machine_Radix = 8 and then  
                   Float'Machine_Mantissa * 3 >= 25) or else  
                  (Float'Machine_Radix = 16 and then  
                   Float'Machine_Mantissa * 4 >= 25) or else  
                  (Float'Machine_Radix = 10 and then  
                   Float'Machine_Mantissa >= 8) then  
                F := ((-0.59450_14419_3246E+00 * G +  
                       0.29058_76237_4859E+01) * G -  
                      0.27516_55529_0596E+01) * G  
                      / (((G - 0.10333_86707_2113E+02) * G +  
                         0.24864_72896_9164E+02) * G -  
                        0.16509_93320_2424E+02);  
            else  
                F := (-0.50440_0557E+00 * G +  
                      0.93393_5835E+00) * G  
                      / ((G - 0.55484_6723E+01) * G +  
                        0.56036_3004E+01);  
            end if;  
            G := Y + Y * F;  
        end if;

----Adjust our answer to a different quadrant if necessary.

        if A < 0.0 then  
            if I = 1 then  
                G := 8#0.62207_73250_42055_06043# +  
                        (8#0.62207_73250_42055_06043# + G);  
            else  
                G := 8#1.44417_66521_04132_14107# +  
                        (8#1.44417_66521_04132_14107# + G);  
            end if;  
        else  
            if I = 1 then  
                G := 8#0.62207_73250_42055_06043# +  
                        (8#0.62207_73250_42055_06043# - G);  
            else  
                G := -G;  
            end if;  
        end if;

----Return the result.

        if not In_Radians then  
            G := G / Degrees_To_Radians;  
        end if;  
        return G;

    end Acos;

--\f

    function Atan (A : Float) return Float is  
        F : Float;  
        G : Float;  
        R : Float;  
        N : Natural32 range 0 .. 3;  
    begin

----Reduce the argument to a small value close to zero.

        F := abs A;  
        if F > 1.0 then  
            F := 1.0 / F;  
            N := 2;  
        else  
            N := 0;  
        end if;  
        if F > 8#0.21114_12136_47546_52614# then  
            F := (((8#0.56663_65641_30231_25163_55# * F - 0.5) - 0.5) + F) /  
                    (8#1.56663_65641_30231_25163_55# + F);  
            N := N + 1;  
        end if;

----Approximate the result.

        if abs F >= Float (Float'Machine_Radix) **  
                       (-Float'Machine_Mantissa / 2) then  
            G := F * F;  
            if (Float'Machine_Radix = 2 and then  
                Float'Machine_Mantissa >= 51) or else  
               (Float'Machine_Radix = 4 and then  
                Float'Machine_Mantissa * 2 >= 51) or else  
               (Float'Machine_Radix = 8 and then  
                Float'Machine_Mantissa * 3 >= 51) or else  
               (Float'Machine_Radix = 16 and then  
                Float'Machine_Mantissa * 4 >= 51) or else  
               (Float'Machine_Radix = 10 and then  
                Float'Machine_Mantissa >= 16) then  
                R := (((-0.83758_29936_81500_59274E+00 * G -  
                        0.84946_24035_13206_83534E+01) * G -  
                       0.20505_85519_58616_51981E+02) * G -  
                      0.13688_76889_41919_26929E+02) * G  
                      / ((((G + 0.15024_00116_00285_76121E+02) * G +  
                          0.59578_43614_25973_44465E+02) * G +  
                         0.86157_34959_71302_42515E+02) * G +  
                        0.41066_30668_25757_81263E+02);  
            elsif (Float'Machine_Radix = 2 and then  
                   Float'Machine_Mantissa >= 33) or else  
                  (Float'Machine_Radix = 4 and then  
                   Float'Machine_Mantissa * 2 >= 33) or else  
                  (Float'Machine_Radix = 8 and then  
                   Float'Machine_Mantissa * 3 >= 33) or else  
                  (Float'Machine_Radix = 16 and then  
                   Float'Machine_Mantissa * 4 >= 33) or else  
                  (Float'Machine_Radix = 10 and then  
                   Float'Machine_Mantissa >= 11) then  
                R := ((-0.79439_12954_08336_251E+00 * G -  
                       0.42744_49853_67930_329E+01) * G -  
                      0.42743_26720_26241_096E+01) * G  
                      / (((G + 0.91978_93648_35039_806E+01) * G +  
                         0.20517_13765_64218_456E+02) * G +  
                        0.12822_98016_07919_841E+02);  
            elsif (Float'Machine_Radix = 2 and then  
                   Float'Machine_Mantissa >= 25) or else  
                  (Float'Machine_Radix = 4 and then  
                   Float'Machine_Mantissa * 2 >= 25) or else  
                  (Float'Machine_Radix = 8 and then  
                   Float'Machine_Mantissa * 3 >= 25) or else  
                  (Float'Machine_Radix = 16 and then  
                   Float'Machine_Mantissa * 4 >= 25) or else  
                  (Float'Machine_Radix = 10 and then  
                   Float'Machine_Mantissa >= 9) then  
                R := (-0.72002_68488_98E+00 * G -  
                      0.14400_83448_74E+01) * G  
                      / ((G + 0.47522_25845_99E+01) * G +  
                        0.43202_50389_19E+01);  
            else  
                R := (-0.50909_58253E-01 * G -  
                      0.47083_25141E+00) * G  
                      / (G + 0.14125_00740E+01);  
            end if;  
            F := F + F * R;  
        end if;

----Correct the result based upon the initial range reduction performed above.

        if N > 1 then  
            F := -F;  
        end if;  
        case N is  
            when 0 =>  
                null;  
            when 1 =>  
                F := F + 0.52359_87755_98298_87308;  
            when 2 =>  
                F := F + 1.57079_63267_94896_61923;  
            when 3 =>  
                F := F + 1.04719_75511_96597_74615;  
        end case;

----Return a result of the correct sign.

        if not In_Radians then  
            F := F / Degrees_To_Radians;  
        end if;  
        if A < 0.0 then  
            F := -F;  
        end if;  
        return F;

    end Atan;

--\f

    function Asec (A : Float) return Float is  
    begin  
        return Acos (1.0 / A);  
    end Asec;

--\f

    function Acsc (A : Float) return Float is  
    begin  
        return Asin (1.0 / A);  
    end Acsc;

--\f

    function Acot (A : Float) return Float is  
    begin  
        return Atan (1.0 / A);  
    end Acot;

--\f

    function Ln (A : Float) return Float is
------------------------------------------------------------------------------
-- Ln  - returns the natural logarithm of its positive argument
--
-- Based upon elements from
-- Software Manual for the Elementary Functions by W. J. Cody, Jr. and
-- W. Waite; Prentice-Hall 1980; pg 35-48
------------------------------------------------------------------------------
        F    : Float;  
        N    : Integer32;  
        Znum : Float;  
        Zden : Float;  
    begin

----Base=10 systems use the log10 function.

        if Float'Machine_Radix = 10 then  
            return Log10 (A) / Log10_E;  
        else

----Negative arguments are errors.  A 0.0 argument has an "easy" answer.

            if A < 0.0 then  
                raise Constraint_Error;  
            end if;  
            if A = 0.0 then  
                raise Numeric_Error;  
            end if;

----Extract the fraction and the exponent from the argument.  Adjust them
--  so that 0.5 <= f < 1.0.  We want the exponent expressed as a power of
--  2.0.

            Fraction_Exponent (A, F, N);  
            case Float'Machine_Radix is  
                when 2 =>  
                    null;  
                when 4 =>  
                    N := N * 2;  
                    while F < 0.5 loop  
                        F := F + F;  
                        N := N - 1;  
                    end loop;  
                when 8 =>  
                    N := N * 3;  
                    while F < 0.5 loop  
                        F := F + F;  
                        N := N - 1;  
                    end loop;  
                when 16 =>  
                    N := N * 4;  
                    while F < 0.5 loop  
                        F := F + F;  
                        N := N - 1;  
                    end loop;  
                when others =>  
                    raise Program_Error;  
            end case;

----Now we adjust N and z if f <= sqrt(0.5) and we just adjust z if
--  f > sqrt(0.5).

            if F > Sqrt_One_Half then  
                Znum := (F - 0.5) - 0.5;  
                Zden := F * 0.5 + 0.5;  
            else  
                N    := N - 1;  
                Znum := F - 0.5;  
                Zden := Znum * 0.5 + 0.5;  
            end if;

----Now we evaluate a rational approximation of exp(z).

            Znum := Znum / Zden;  
            F    := Znum * Znum;  
            if (Float'Machine_Radix = 2 and then  
                Float'Machine_Mantissa >= 49) or else  
               (Float'Machine_Radix = 4 and then  
                Float'Machine_Mantissa * 2 >= 49) or else  
               (Float'Machine_Radix = 8 and then  
                Float'Machine_Mantissa * 3 >= 49) or else  
               (Float'Machine_Radix = 16 and then  
                Float'Machine_Mantissa * 4 >= 49) then  
                Zden :=  
                   ((-0.78956_11288_74912_57267E+00 * F +  
                     0.16383_94356_30215_34222E+02)  
                     * F - 0.64124_94342_37455_81147E+02)  
                    / (((F - 0.35667_97773_90346_46171E+02)  
                        * F + 0.31203_22209_19245_32844E+03)  
                       * F - 0.76949_93210_84948_79777E+03);  
            elsif (Float'Machine_Radix = 2 and then  
                   Float'Machine_Mantissa >= 33) or else  
                  (Float'Machine_Radix = 4 and then  
                   Float'Machine_Mantissa * 2 >= 33) or else  
                  (Float'Machine_Radix = 8 and then  
                   Float'Machine_Mantissa * 3 >= 33) or else  
                  (Float'Machine_Radix = 16 and then  
                   Float'Machine_Mantissa * 4 >= 33) then  
                Zden :=  
                   ((F * 0.44445_51510_98033_23E-02 -  
                     0.63260_86623_38596_65E+00)  
                     * F + 0.37339_16896_31608_66E+01) /  
                   ((F - 0.14312_35435_58853_24E+02)  
                     * F + 0.44807_00275_57364_36E+02);  
            elsif (Float'Machine_Radix = 2 and then  
                   Float'Machine_Mantissa >= 25) or else  
                  (Float'Machine_Radix = 4 and then  
                   Float'Machine_Mantissa >= 25) or else  
                  (Float'Machine_Radix = 8 and then  
                   Float'Machine_Mantissa >= 25) or else  
                  (Float'Machine_Radix = 16 and then  
                   Float'Machine_Mantissa >= 25) then  
                Zden := (F * 0.13600_95468_621E-01 - 0.46490_62303_464E+00) /  
                           (F - 0.55788_73750_242E+01);  
            else  
                Zden := (-0.555270_74855E+00) / (F - 0.66327_18214E+01);  
            end if;  
            F := Znum + Znum * (F * Zden);

----Now compute the result.

            Znum := Float (N);  
            F    := (-2.12194_44005_46905_827679E-04 * Znum + F);  
            F    := F + Znum * 8#0.543#;  
            return F;

        end if;

    end Ln;

--\f

    function Log10 (A : Float) return Float is
------------------------------------------------------------------------------
-- Log10  - returns the base-10 logarithm of its positive argument
--
-- Based upon elements from
-- Software Manual for the Elementary Functions by W. J. Cody, Jr. and
-- W. Waite; Prentice-Hall 1980; pg 35-48
------------------------------------------------------------------------------
        F  : Float;  
        N  : Integer32;  
        S  : Float;  
        S2 : Float;  
    begin

----Base/=10 systems use the ln function.

        if Float'Machine_Radix /= 10 then  
            return Ln (A) / Ln_10;  
        else

----Negative arguments are errors.  A 0.0 argument has an "easy" answer.

            if A < 0.0 then  
                raise Constraint_Error;  
            end if;  
            if A = 0.0 then  
                raise Numeric_Error;  
            end if;

----Extract the fraction and the exponent from the argument.

            Fraction_Exponent (A, F, N);

----Now we adjust N and z if f <= sqrt(0.5) and we just adjust z if
--  f > sqrt(0.5).

            if F > Sqrt_One_Tenth then  
                null;  
            else  
                N := N - 1;  
                F := F * 10.0;  
            end if;

----Now we evaluate a rational approximation of exp(z).

            S  := ((F - 0.5) - 0.5) / (F + 1.0);  
            S2 := S * S;  
            if Float'Machine_Mantissa >= 16 then  
                F := ((((-0.74101_07841_61919_23924E+00 * S2 +  
                         0.10333_85715_14793_86456E+02)  
                         * S2 - 0.39273_74102_03156_25018E+02)  
                        * S2 + 0.55408_59120_41205_93097E+02)  
                       * S2 - 0.26044_70024_05557_63612E+02)  
                      / (((((S2 - 0.19373_23458_32854_78594E+2)  
                            * S2 + 0.10710_97891_15668_00857E+03)  
                           * S2 - 0.24430_30353_41829_54205E+03)  
                          * S2 + 0.24534_76188_68489_34790E+3)  
                         * S2 - 0.89955_20778_81033_11704E+02);  
            elsif Float'Machine_Mantissa >= 11 then  
                F := (((-0.71433_38215_32264_273E+00 * S2 +  
                        0.62503_65112_79083_731E+01)  
                        * S2 - 0.13682_37024_15026_896E+02)  
                       * S2 + 0.85167_31987_23885_403E+01)  
                      / ((((S2 - 0.13210_47835_01562_817E+02)  
                           * S2 + 0.47925_25604_38733_968E+02)  
                          * S2 - 0.64906_68274_09428_483E+02)  
                         * S2 + 0.29415_75017_23226_173E+02);  
            elsif Float'Machine_Mantissa >= 8 then  
                F := ((-0.67358_16014_777E+00 * S2 + 0.31630_34915_570E+01)  
                       * S2 - 0.29156_81437_901E+01)  
                      / (((S2 - 0.81908_00454_670E+01)  
                          * S2 + 0.16966_98140_210E+02)  
                         * S2 - 0.10070_40695_423E+02);  
            else  
                F := (-0.60368_24627E+00 * S2 + 0.10756_13712E+01) /  
                        ((S2 - 0.43144_78001E+01) * S2 + 0.37150_53570E+01);  
            end if;  
            S2 := S2 * F;  
            F  := S * (Log10_E + S2);

----Now compute the result.

            F := F + Float (N);  
            return F;

        end if;

    end Log10;

--\f

    function Logx (X : Float; Y : Float) return Float is  
    begin  
        if Float'Machine_Radix /= 10 then  
            return Ln (Y) / Ln (X);  
        else  
            return Log10 (Y) / Log10 (X);  
        end if;  
    end Logx;

--\f

    function Exp (A : Float) return Float is
------------------------------------------------------------------------------
-- Exp - returns E to the power represented by its argument
--
-- Based upon elements from Computer Approximations by John F. Hart et. al.
-- Krieger Publishing 1968
-- pg 97-105
------------------------------------------------------------------------------
        N  : Integer32;  
        Xn : Float;  
        G  : Float;  
        Xi : Integer32;  
        X1 : Float;  
        X2 : Float;  
        Pz : Float;  
        Qz : Float;  
    begin

----Make sure that it is possible to represent the result.  (too big/small)

        if A > Ln_Float_Large then  
            raise Numeric_Error;  
        end if;  
        if A < Ln_Float_Small then  
            return 0.0;  
        end if;

----If the answer will inevitably be 1.0 then return 1.0.

        if abs A <  
           Float (Float'Machine_Radix) ** (-Float'Machine_Mantissa) / 2.0 then  
            return 1.0;  
        end if;

----Handle non-base-10 machines here.

        if Float'Machine_Radix /= 10 then

----Perform agrument range reduction.

            N  := Integer32  
                     (A * 1.44269504088896340735992468100189213742664595602385);  
            Xn := Float (N);  
            Xi := Integer32 (abs A);  
            if Float (Xi) > abs A then  
                Xi := Xi - 1;  
            end if;  
            if A < 0.0 then  
                X1 := -Float (Xi);  
            else  
                X1 := Float (Xi);  
            end if;  
            X2 := A - X1;  
            G  := ((X1 - Xn * 8#0.543#) + X2) +  
                     Xn * 2.1219_44400_54690_58277E-04;

----Now approximate exp(g)

            X1 := G * G;  
            if (Float'Machine_Radix = 2 and then  
                Float'Machine_Mantissa >= 57) or else  
               (Float'Machine_Radix = 4 and then  
                Float'Machine_Mantissa * 2 >= 57) or else  
               (Float'Machine_Radix = 8 and then  
                Float'Machine_Mantissa * 3 >= 57) or else  
               (Float'Machine_Radix = 16 and then  
                Float'Machine_Mantissa * 4 >= 57) then  
                Pz := (0.31555_19276_56846_46356E-04 * X1 +  
                       0.75753_18015_94227_76666E-02) * X1 + 0.25;  
                Qz := ((0.75104_02839_98700_46114E-06 * X1 +  
                        0.63121_89437_43985_03557E-03) * X1 +  
                       0.56817_30269_85512_21787E-01) * X1 + 0.5;  
            elsif (Float'Machine_Radix = 2 and then  
                   Float'Machine_Mantissa >= 43) or else  
                  (Float'Machine_Radix = 4 and then  
                   Float'Machine_Mantissa * 2 >= 43) or else  
                  (Float'Machine_Radix = 8 and then  
                   Float'Machine_Mantissa * 3 >= 43) or else  
                  (Float'Machine_Radix = 16 and then  
                   Float'Machine_Mantissa * 4 >= 43) then  
                Pz := (0.16520_33002_68279_130E-04 * X1 +  
                       0.69436_00015_11792_852E-02) * X1 +  
                      0.24999_99999_99999_993E+00;  
                Qz := (0.49586_28849_05441_294E-03 * X1 +  
                       0.55553_86669_69001_118E-01) * X1 + 0.5;  
            elsif (Float'Machine_Radix = 2 and then  
                   Float'Machine_Mantissa >= 30) or else  
                  (Float'Machine_Radix = 4 and then  
                   Float'Machine_Mantissa * 2 >= 30) or else  
                  (Float'Machine_Radix = 8 and then  
                   Float'Machine_Mantissa * 3 >= 30) or else  
                  (Float'Machine_Radix = 16 and then  
                   Float'Machine_Mantissa * 4 >= 30) then  
                Pz := 0.59504_25497_7591E-02 * X1 + 0.24999_99999_9992E+00;  
                Qz :=  
                   (0.29729_36368_2238E-03 * X1 + 0.53567_51764_522E-01) * X1 +  
                      0.5;  
            else  
                Pz := 0.41602_88626_8E-02 * X1 + 0.24999_99995_0E+00;  
                Qz := 0.49987_17877_8E-01 * X1 + 0.5;  
            end if;  
            Pz := Pz * G;  
            G  := 0.5 + Pz / (Qz - Pz);

----Compute the final result.

            N  := N + 1;  
            Xi := N / 2;  
            N  := N - Xi;  
            G  := G * 2.0 ** Integer (Xi);  
            G  := G * 2.0 ** Integer (N);  
            return G;

----Handle base-10 machines here.

        else

----Perform agrument range reduction.

            N := Integer32  
                    (A * 0.868588963806503655302257837833210164588794011255454);  
            Xn := Float (N);  
            Xi := Integer32 (abs A);  
            if Float (Xi) > abs A then  
                Xi := Xi - 1;  
            end if;  
            if A < 0.0 then  
                X1 := -Float (Xi);  
            else  
                X1 := Float (Xi);  
            end if;  
            X2 := A - X1;  
            G := ((X1 - Xn * 1.151E+00) + X2) - Xn * 2.9254_64970_22842_009E-04;

----Now approximate exp(g)

            X1 := G * G;  
            if Float'Machine_Mantissa >= 15 then  
                Pz := (0.42041_42681_37450_31524E+02 * X2 +  
                       0.10097_41487_24273_91798E+05) * X2 +  
                      0.33326_70292_26801_61137E+06;  
                Qz := ((X2 + 0.84124_35845_14154_54514E+03) * X2 +  
                       0.75739_33461_59883_44395E+05) * X2 +  
                      0.66653_40584_53603_22323E+06;  
            elsif Float'Machine_Mantissa >= 11 then  
                Pz := (0.33287_36465_16410_3E-01 * X1 +  
                       0.14008_29975_62819_6E+02) * X1 +  
                      0.50446_48895_05869_7E+03;  
                Qz := (X1 + 0.11209_40810_96616_9E+03) * X1 +  
                         0.50445_48895_05869_7E+03;  
            elsif Float'Machine_Mantissa >= 8 then  
                Pz := 0.20041_43275_526E+02 * X1 + 0.84249_03867_900E+03;  
                Qz := (X2 + 0.18049_79288_462E+03) * X2 + 0.16849_80773_608E+04;  
            else  
                Pz := 0.83046_5413E-01 * X1 + 0.50034_9857E+01;  
                Qz := X1 + 0.10006_9975E+02;  
            end if;  
            Pz := Pz * G;  
            G  := 0.5 + Pz / (Qz - Pz);  
            G  := G + G;

----Compute the final result.

            if N rem 2 /= 0 then  
                if N > 0 then  
                    G :=  
                       G * 3.16227766016837933199889354443271853371955514008093;  
                else  
                    G :=  
                       G / 3.16227766016837933199889354443271853371955514008093;  
                end if;  
            end if;  
            N  := N / 2;  
            Xi := N / 2;  
            N  := N - Xi;  
            G  := G * 2.0 ** Integer (Xi);  
            G  := G * 2.0 ** Integer (N);  
            return G;

        end if;

    end Exp;

--\f

    function Exp10 (A : Float) return Float is  
    begin  
        return Expx (10.0, A);  
    end Exp10;

--\f

    function Expx (X : Float; Y : Float) return Float is  
        M   : Integer32;  
        F   : Float;  
        R   : Integer32;  
        P   : Integer32;  
        Z   : Float;  
        Rv  : Float;  
        U1  : Float;  
        U2  : Float;  
        Y1  : Float;  
        Y2  : Float;  
        Iw1 : Integer32;  
        W   : Float;  
        W1  : Float;  
        W2  : Float;  
        I   : Integer32;  
        N   : Integer32;  
    begin

----We use the conventions where (-f)**g is illegal and f**0.0 is 1.0.

        if Y = 0.0 then  
            return 1.0;  
        elsif X < 0.0 then  
            raise Constraint_Error;  
        end if;

----Here we handle non-base-10 floating point machines.

        if Float'Machine_Radix /= 10 then  
            Fraction_Exponent (X, F, M);  
            case Float'Machine_Radix is  
                when 2 =>  
                    null;  
                when 4 =>  
                    if F < 0.5 then  
                        R := 1;  
                        F := F + F;  
                    else  
                        R := 0;  
                    end if;  
                when 8 =>  
                    if F < 0.25 then  
                        R := 2;  
                        F := F + F;  
                        F := F + F;  
                    elsif F < 0.5 then  
                        R := 1;  
                        F := F + F;  
                    else  
                        R := 0;  
                    end if;  
                when 16 =>  
                    if F < 0.125 then  
                        R := 3;  
                        F := F + F;  
                        F := F + F;  
                        F := F + F;  
                    elsif F < 0.25 then  
                        R := 2;  
                        F := F + F;  
                        F := F + F;  
                    elsif F < 0.5 then  
                        R := 1;  
                        F := F + F;  
                    else  
                        R := 0;  
                    end if;  
                when others =>  
                    raise Program_Error;  
            end case;

----Compute the P index.

            P := 1;  
            if F <= Expx_A1 (9) then  
                P := 9;  
            end if;  
            if F <= Expx_A1 (P + 4) then  
                P := P + 4;  
            end if;  
            if F <= Expx_A1 (P + 2) then  
                P := P + 2;  
            end if;

----Compute z.

            if Float'Machine_Radix /= 2 then  
                Z := F - Expx_A1 (P + 1);  
                Z := Z - Expx_A2 ((P + 1) / 2);  
                Z := Z / (F * 0.5 + Expx_A1 (P + 1) * 0.5);  
            else  
                Z := F - Expx_A1 (P + 1);  
                Z := Z - Expx_A2 ((P + 1) / 2);  
                Z := Z / (F + Expx_A1 (P + 1));  
                Z := Z + Z;  
            end if;

----Now compute the z*P(z**2) approximation.

            F := Z * Z;  
            if (Float'Machine_Radix = 2 and then  
                Float'Machine_Mantissa >= 51) or else  
               (Float'Machine_Radix = 4 and then  
                Float'Machine_Mantissa * 2 >= 51) or else  
               (Float'Machine_Radix = 8 and then  
                Float'Machine_Mantissa * 3 >= 51) or else  
               (Float'Machine_Radix = 16 and then  
                Float'Machine_Mantissa * 4 >= 51) then  
                Rv := ((0.43445_77567_21631_19635E-03 * F +  
                        0.22321_42128_59242_58967E-02) * F +  
                       0.12500_00000_05037_99174E-01) * F +  
                      0.83333_33333_33332_11405E-01;  
            elsif (Float'Machine_Radix = 2 and then  
                   Float'Machine_Mantissa >= 37) or else  
                  (Float'Machine_Radix = 4 and then  
                   Float'Machine_Mantissa * 2 >= 37) or else  
                  (Float'Machine_Radix = 8 and then  
                   Float'Machine_Mantissa * 3 >= 37) or else  
                  (Float'Machine_Radix = 16 and then  
                   Float'Machine_Mantissa * 4 >= 37) then  
                Rv :=  
                   (0.22338_24352_81541_8E-02 * F + 0.12499_99796_50060_8E-01) *  
                      F +  
                   0.83333_33333_41213_6E-01;  
            elsif (Float'Machine_Radix = 2 and then  
                   Float'Machine_Mantissa >= 25) or else  
                  (Float'Machine_Radix = 4 and then  
                   Float'Machine_Mantissa * 2 >= 25) or else  
                  (Float'Machine_Radix = 8 and then  
                   Float'Machine_Mantissa * 3 >= 25) or else  
                  (Float'Machine_Radix = 16 and then  
                   Float'Machine_Mantissa * 4 >= 25) then  
                Rv := 0.12506_48500_52E-01 * F + 0.83333_32862_45E-01;  
            else  
                Rv := 0.83357_541E-01;  
            end if;  
            Rv := Z * F * Rv;  
            Rv := Rv + 0.44269_50408_88963_40736 * Rv;  
            U2 := (Rv + Z * 0.44269_50408_88963_40736) + Z;

----Now create W, our partial answer.

            case Float'Machine_Radix is  
                when 2 =>  
                    U1 := Float (M * 16 - P) * 0.0625;  
                when 4 =>  
                    U1 := Float ((2 * M - R) * 16 - P) * 0.0625;  
                when 8 =>  
                    U1 := Float ((3 * M - R) * 16 - P) * 0.0625;  
                when 16 =>  
                    U1 := Float ((4 * M - R) * 16 - P) * 0.0625;  
                when others =>  
                    raise Program_Error;  
            end case;  
            Y1  := Truncate (Y * 16.0) * 0.0625;  
            Y2  := Y - Y1;  
            W   := U2 * Y + U1 * Y2;  
            W1  := Truncate (W * 16.0) * 0.0625;  
            W2  := W - W1;  
            W   := W1 + U1 * Y1;  
            W1  := Truncate (W * 16.0) * 0.0625;  
            W2  := W2 + (W - W1);  
            W   := Truncate (W2 * 16.0) * 0.0625;  
            Iw1 := Truncate (16.0 * (W1 + W));  
            W2  := W2 - W;  
            if W2 > 0.0 then  
                Iw1 := Iw1 + 1;  
                W2  := W2 - 1.0 / 16.0;  
            end if;  
            if Iw1 < 0 then  
                I := 0;  
            else  
                I := 1;  
            end if;

----Compute the new p,r,m values.

            case Float'Machine_Radix is  
                when 2 =>  
                    M := Iw1 / 16 + I;  
                    P := 16 * M - Iw1;  
                when 4 =>  
                    N := Iw1 / 16 + I;  
                    P := 16 * N - Iw1;  
                    M := N / 2 + I;  
                    R := 2 * M - N;  
                when 8 =>  
                    N := Iw1 / 16 + I;  
                    P := 16 * N - Iw1;  
                    M := N / 3 + I;  
                    R := 3 * M - N;  
                when 16 =>  
                    N := Iw1 / 16 + I;  
                    P := 16 * N - Iw1;  
                    M := N / 4 + I;  
                    R := 4 * M - N;  
                when others =>  
                    raise Program_Error;  
            end case;

----Now do Z.

            if (Float'Machine_Radix = 2 and then  
                Float'Machine_Mantissa >= 53) or else  
               (Float'Machine_Radix = 4 and then  
                Float'Machine_Mantissa * 2 >= 53) or else  
               (Float'Machine_Radix = 8 and then  
                Float'Machine_Mantissa * 3 >= 53) or else  
               (Float'Machine_Radix = 16 and then  
                Float'Machine_Mantissa * 4 >= 53) then  
                Z := ((((((0.14928_85268_05956_08186E-04 * W2 +  
                           0.15400_29044_09897_64601E-03) * W2 +  
                          0.13333_54131_35857_84703E-02) * W2 +  
                         0.96181_29059_51724_16964E-02) * W2 +  
                        0.55504_10866_40855_95326E-01) * W2 +  
                       0.24022_65069_59095_37056E+00) * W2 +  
                      0.69314_71805_59945_29629E+00) * W2;  
            elsif (Float'Machine_Radix = 2 and then  
                   Float'Machine_Mantissa >= 46) or else  
                  (Float'Machine_Radix = 4 and then  
                   Float'Machine_Mantissa * 2 >= 46) or else  
                  (Float'Machine_Radix = 8 and then  
                   Float'Machine_Mantissa * 3 >= 46) or else  
                  (Float'Machine_Radix = 16 and then  
                   Float'Machine_Mantissa * 4 >= 46) then  
                Z := (((((0.15077_40617_88142_382E-03 * W2 +  
                          0.13330_81011_34082_075E-02) * W2 +  
                         0.96181_17691_38724_104E-02) * W2 +  
                        0.55504_10842_47568_661E-01) * W2 +  
                       0.24022_65069_56777_522E+00) * W2 +  
                      0.69314_71805_59937_815E+00) * W2;  
            elsif (Float'Machine_Radix = 2 and then  
                   Float'Machine_Mantissa >= 25) or else  
                  (Float'Machine_Radix = 4 and then  
                   Float'Machine_Mantissa * 2 >= 25) or else  
                  (Float'Machine_Radix = 8 and then  
                   Float'Machine_Mantissa * 3 >= 25) or else  
                  (Float'Machine_Radix = 16 and then  
                   Float'Machine_Mantissa * 4 >= 25) then  
                Z := ((((0.13052_55159_42810E-02 * W2 +  
                         0.96162_06595_83789E-02) * W2 +  
                        0.55504_04881_30765E-01) * W2 +  
                       0.24022_65061_44710E+00) * W2 +  
                      0.69314_71805_56341E+00) * W2;  
            else  
                Z := ((0.54360_383E-01 * W2 +  
                       0.24018_510E+00) * W2 +  
                      0.69314_675E+00) * W2;  
            end if;

----Nearly there; one more step.

            Z := Expx_A1 (P + 1) + Expx_A1 (P + 1) * Z;  
            if Float'Machine_Radix /= 2 then  
                Z := Z * 2.0 ** Integer (R);  
            end if;  
            I := M / 2;  
            M := M - I;  
            Z := Z * Float (Float'Machine_Radix) ** Integer (I);  
            Z := Z * Float (Float'Machine_Radix) ** Integer (M);  
            return Z;

----Here we handle base-10 floating point machines.

        else  
            Fraction_Exponent (X, F, M);

----Compute the P index.

            P := 1;  
            while F <= Expx_A1 (P + 2) loop  
                P := P + 2;  
            end loop;

----Compute z.

            Z := F - Expx_A1 (P + 1);  
            Z := Z - Expx_A2 ((P + 1) / 2);  
            Z := Z / (F + Expx_A1 (P + 1));

----Now compute the z*P(z**2) approximation.

            F := Z * Z;  
            if Float'Machine_Mantissa >= 15 then  
                Rv := ((-0.68650_74087_13682_76474E+00 * F +  
                        0.35829_52156_67236_56966E+01) * F -  
                       0.35184_51552_93711_75192E+01) * F  
                       / (((F - 0.89628_58940_90517_72560E+01) * F +  
                          0.19666_45902_38905_35938E+02) * F -  
                         0.12152_30114_43221_36016E+02);  
            elsif Float'Machine_Mantissa >= 9 then  
                Rv := (-0.62077_57571_87918_2E+00 * F +  
                       0.12095_10556_17031_7E+01) * F  
                       / ((F - 0.46505_84419_61834_9E+01) * F +  
                         0.41775_01464_74140_2E+01);  
            else  
                Rv := -0.48144_78516E+00 * F / (F - 0.16628_73334E+01);  
            end if;  
            U2 := Z * (Rv + 0.86858_89638_06503_65530);

----Now create W, our partial answer.

            U1  := Float (M * 16 - P) * 0.0625;  
            Y1  := Truncate (Y * 10.0) * 0.1;  
            Y2  := Y - Y1;  
            W   := U2 * Y + U1 * Y2;  
            W1  := Truncate (W * 10.0) * 0.1;  
            W2  := W - W1;  
            W   := W1 + U1 * Y1;  
            W1  := Truncate (W * 10.0) * 0.1;  
            W2  := W2 + (W - W1);  
            W   := Truncate (W2 * 10.0) * 0.1;  
            Iw1 := Truncate (10.0 * (W1 + W));  
            W2  := W2 - W;  
            if Iw1 > 10 * Float'Machine_Emax - 1 then  
                raise Numeric_Error;                -- Result far too large
            end if;  
            if Iw1 < 10 * Float'Machine_Emin + 1 then  
                return 0.0;                         -- Result far too small
            end if;  
            if W2 > 0.0 then  
                Iw1 := Iw1 + 1;  
                W2  := W2 - 0.1;  
            end if;  
            if Iw1 < 0 then  
                I := 0;  
            else  
                I := 1;  
            end if;

----Compute the new p,m values.

            M := Iw1 / 10 + I;  
            P := 10 * M - Iw1;

----Now do Z.

            if Float'Machine_Mantissa >= 15 then  
                Z := ((((0.41071_57590_70777_89734E+0 * W2 +  
                         0.35616_93889_90563_55495E+01) * W2 +  
                        0.34363_74379_23829_31895E+02) * W2 +  
                       0.28195_17958_42955_14217E+02) * W2 +  
                      0.23539_70807_49261_19272E+03) * W2  
                      / ((((W2 - 0.10224_89787_60959_10327E+02) * W2 +  
                          0.45994_87860_33907_54883E+02) * W2 -  
                         0.10545_35294_64899_79088E+03) * W2 +  
                        0.10223_16532_25538_32434E+03);  
            elsif Float'Machine_Mantissa >= 11 then  
                Z := (((-0.51359_30940_15587_9E+0 * W2 -  
                        0.46710_07431_07040_0E+01) * W2 -  
                       0.57975_91698_28398_6E+01) * W2 -  
                      0.37713_99188_31947_3E+02) * W2  
                      / (((W2 - 0.63664_36505_45209_5E+01) * W2 +  
                         0.16339_13385_89750_1E+02) * W2 -  
                        0.16378_97856_54144_4E+02);  
            elsif Float'Machine_Mantissa >= 8 then  
                Z := ((0.68536_85048_84E+00 * W2 +  
                       0.17805_77775_84E+01) * W2 +  
                      0.82998_14853_06E+01) * W2  
                      / ((W2 - 0.33766_11203_26E+01) * W2 +  
                        0.36045_63799_45E+01);  
            else  
                Z := ((-0.78964_7014E+00 * W2 -  
                       0.21742_4847E+01) * W2 -  
                      0.38883_9832E+01) * W2  
                      / (W2 - 0.16887_0925E+01);  
            end if;

----Nearly there; one more step.

            Z := Expx_A1 (P + 1) + Expx_A1 (P + 1) * Z;  
            I := M / 2;  
            M := M - I;  
            Z := Z * 10.0 ** Integer (I);  
            Z := Z * 10.0 ** Integer (M);  
            return Z;

        end if;

    end Expx;

--\f

    function Sinh (A : Float) return Float is
------------------------------------------------------------------------------
-- Sinh - returns the hyperbolic sine of its argument
------------------------------------------------------------------------------
        Y : Float := abs A;  
        W : Float;  
        Z : Float;  
    begin

----SINH for "large" arguments is computed this way.

        if Y > 1.0 then

----Base-2 machines can compute things the easy way for "small" arguments.

            if Float'Machine_Radix = 2 and then Y <= Ln_Float_Large then  
                Z := Exp (Y);  
                Z := (Z * 0.5 - 0.5 / Z);  
                if A < 0.0 then  
                    return -Z;  
                else  
                    return Z;  
                end if;  
            end if;

----Reduce the argument.

            if Float'Machine_Radix /= 10 then  
                W := Y - 8#0.542714#;  
            else  
                W := Y - 0.6932;  
            end if;  
            if W > Ln_Float_Large then  
                raise Numeric_Error;  
            end if;

----Now take the EXP of the argument and compute the result with corrections
--  based on the machine radix.

            Z := Exp (W);  
            if (Float'Machine_Radix = 10 and then  
                W <= 1.16 * Float (Float'Machine_Mantissa + 1)) then  
                Z := Z - 0.24997_35916_74870_15965 / Z;  
            elsif (Float'Machine_Radix /= 10 and then  
                   Float'Machine_Radix /= 2 and then  
                   W <= 0.35 * Float (Float'Machine_Mantissa + 1)) then  
                Z := Z - 0.24999_30850_04514_99336 / Z;  
            end if;  
            if Float'Machine_Radix /= 10 then  
                Z := Z + 0.13830_27787_96019_02638E-04 * Z;  
            else  
                Z := Z + 0.52820_83502_58748_52469E-04 * Z;  
            end if;

----Return the result with the correct sign.

            if A < 0.0 then  
                return -Z;  
            else  
                return Z;  
            end if;

----Compute the result for "small" arguments here.

        else  
            if Y < Float (Float'Machine_Radix) **  
                      (-Float'Machine_Mantissa / 2) then  
                return A;  
            end if;  
            Y := A * A;  
            if (Float'Machine_Radix = 2 and then  
                Float'Machine_Mantissa >= 51) or else  
               (Float'Machine_Radix = 4 and then  
                Float'Machine_Mantissa * 2 >= 51) or else  
               (Float'Machine_Radix = 8 and then  
                Float'Machine_Mantissa * 3 >= 51) or else  
               (Float'Machine_Radix = 16 and then  
                Float'Machine_Mantissa * 4 >= 51) or else  
               (Float'Machine_Radix = 10 and then  
                Float'Machine_Mantissa >= 16) then  
                W := (((-0.78966_12741_73570_99479E+00 * Y -  
                        0.16375_79820_26307_51372E+03) * Y -  
                       0.11563_52119_68517_68270E+05) * Y -  
                      0.35181_28343_01771_17881E+06) * Y  
                      / (((Y - 0.27773_52311_96507_01667E+03) * Y +  
                         0.36162_72310_94218_36460E+05) * Y -  
                        0.21108_77005_81062_71242E+07);  
            elsif (Float'Machine_Radix = 2 and then  
                   Float'Machine_Mantissa >= 41) or else  
                  (Float'Machine_Radix = 4 and then  
                   Float'Machine_Mantissa * 2 >= 41) or else  
                  (Float'Machine_Radix = 8 and then  
                   Float'Machine_Mantissa * 3 >= 41) or else  
                  (Float'Machine_Radix = 16 and then  
                   Float'Machine_Mantissa * 4 >= 41) or else  
                  (Float'Machine_Radix = 10 and then  
                   Float'Machine_Mantissa >= 13) then  
                W := (((0.77239_39820_29419_23E-02 * Y +  
                        0.13286_42866_92242_29E+01) * Y +  
                       0.85943_28483_85490_10E+02) * Y +  
                      0.23941_43592_30500_69E+04) * Y  
                      / ((Y - 0.20258_33686_64278_69E+03) * Y +  
                        0.14364_86155_38302_92E+05);  
            elsif (Float'Machine_Radix = 2 and then  
                   Float'Machine_Mantissa >= 25) or else  
                  (Float'Machine_Radix = 4 and then  
                   Float'Machine_Mantissa * 2 >= 25) or else  
                  (Float'Machine_Radix = 8 and then  
                   Float'Machine_Mantissa * 3 >= 25) or else  
                  (Float'Machine_Radix = 16 and then  
                   Float'Machine_Mantissa * 4 >= 25) or else  
                  (Float'Machine_Radix = 10 and then  
                   Float'Machine_Mantissa >= 8) then  
                W := ((0.34364_14035_8506E+00 * Y +  
                       0.31359_75645_6058E+02) * Y +  
                      0.10622_28883_7151E+04) * Y  
                      / ((Y - 0.13051_01250_9199E+03) * Y +  
                        0.63733_73302_1822E+04);  
            else  
                W := (-0.19033_3399E+00 * Y - 0.71379_3159E+01) * Y  
                         / (Y - 0.42827_7109E+02);  
            end if;  
            W := A + A * W;  
            return W;  
        end if;

    end Sinh;

--\f

    function Cosh (A : Float) return Float is
------------------------------------------------------------------------------
--  Cosh - returns the hyperbolic cosine of its argument
------------------------------------------------------------------------------
        Y : Float := abs A;  
        W : Float;  
        Z : Float;  
    begin

----Base-2 machines can compute things the easy way for "small" arguments.

        if (Float'Machine_Radix = 2 and then Y <= Ln_Float_Large) or else  
           (Float'Machine_Radix /= 2 and then Y <= 1.0) then  
            Z := Exp (Y);  
            Z := (Z * 0.5 + 0.5 / Z);  
            return Z;  
        end if;

----Reduce the argument.

        if Float'Machine_Radix /= 10 then  
            W := Y - 8#0.542714#;  
        else  
            W := Y - 0.6932;  
        end if;  
        if W > Ln_Float_Large then  
            raise Numeric_Error;  
        end if;

----Now take the EXP of the argument and compute the result with corrections
--  based on the machine radix.

        Z := Exp (W);  
        if (Float'Machine_Radix = 10 and then  
            W <= 1.16 * Float (Float'Machine_Mantissa + 1)) then  
            Z := Z + 0.24997_35916_74870_15965 / Z;  
        elsif (Float'Machine_Radix /= 10 and then  
               Float'Machine_Radix /= 2 and then  
               W <= 0.35 * Float (Float'Machine_Mantissa + 1)) then  
            Z := Z + 0.24999_30850_04514_99336 / Z;  
        end if;  
        if Float'Machine_Radix /= 10 then  
            Z := Z + 0.13830_27787_96019_02638E-04 * Z;  
        else  
            Z := Z + 0.52820_83502_58748_52469E-04 * Z;  
        end if;

----Return the result.

        return Z;

    end Cosh;

--\f

    function Tanh (A : Float) return Float is
------------------------------------------------------------------------------
--  Tanh - returns the hyperbolic tangent of its argument
------------------------------------------------------------------------------
        F      : Float := abs A;  
        G      : Float;  
        Result : Float;  
    begin

----Bigger arguments tend toward 1.0.

        if F > Tanh_Xbig then  
            Result := 1.0;

----Middling big arguments can be best approximated this way.

        elsif F > 0.54930_61443_34054_84570 then  
            Result := 0.5 - 1.0 / (Exp (F + F) + 1.0);  
            Result := Result + Result;

----Very small arguments result in themselves.

        elsif F < Float (Float'Machine_Radix) **  
                     (-Float'Machine_Mantissa / 2) then  
            Result := F;

----All other arguments get approximated here.

        else  
            G := F * F;  
            if (Float'Machine_Radix = 2 and then  
                Float'Machine_Mantissa >= 49) or else  
               (Float'Machine_Radix = 4 and then  
                Float'Machine_Mantissa * 2 >= 49) or else  
               (Float'Machine_Radix = 8 and then  
                Float'Machine_Mantissa * 3 >= 49) or else  
               (Float'Machine_Radix = 16 and then  
                Float'Machine_Mantissa * 4 >= 49) or else  
               (Float'Machine_Radix = 10 and then  
                Float'Machine_Mantissa >= 15) then  
                Result := ((-0.96437_49277_72254_69787E+00 * G -  
                            0.99225_92967_22360_83313E+02) * G -  
                           0.16134_11902_39962_28058E+04) * G  
                           / (((G + 0.11274_47438_05349_49335E+03) * G +  
                              0.22337_72071_89623_12926E+04) * G +  
                             0.48402_35707_19886_88686E+04);  
            elsif (Float'Machine_Radix = 2 and then  
                   Float'Machine_Mantissa >= 37) or else  
                  (Float'Machine_Radix = 4 and then  
                   Float'Machine_Mantissa * 2 >= 37) or else  
                  (Float'Machine_Radix = 8 and then  
                   Float'Machine_Mantissa * 3 >= 37) or else  
                  (Float'Machine_Radix = 16 and then  
                   Float'Machine_Mantissa * 4 >= 37) or else  
                  (Float'Machine_Radix = 10 and then  
                   Float'Machine_Mantissa >= 12) then  
                Result := ((-0.36242_42193_46421_73E-03 * G -  
                            0.92318_68945_14261_77E+00) * G -  
                           0.19059_52242_69822_92E+02) * G  
                           / ((G + 0.25640_98759_51789_75E+02) * G +  
                             0.57178_56728_09658_17E+02);  
            elsif (Float'Machine_Radix = 2 and then  
                   Float'Machine_Mantissa >= 25) or else  
                  (Float'Machine_Radix = 4 and then  
                   Float'Machine_Mantissa * 2 >= 25) or else  
                  (Float'Machine_Radix = 8 and then  
                   Float'Machine_Mantissa * 3 >= 25) or else  
                  (Float'Machine_Radix = 16 and then  
                   Float'Machine_Mantissa * 4 >= 25) or else  
                  (Float'Machine_Radix = 10 and then  
                   Float'Machine_Mantissa >= 9) then  
                Result := (-0.93363_47565_2401E+00 * G -  
                           0.21063_95800_0245E+02)  
                           / ((G + 0.28077_65347_0471E+02) * G +  
                             0.63191_87401_5582E+02);  
            else  
                Result := (-0.38310_10665E-02 * G - 0.82377_28127E+00) * G  
                              / (G + 0.24713_19654E+01);  
            end if;  
            Result := F + F * Result;  
        end if;

----Return a result with the correct sign.

        if A < 0.0 then  
            return -Result;  
        else  
            return Result;  
        end if;

    end Tanh;

--\f

    function Sech (A : Float) return Float is  
    begin  
        return 1.0 / Cosh (A);  
    end Sech;

--\f

    function Csch (A : Float) return Float is  
    begin  
        return 1.0 / Sinh (A);  
    end Csch;

--\f

    function Coth (A : Float) return Float is  
    begin  
        return 1.0 / Tanh (A);  
    end Coth;

--\f

    function Asinh (A : Float) return Float is  
    begin  
        return Ln (A + Sqrt (A * A + 1.0));  
    end Asinh;

--\f

    function Acosh (A : Float) return Float is  
    begin  
        return Ln (A + Sqrt (A * A - 1.0));  
    end Acosh;

--\f

    function Atanh (A : Float) return Float is  
    begin  
        return Ln ((1.0 + A) / (1.0 - A)) / 2.0;  
    end Atanh;

--\f

    function Asech (A : Float) return Float is  
    begin  
        return Acosh (1.0 / A);  
    end Asech;

--\f

    function Acsch (A : Float) return Float is  
    begin  
        return Asinh (1.0 / A);  
    end Acsch;

--\f

    function Acoth (A : Float) return Float is  
    begin  
        return Atanh (1.0 / A);  
    end Acoth;

--\f

    function Hex (Dig : Integer32) return Character is
------------------------------------------------------------------------------
-- Returns '0'..'9' or 'A'..'F'; Dig in 0..15.
------------------------------------------------------------------------------
    begin  
        if Dig <= 9 then  
            return Character'Val (Dig + Character'Pos ('0'));  
        else  
            return Character'Val (Dig - 10 + Character'Pos ('A'));  
        end if;  
    end Hex;

--\f

    function Cat (Frac : Float; Bits : Natural32) return String is
------------------------------------------------------------------------------
-- Convert a fraction that is in the range 1.0/Float'Machine_Radix <= X < 1.0
-- into a bunch of Hex digits with an implicit initial decimal point.
-- This routine is only called when Float'Machine_Radix /= 10
------------------------------------------------------------------------------
        Valstr : String (1 .. 4);  
        I      : Integer32;  
        Nxt    : Float;  
    begin

----Get the next 16 bits off of the top of the fraction.

        Nxt := Frac * 16.0 ** 4;  
        I   := Integer32 (Nxt);  
        if Float (I) > Nxt then  
            I := I - 1;  
        end if;

----Remove those 16 bits from the fraction.  Do it this "strange" way so that
--  we don't get wobbling-precision problems on base-4/8/16 machines.

        Nxt := Frac - Float (I) / 16.0 ** 4;  
        Nxt := Nxt * 16.0 ** 4;

----Convert the 4 hex digits to a string.

        Valstr (1) := Hex (I / 16 ** 3);  
        Valstr (2) := Hex (I / 16 ** 2 rem 16);  
        Valstr (3) := Hex (I / 16 rem 16);  
        Valstr (4) := Hex (I rem 16);

----Get the rest of the digits and return.

        if Bits > 16 then  
            return Valstr & Cat (Nxt, Bits - 16);  
        else  
            return Valstr;  
        end if;

    end Cat;

--\f

--/ if R1000 then

    function Debug_Tools_Float_Image  
                (Value           : Float;  
                 Level           : Natural;  
                 Prefix          : String;  
                 Expand_Pointers : Boolean) return String is
------------------------------------------------------------------------------
-- Return a string containing the normal base-10 printed representation of the
-- string and (if Machine_Radix /= 10) the Machine_Radix representation.
------------------------------------------------------------------------------
        Frac : Float;  
        Expo : Integer32;  
        Str  : String (1 .. Float'Digits + 12);  
        Bits : Natural32;  
    begin

----Get the base 10 string.

        Float_Text_Io.Put (Str, Value, Aft => Float'Digits + 3, Exp => 5);

----Compute the number of bits in the mantissa.

        case Float'Machine_Radix is  
            when 2 =>  
                Bits := Float'Machine_Mantissa;  
            when 4 =>  
                Bits := Float'Machine_Mantissa * 2;  
            when 8 =>  
                Bits := Float'Machine_Mantissa * 3;  
            when 16 =>  
                Bits := Float'Machine_Mantissa * 4;  
            when others =>  
                return Str;  
        end case;

----Get the fraction and the exponent.

        Fraction_Exponent (Value, Frac, Expo);

----Create the string containing the Machine_Radix version of the string and
--  return it plus the base-10 string.

        return Str & "  0."  
                   & Cat (Frac, Bits)  
                   & " *" & Natural32'Image (Float'Machine_Radix)  
                   & " **" & Integer32'Image (Expo);

    end Debug_Tools_Float_Image;

    procedure Register_Float is  
       new Debug_Tools.Register (Float,  
                                 Debug_Tools_Float_Image);

--/ end if; -- R1000

--\f

begin

----Register our preferred debug image routine.

--/ if R1000 then
    Register_Float;
--/ end if;

----Ln of the largest and of the smallest numbers that exist.

    Ln_Float_Large := Ln (Float'Last);  
    Ln_Float_Small := Ln (2.0 ** (Float'Machine_Emin - 1));

----Used in argument reduciton in Tanh.

    Tanh_Xbig :=  
       (Ln_2 +  
        Float (Float'Machine_Mantissa + 1) * Ln (Float (Float'Machine_Radix)))  
        / 2.0;

----Cot_Eps1 is smallest positive number that can be reciprocated without
--  Numeric_Error.

    if Float'Machine_Emax < -Float'Machine_Emin then  
        Cot_Eps1 := Float (Float'Machine_Radix) ** (-Float'Machine_Emax + 1);  
    else  
        Cot_Eps1 := Float (Float'Machine_Radix) ** (Float'Machine_Emin - 1);  
    end if;

----Set up arrays of values for the Expx routine.  We must do this dynamicly
--  because this package is intended to be portable across architectures and
--  Ada provides us with no way to do this staticly.  See pages 98 and 106 of
--  the William J. Cody book for details of what and why.

    declare  
        type Natural32_Array is array (Positive range <>) of Natural32;  
        procedure Expx_Stuff (A1 : in out Float;  
                              A2 : in out Float;  
                              N1 :        Natural32_Array) is
            ----N1/2/3 represent a floating point value.  There are 3*30=90
            --  bits of fraction in there.  We take those 90 bits and we
            --  place the first Float'Machine_Mantissa worth of them into A1
            --  and we stick "the rest" into A2.  No way to do this at
            --  compile time.
            Bits  : Integer32;  
            Mask  : Integer32;  
            Shift : Float;  
            Tmp   : Integer;  
        begin  
            if Float'Machine_Radix /= 10 then
                ----Handle all 2**N-base floating points here.
                case Float'Machine_Radix is  
                    when 2 =>  
                        Bits := Float'Machine_Mantissa;  
                    when 4 =>  
                        Bits := Float'Machine_Mantissa * 2;  
                    when 8 =>  
                        Bits := Float'Machine_Mantissa * 3;  
                    when 16 =>  
                        Bits := Float'Machine_Mantissa * 4;  
                    when others =>  
                        raise Program_Error;  
                end case;  
                if Bits >= 60 then  
                    A1   := Float (N1 (1)) / 2.0 ** 30 +  
                               Float (N1 (2)) / 2.0 ** 60;  
                    Mask := 2 ** Integer (90 - Bits);  
                    A1   := A1 + Float (N1 (3) - N1 (3) rem Mask) / 2.0 ** 90;  
                    A2   := Float (N1 (3) rem Mask) / 2.0 ** 90;  
                elsif Bits >= 30 then  
                    A1   := Float (N1 (1)) / 2.0 ** 30;  
                    Mask := 2 ** Integer (60 - Bits);  
                    A1   := A1 + Float (N1 (2) - N1 (2) rem Mask) / 2.0 ** 60;  
                    A2   := Float (N1 (2) rem Mask) / 2.0 ** 60 +  
                               Float (N1 (3)) / 2.0 ** 90;  
                else  
                    Mask := 2 ** Integer (30 - Bits);  
                    A1   := Float (N1 (1) - N1 (1) rem Mask) / 2.0 ** 30;  
                    A2   := Float (N1 (1) rem Mask) / 2.0 ** 30 +  
                               Float (N1 (2)) / 2.0 ** 60;  
                end if;  
            else
                ----Handle base-10 floating-point here.
                Shift := 1.0 / 10.0 ** 5;  
                Bits  := 0;  
                A1    := 0.0;  
                for I in 1 .. 6 loop  
                    A1    := A1 + Float (N1 (I)) * Shift;  
                    Shift := Shift / 10.0 ** 5;  
                    Bits  := Bits + 5;  
                    if Bits + 5 > Float'Machine_Mantissa then  
                        Tmp := I + 1;  
                        goto Do_A2;  
                    end if;  
                end loop;  
                A2 := 0.0;  
                return;  
                <<Do_A2>> null;  
                Bits := Float'Machine_Mantissa - Bits;  
                Mask := 10 ** Integer (5 - Bits);  
                A1   := A1 + Float (N1 (Tmp) - N1 (Tmp) rem Mask) * Shift;  
                A2   := Float (N1 (Tmp) rem Mask) * Shift;  
                for I in Tmp + 1 .. 6 loop  
                    A2    := A2 + Float (N1 (I)) * Shift;  
                    Shift := Shift / 10.0 ** 5;  
                end loop;  
            end if;  
        end Expx_Stuff;  
    begin  
        if Float'Machine_Radix /= 10 then  
            Expx_A1 (1)  := 1.0;  
            Expx_A1 (3)  := 8#0.72540_30671_75644_41622_73201_32727#;  
            Expx_A1 (5)  := 8#0.65642_37462_55323_53257_20715_15057#;  
            Expx_A1 (7)  := 8#0.61263_45204_25240_66655_64761_25503#;  
            Expx_A1 (9)  := 8#0.55202_36314_77473_63110_21313_73047#;  
            Expx_A1 (11) := 8#0.51377_32652_33052_11616_50345_72776#;  
            Expx_A1 (13) := 8#0.46033_76024_30667_05226_75065_32214#;  
            Expx_A1 (15) := 8#0.42712_70170_76521_36556_33737_10612#;  
            Expx_A1 (17) := 8#0.4#;

            Expx_Stuff (Expx_A1 (2), Expx_A2 (1),  
                        (8#75222_57505#, 8#22220_66302#, 8#61734_72062#));  
            Expx_Stuff (Expx_A1 (4), Expx_A2 (2),  
                        (8#70146_33673#, 8#02522_47021#, 8#04062_61124#));  
            Expx_Stuff (Expx_A1 (6), Expx_A2 (3),  
                        (8#63422_21405#, 8#21760_44016#, 8#17421_53016#));  
            Expx_Stuff (Expx_A1 (8), Expx_A2 (4),  
                        (8#57204_24347#, 8#65401_41553#, 8#72504_02177#));  
            Expx_Stuff (Expx_A1 (10), Expx_A2 (5),  
                        (8#53254_07672#, 8#44124_12254#, 8#31114_01243#));  
            Expx_Stuff (Expx_A1 (12), Expx_A2 (6),  
                        (8#47572_46230#, 8#11064_10432#, 8#66404_42174#));  
            Expx_Stuff (Expx_A1 (14), Expx_A2 (7),  
                        (8#44341_72334#, 8#72542_16063#, 8#30176_55544#));  
            Expx_Stuff (Expx_A1 (16), Expx_A2 (8),  
                        (8#41325_30331#, 8#74611_03661#, 8#23056_22556#));  
        else  
            Expx_A1 (1)  := 1.0;  
            Expx_A1 (3)  := 0.63095_73444_80193_24943_43601_366;  
            Expx_A1 (5)  := 0.39810_71705_53497_25077_02523_051;  
            Expx_A1 (7)  := 0.25118_86431_50958_01110_85032_068;  
            Expx_A1 (9)  := 0.15848_93192_46111_34852_02101_373;  
            Expx_A1 (11) := 0.1;

            Expx_Stuff (Expx_A1 (2), Expx_A2 (1),  
                        (79432, 82347, 24281, 50206, 59182, 82800));  
            Expx_Stuff (Expx_A1 (4), Expx_A2 (2),  
                        (50118, 72336, 27272, 28500, 15541, 86900));  
            Expx_Stuff (Expx_A1 (6), Expx_A2 (3),  
                        (31622, 77660, 16837, 93319, 98893, 54400));  
            Expx_Stuff (Expx_A1 (8), Expx_A2 (4),  
                        (19952, 62314, 96887, 96013, 52455, 39700));  
            Expx_Stuff (Expx_A1 (10), Expx_A2 (5),  
                        (12589, 25411, 79416, 72104, 23954, 10600));  
        end if;  
    end;

----If this is not true then this module will be producing values that are
--  accurate to at-most 60-bits/18-digits.  There aren't too many machines
--  with larger floating-point than that but we definitely won't function
--  properly (ie. with the desired accuracy).  If you don't care about this
--  then comment out the Assert.

    if ((Float'Machine_Radix = 2 and then Float'Machine_Mantissa <= 60) or else  
        (Float'Machine_Radix = 4 and then  
         Float'Machine_Mantissa * 2 <= 60) or else  
        (Float'Machine_Radix = 8 and then  
         Float'Machine_Mantissa * 3 <= 60) or else  
        (Float'Machine_Radix = 16 and then  
         Float'Machine_Mantissa * 4 <= 60) or else  
        (Float'Machine_Radix = 10 and then Float'Machine_Mantissa <= 18)) then  
        null;  
    else  
        raise Program_Error;  
    end if;

end Trig;