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

⟦fed769ff7⟧ TextFile

    Length: 16052 (0x3eb4)
    Types: TextFile
    Names: »B«

Derivation

└─⟦a7d1ea751⟧ Bits:30000550 8mm tape, Rational 1000, !users!projects 94_04_11
    └─ ⟦129cab021⟧ »DATA« 
        └─⟦9b477e385⟧ 
└─⟦f64eaa120⟧ Bits:30000752 8mm tape, Rational 1000, !projects 93 02 16
    └─ ⟦6f12a12be⟧ »DATA« 
        └─⟦9b477e385⟧ 
└─⟦2f6cfab89⟧ Bits:30000547 8mm tape, Rational 1000, !projects 94-01-04
    └─ ⟦d65440be7⟧ »DATA« 
        └─⟦9b477e385⟧ 
            └─⟦this⟧ 

TextFile

with Text_Io;
use Text_Io;
with Floating_Characteristics;
use Floating_Characteristics;
with Numeric_Io;
use Numeric_Io;
with Numeric_Primitives;
use Numeric_Primitives;

package body Core_Functions is

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

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

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

    function Sqrt (X : Float) return Float is
        M, N : Exponent_Type;
        F, Y : Mantissa_Type;
        Result : Float;
        subtype Index is Integer range 0 ..
                                          100;    --  #########################
        Sqrt_L1 : Index := 3;
        --  Could get away with SQRT_L1 := 2 for 28 bits
        --  Using the better Cody-Waite coeficients overflows MANTISSA_TYPE
        Sqrt_C1 : Mantissa_Type := 8#0.3317777777#;
        Sqrt_C2 : Mantissa_Type := 8#0.4460000000#;
        Sqrt_C3 : Mantissa_Type := 8#0.55202_36314_77747_36311_0#;
    begin
        if X = Zero then
            Result := Zero;
            return Result;
        elsif X = One then
            --  To get exact SQRT(1.0)
            Result := One;
            return Result;
        elsif X < Zero then
            New_Line;
            Put ("CALLED SQRT FOR NEGATIVE ARGUMENT   ");
            raise Program_Error;
            --Put (X);
            Put ("   USED ABSOLUTE VALUE");
            New_Line;
            Result := Sqrt (abs (X));
            return Result;
        else
            Defloat (X, N, F);
            Y := Sqrt_C1 + Mantissa_Type (Sqrt_C2 * F);

            for J in 1 .. Sqrt_L1 loop
                Y := Y / Mantissa_Divisor_2 + Mantissa_Type
                                                 ((F / Mantissa_Divisor_2) / Y);
            end loop;

            if (N mod 2) /= 0 then
                Y := Mantissa_Type (Sqrt_C3 * Y);
                N := N + 1;
            end if;

            M := N / 2;
            Refloat (M, Y, Result);
            return Result;
        end if;

    exception
        when others =>
            New_Line;
            Put (" EXCEPTION IN SQRT, X = ");
            raise Program_Error;

            --Put (X);
            Put ("  RETURNED 1.0");
            New_Line;
            return One;
    end Sqrt;

    function Cbrt (X : Float) return Float is
        M, N : Exponent_Type;
        F, Y : Mantissa_Type;
        Result : Float;
        subtype Index is Integer range 0 ..
                                          100;    --  #########################
        Cbrt_L1 : Index := 3;
        Cbrt_C1 : Mantissa_Type := 0.5874009;
        Cbrt_C2 : Mantissa_Type := 0.4125990;
        Cbrt_C3 : Mantissa_Type := 0.62996_05249;
        Cbrt_C4 : Mantissa_Type := 0.79370_05260;
    begin
        if X = Zero then
            Result := Zero;
            return Result;
        else
            Defloat (X, N, F);
            F := abs (F);
            Y := Cbrt_C1 + Mantissa_Type (Cbrt_C2 * F);

            for J in 1 .. Cbrt_L1 loop
                Y := Y - (Y / Mantissa_Divisor_3 -
                          Mantissa_Type ((F / Mantissa_Divisor_3) /
                                         Mantissa_Type (Y * Y)));
            end loop;

            case (N mod 3) is
                when 0 =>
                    null;
                when 1 =>
                    Y := Mantissa_Type (Cbrt_C3 * Y);
                    N := N + 2;
                when 2 =>
                    Y := Mantissa_Type (Cbrt_C4 * Y);
                    N := N + 1;
                when others =>
                    null;
            end case;

            M := N / 3;

            if X < Zero then
                Y := -Y;
            end if;

            Refloat (M, Y, Result);
            return Result;
        end if;

    exception
        when others =>
            Result := One;

            if X < Zero then
                Result := -One;
            end if;

            New_Line;
            Put ("EXCEPTION IN CBRT, X = ");
            raise Program_Error;
            --Put (X);
            Put ("  RETURNED  ");
            raise Program_Error;
            --Put (Result);
            New_Line;
            return Result;
    end Cbrt;

    function Log (X : Float) return Float is
        --  Uses fixed formulation for generality
        Result : Float;
        N : Exponent_Type;
        Xn : Float;
        Y : Float;
        F : Mantissa_Type;
        Z, Zden, Znum : Mantissa_Type;
        C0 : constant Mantissa_Type := 0.20710_67811_86547_52440;
        --  SQRT(0.5) - 0.5
        C1 : constant Float := 8#0.543#;
        C2 : constant Float := -2.12194_44005_46905_82767_9E-4;

        function R (Z : Mantissa_Type) return Mantissa_Type is
            --  Use fixed formulation here because the float coeficents are > 1.0
            --  and would exceed the limits on a MANTISSA_TYPE
            A0 : constant Mantissa_Type := 0.04862_85276_587;
            B0 : constant Mantissa_Type := 0.69735_92187_803;
            B1 : constant Mantissa_Type := -0.125;
            C : constant Mantissa_Type := 0.01360_09546_862;
        begin
            return Z + Mantissa_Type
                          (Z * Mantissa_Type
                                  (Mantissa_Type (Z * Z) *
                                   (C + Mantissa_Type
                                           (A0 / (B0 +
                                                  Mantissa_Type
                                                     (B1 * Mantissa_Type
                                                              (Z * Z)))))));
        end R;
    begin
        if X < Zero then
            New_Line;
            Put ("CALLED LOG FOR NEGATIVE ");
            raise Program_Error; --put (X);;
            Put ("   USE ABS => ");
            Result := Log (abs (X));
            raise Program_Error;
            --Put (Result);
            New_Line;
        elsif X = Zero then
            New_Line;
            Put ("CALLED LOG FOR ZERO ARGUMENT, RETURNED ");
            Result := -Xmax;      --  SUPPOSED TO BE -LARGE
            raise Program_Error;--Put (Result);
            New_Line;
        else
            Defloat (X, N, F);
            Znum := F - Mantissa_Half;
            Y := Convert_To_Float (Znum);
            Zden := Znum / Mantissa_Divisor_2 + Mantissa_Half;

            if Znum > C0 then
                Y := Y - Mantissa_Half;
                Znum := Znum - Mantissa_Half;
                Zden := Zden + Mantissa_Half / Mantissa_Divisor_2;
            else
                N := N - 1;
            end if;

            Z := Mantissa_Type (Znum / Zden);
            Result := Convert_To_Float (R (Z));

            if N /= 0 then
                Xn := Convert_To_Float (N);
                Result := (Xn * C2 + Result) + Xn * C1;
            end if;
        end if;

        return Result;

    exception
        when others =>
            New_Line;
            Put (" EXCEPTION IN LOG, X = ");
            raise Program_Error; --put (X);;
            Put ("  RETURNED 0.0");
            New_Line;
            return Zero;
    end Log;

    function Log10 (X : Float) return Float is
        Log_10_Of_2 : constant Float :=
           Convert_To_Float (Mantissa_Type (

                                            8#0.33626_75425_11562_41615#));
    begin
        return Log (X) * Log_10_Of_2;
    end Log10;

    function Exp (X : Float) return Float is
        Result : Float;
        N : Exponent_Type;
        Xg, Xn, X1, X2 : Float;
        F, G : Mantissa_Type;
        Bigx : Float := Exp_Large;
        Smallx : Float := Exp_Small;
        One_Over_Log_2 : constant Float := 1.4426_95040_88896_34074;
        C1 : constant Float := 0.69335_9375;
        C2 : constant Float := -2.1219_44400_54690_58277E-4;

        function R (G : Mantissa_Type) return Mantissa_Type is
            Z, Gp, Q : Mantissa_Type;
            P0 : constant Mantissa_Type := 0.24999_99999_9992;
            P1 : constant Mantissa_Type := 0.00595_04254_9776;
            Q0 : constant Mantissa_Type := 0.5;
            Q1 : constant Mantissa_Type := 0.05356_75176_4522;
            Q2 : constant Mantissa_Type := 0.00029_72936_3682;
        begin
            Z := Mantissa_Type (G * G);
            Gp := Mantissa_Type ((Mantissa_Type (P1 * Z) + P0) * G);
            Q := Mantissa_Type ((Mantissa_Type (Q2 * Z) + Q1) * Z) + Q0;
            return Mantissa_Half + Mantissa_Type (Gp / (Q - Gp));
        end R;
    begin
        if X > Bigx then
            New_Line;
            Put ("  EXP CALLED WITH TOO BIG A POSITIVE ARGUMENT, ");
            raise Program_Error; --put (X);;
            Put ("   RETURNED XMAX");
            New_Line;
            Result := Xmax;
        elsif X < Smallx then
            New_Line;
            Put ("  EXP CALLED WITH TOO BIG A NEGATIVE ARGUMENT,  ");
            raise Program_Error; --put (X);;
            Put ("    RETURNED ZERO");
            New_Line;
            Result := Zero;
        elsif abs (X) < Eps then
            Result := One;
        else
            N := Exponent_Type (X * One_Over_Log_2);
            Xn := Convert_To_Float (N);
            X1 := Round (X);
            X2 := X - X1;
            Xg := ((X1 - Xn * C1) + X2) - Xn * C2;
            G := Mantissa_Type (Xg);
            N := N + 1;
            F := R (G);
            Refloat (N, F, Result);
        end if;

        return Result;

    exception
        when others =>
            New_Line;
            Put (" EXCEPTION IN EXP, X = ");
            raise Program_Error; --put (X);;
            Put ("  RETURNED 1.0");
            New_Line;
            return One;
    end Exp;

    function "**" (X, Y : Float) return Float is
        --  This is the last function to be coded since it appeared that it really
        --  was un-Ada-like and ought not be in the regular package
        --  Nevertheless it was included in this version
        --  It is specific for FLOAT and does not have the MANTISSA_TYPE generality
        M, N : Exponent_Type;
        G : Mantissa_Type;
        P, Temp, Iw1, I : Integer;
        Result, Z, V, R, U1, U2, W, W1, W2, W3, Y1, Y2 : Float;
        K : constant Float := 0.44269_50408_88963_40736;
        Ibigx : constant Integer := Integer
                                       (Truncate (16.0 * Log (Xmax) - 1.0));
        Ismallx : constant Integer := Integer
                                         (Truncate (16.0 * Log (Xmin) + 1.0));
        P1 : constant Float := 0.83333_32862_45E-1;
        P2 : constant Float := 0.12506_48500_52E-1;
        Q1 : constant Float := 0.69314_71805_56341;
        Q2 : constant Float := 0.24022_65061_44710;
        Q3 : constant Float := 0.55504_04881_30765E-1;
        Q4 : constant Float := 0.96162_06595_83789E-2;
        Q5 : constant Float := 0.13052_55159_42810E-2;
        A1 : array (1 .. 17) of Float :=
           (8#1.00000_0000#, 8#0.75222_5750#, 8#0.72540_3067#, 8#0.70146_3367#,
            8#0.65642_3746#, 8#0.63422_2140#, 8#0.61263_4520#, 8#0.57204_2434#,
            8#0.55202_3631#, 8#0.53254_0767#, 8#0.51377_3265#,
            8#0.47572_4623#, 8#0.46033_7602#, 8#0.44341_7233#,
            8#0.42712_7017#, 8#0.41325_3033#, 8#0.40000_0000#);
        A2 : array (1 .. 8) of Float :=
           (8#0.00000_00005_22220_66302_61734_72062#,
            8#0.00000_00003_02522_47021_04062_61124#,
            8#0.00000_00005_21760_44016_17421_53016#,
            8#0.00000_00007_65401_41553_72504_02177#,
            8#0.00000_00002_44124_12254_31114_01243#,
            8#0.00000_00000_11064_10432_66404_42174#,
            8#0.00000_00004_72542_16063_30176_55544#,
            8#0.00000_00001_74611_03661_23056_22556#);

        function Reduce (V : Float) return Float is
        begin
            return Float (Integer (16.0 * V)) * 0.0625;
        end Reduce;
    begin
        if X <= Zero then
            if X < Zero then
                Result := (abs (X)) ** Y;
                New_Line;
                Put ("X**Y CALLED WITH X = ");
                raise Program_Error; --put (X);
                New_Line;
                Put ("USED ABS, RETURNED ");
                raise Program_Error;--Put (Result);
                New_Line;
            else
                if Y <= Zero then
                    if Y = Zero then
                        Result := Zero;
                    else
                        Result := Xmax;
                    end if;

                    New_Line;
                    Put ("X**Y CALLED WITH X = 0, Y = ");
                    raise Program_Error; --put (Y);
                    New_Line;
                    Put ("RETURNED ");
                    raise Program_Error;--Put (Result);
                    New_Line;
                else
                    Result := Zero;
                end if;
            end if;
        else
            Defloat (X, M, G);
            P := 1;

            if G <= A1 (9) then
                P := 9;
            end if;

            if G <= A1 (P + 4) then
                P := P + 4;
            end if;

            if G <= A1 (P + 2) then
                P := P + 2;
            end if;

            Z := ((G - A1 (P + 1)) - A2 ((P + 1) / 2)) / (G + A1 (P + 1));
            Z := Z + Z;
            V := Z * Z;
            R := (P2 * V + P1) * V * Z;
            R := R + K * R;
            U2 := (R + Z * K) + Z;
            U1 := Float (Integer (M) * 16 - P) * 0.0625;
            Y1 := Reduce (Y);
            Y2 := Y - Y1;
            W := U2 * Y + U1 * Y2;
            W1 := Reduce (W);
            W2 := W - W1;
            W := W1 + U1 * Y1;
            W1 := Reduce (W);
            W2 := W2 + (W - W1);
            W3 := Reduce (W2);
            Iw1 := Integer (Truncate (16.0 * (W1 + W3)));
            W2 := W2 - W3;

            if W > Float (Ibigx) then
                Result := Xmax;
                Put ("X**Y CALLED  X =");
                raise Program_Error; --put (X);
                Put ("   Y =");
                raise Program_Error; --put (Y);
                Put ("   TOO LARGE  RETURNED ");
                raise Program_Error;--Put (Result);
                New_Line;
            elsif W < Float (Ismallx) then
                Result := Zero;
                Put ("X**Y CALLED  X =");
                raise Program_Error;
                --put(X);
                Put ("   Y =");
                raise Program_Error;
                Put ("   TOO SMALL  RETURNED ");
                raise Program_Error;--Put (Result);
                New_Line;
            else
                if W2 > Zero then
                    W2 := W2 - 0.0625;
                    Iw1 := Iw1 + 1;
                end if;

                if Iw1 < Integer (Zero) then
                    I := 0;
                else
                    I := 1;
                end if;

                M := Exponent_Type (I + Iw1 / 16);
                P := 16 * Integer (M) - Iw1;
                Z := ((((Q5 * W2 + Q4) * W2 + Q3) * W2 + Q2) * W2 + Q1) * W2;
                Z := A1 (P + 1) + (A1 (P + 1) * Z);
                Refloat (M, Z, Result);
            end if;
        end if;

        return Result;
    end "**";
begin
    Exp_Large := Log (Xmax) * (One - Eps);
    Exp_Small := Log (Xmin) * (One - Eps);
end Core_Functions;