DataMuseum.dk Presents historical artifacts from the history of: Rational R1000/400

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 Excavated with: AutoArchaeologist - Free & Open Source Software.

top - download

⟦1329b5ea7⟧ Ada Source

    Length: 26624 (0x6800)
    Types: Ada Source
    Notes: 03_class, FILE, R1k_Segment, e3_tag, package body Trig_Lib, seg_00eb17

Derivation

└─⟦8527c1e9b⟧ Bits:30000544 8mm tape, Rational 1000, Arrival backup of disks in PAM's R1000
    └─ ⟦cfc2e13cd⟧ »Space Info Vol 2« 
        └─⟦this⟧ 

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

package body Trig_Lib is

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

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

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

    function Sin (X : Float) return Float is
        Sgn, Y       : Float;
        N            : Integer;
        Xn           : Float;
        F, G, X1, X2 : Float;
        Result       : Float;
        Ymax         : Float          := Float (Integer (Pi * Two ** (It / 2)));
        Beta         : Float          := Convert_To_Float (Ibeta);
        Epsilon      : Float          := Beta ** (-It / 2);
        C1           : constant Float := 3.140625;
        C2           : constant Float := 9.6765_35897_93E-4;

        function R (G : Float) return Float is
            R1 : constant Float := -0.16666_66660_883;
            R2 : constant Float := 0.83333_30720_556E-2;
            R3 : constant Float := -0.19840_83282_313E-3;
            R4 : constant Float := 0.27523_97106_775E-5;
            R5 : constant Float := -0.23868_34640_601E-7;
        begin
            return ((((R5 * G + R4) * G + R3) * G + R2) * G + R1) * G;
        end R;
    begin
        if X < Zero then
            Sgn := -One;
            Y   := -X;
        else
            Sgn := One;
            Y   := X;
        end if;

        if Y > Ymax then
            New_Line;
            Put (" SIN CALLED WITH ARGUMENT TOO LARGE FOR ACCURACY ");
            New_Line;
        end if;

        N  := Integer (Y * One_Over_Pi);
        Xn := Convert_To_Float (N);

        if N mod 2 /= 0 then
            Sgn := -Sgn;
        end if;

        X1 := Truncate (abs (X));
        X2 := abs (X) - X1;
        F  := ((X1 - Xn * C1) + X2) - Xn * C2;

        if abs (F) < Epsilon then
            Result := F;
        else
            G      := F * F;
            Result := F + F * R (G);
        end if;

        return (Sgn * Result);
    end Sin;

    function Cos (X : Float) return Float is
        Sgn, Y       : Float;
        N            : Integer;
        Xn           : Float;
        F, G, X1, X2 : Float;
        Result       : Float;
        Ymax         : Float          := Float (Integer (Pi * Two ** (It / 2)));
        Beta         : Float          := Convert_To_Float (Ibeta);
        Epsilon      : Float          := Beta ** (-It / 2);
        C1           : constant Float := 3.140625;
        C2           : constant Float := 9.6765_35897_93E-4;

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

        if Y > Ymax then
            New_Line;
            Put (" COS CALLED WITH ARGUMENT TOO LARGE FOR ACCURACY ");
            New_Line;
        end if;

        N  := Integer (Y * One_Over_Pi);
        Xn := Convert_To_Float (N);

        if N mod 2 /= 0 then
            Sgn := -Sgn;
        end if;

        Xn := Xn - 0.5;          -- TO FORM COS INSTEAD OF SIN
        X1 := Truncate (abs (X));
        X2 := abs (X) - X1;
        F  := ((X1 - Xn * C1) + X2) - Xn * C2;

        if abs (F) < Epsilon then
            Result := F;
        else
            G      := F * F;
            Result := F + F * R (G);
        end if;

        return (Sgn * Result);
    end Cos;

    function Tan (X : Float) return Float is
        Sgn, Y       : Float;
        N            : Integer;
        Xn           : Float;
        F, G, X1, X2 : Float;
        Result       : Float;
        Ymax         : Float := Float (Integer (Pi * Two ** (It / 2))) / 2.0;
        Beta         : Float := Convert_To_Float (Ibeta);
        Epsilon      : Float := Beta ** (-It / 2);
        C1           : constant Float := 8#1.444#;
        C2           : constant Float := 4.8382_67948_97E-4;

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

        if Y > Ymax then
            New_Line;
            Put (" TAN CALLED WITH ARGUMENT TOO LARGE FOR ACCURACY ");
            New_Line;
        end if;

        N  := Integer (X * Two_Over_Pi);
        Xn := Convert_To_Float (N);
        X1 := Truncate (X);
        X2 := X - X1;
        F  := ((X1 - Xn * C1) + X2) - Xn * C2;

        if abs (F) < Epsilon then
            Result := F;
        else
            G      := F * F;
            Result := R (G);
        end if;

        if N mod 2 = 0 then
            return Result;
        else
            return -1.0 / Result;
        end if;
    end Tan;

    function Cot (X : Float) return Float is
        Sgn, Y       : Float;
        N            : Integer;
        Xn           : Float;
        F, G, X1, X2 : Float;
        Result       : Float;
        Ymax         : Float := Float (Integer (Pi * Two ** (It / 2))) / 2.0;
        Beta         : Float := Convert_To_Float (Ibeta);
        Epsilon      : Float := Beta ** (-It / 2);
        Epsilon1     : Float := 1.0 / Xmax;
        C1           : constant Float := 8#1.444#;
        C2           : constant Float := 4.8382_67948_97E-4;

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

        if Y < Epsilon1 then
            New_Line;
            Put (" COT CALLED WITH ARGUMENT TOO NEAR ZERO ");
            New_Line;

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

        if Y > Ymax then
            New_Line;
            Put (" COT CALLED WITH ARGUMENT TOO LARGE FOR ACCURACY ");
            New_Line;
        end if;

        N  := Integer (X * Two_Over_Pi);
        Xn := Convert_To_Float (N);
        X1 := Truncate (X);
        X2 := X - X1;
        F  := ((X1 - Xn * C1) + X2) - Xn * C2;

        if abs (F) < Epsilon then
            Result := F;
        else
            G      := F * F;
            Result := R (G);
        end if;

        if N mod 2 /= 0 then
            return -Result;
        else
            return 1.0 / Result;
        end if;
    end Cot;

    function Asin (X : Float) return Float is
        G, Y    : Float;
        Result  : Float;
        Beta    : Float := Convert_To_Float (Ibeta);
        Epsilon : Float := Beta ** (-It / 2);

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

    function Acos (X : Float) return Float is
        G, Y    : Float;
        Result  : Float;
        Beta    : Float := Convert_To_Float (Ibeta);
        Epsilon : Float := Beta ** (-It / 2);

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

    function Atan (X : Float) return Float is
        F, G : Float;
        subtype Region is Integer range 0 .. 3;    --  ##########
        N                : Region;
        Result           : Float;
        Beta             : Float          := Convert_To_Float (Ibeta);
        Epsilon          : Float          := Beta ** (-It / 2);
        Sqrt_3           : constant Float := 1.73205_08075_68877_29353;
        Sqrt_3_Minus_1   : constant Float := 0.73205_08075_68877_29353;
        Two_Minus_Sqrt_3 : constant Float := 0.26794_91924_31122_70647;

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

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

        if F > Two_Minus_Sqrt_3 then
            F := (((Sqrt_3_Minus_1 * F - 0.5) - 0.5) + F) / (Sqrt_3 + F);
            N := N + 1;
        end if;

        if abs (F) < Epsilon then
            Result := F;
        else
            G      := F * F;
            Result := F + F * R (G);
        end if;

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

        case N is
            when 0 =>
                Result := Result;
            when 1 =>
                Result := Pi_Over_Six + Result;
            when 2 =>
                Result := Pi_Over_Two + Result;
            when 3 =>
                Result := Pi_Over_Three + Result;
        end case;

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

        return Result;
    end Atan;

    function Atan2 (V, U : Float) return Float is
        X, Result : Float;
    begin
        if U = 0.0 then
            if V = 0.0 then
                Result := 0.0;
                New_Line;
                Put (" ATAN2 CALLED WITH 0/0   RETURNED ");
                New_Line;
            elsif V > 0.0 then
                Result := Pi_Over_Two;
            else
                Result := -Pi_Over_Two;
            end if;
        else
            X := abs (V / U);
            --  If underflow or overflow is detected, go to the exception
            Result := Atan (X);

            if U < 0.0 then
                Result := Pi - Result;
            end if;

            if V < 0.0 then
                Result := -Result;
            end if;
        end if;

        return Result;

    exception
        when Numeric_Error =>
            if abs (V) > abs (U) then
                Result := Pi_Over_Two;

                if V < 0.0 then
                    Result := -Result;
                end if;
            else
                Result := 0.0;

                if U < 0.0 then
                    Result := Pi - Result;
                end if;
            end if;

            return Result;
    end Atan2;

    function Sinh (X : Float) return Float is
    begin
        return X;
    end Sinh;

    function Cosh (X : Float) return Float is
    begin
        return X;
    end Cosh;

    function Tanh (X : Float) return Float is
    begin
        return X;
    end Tanh;
begin
    null;
end Trig_Lib;

    nblk1=19
    nid=0
    hdr6=32
        [0x00] rec0=1b rec1=00 rec2=01 rec3=022
        [0x01] rec0=16 rec1=00 rec2=02 rec3=008
        [0x02] rec0=03 rec1=00 rec2=19 rec3=006
        [0x03] rec0=28 rec1=00 rec2=03 rec3=010
        [0x04] rec0=00 rec1=00 rec2=17 rec3=016
        [0x05] rec0=17 rec1=00 rec2=18 rec3=028
        [0x06] rec0=04 rec1=00 rec2=04 rec3=046
        [0x07] rec0=21 rec1=00 rec2=05 rec3=016
        [0x08] rec0=01 rec1=00 rec2=16 rec3=04e
        [0x09] rec0=21 rec1=00 rec2=06 rec3=01c
        [0x0a] rec0=00 rec1=00 rec2=15 rec3=00e
        [0x0b] rec0=19 rec1=00 rec2=07 rec3=062
        [0x0c] rec0=03 rec1=00 rec2=14 rec3=002
        [0x0d] rec0=27 rec1=00 rec2=08 rec3=012
        [0x0e] rec0=00 rec1=00 rec2=13 rec3=00e
        [0x0f] rec0=1d rec1=00 rec2=09 rec3=016
        [0x10] rec0=00 rec1=00 rec2=12 rec3=01c
        [0x11] rec0=19 rec1=00 rec2=0a rec3=012
        [0x12] rec0=01 rec1=00 rec2=11 rec3=00e
        [0x13] rec0=1d rec1=00 rec2=0b rec3=018
        [0x14] rec0=00 rec1=00 rec2=10 rec3=018
        [0x15] rec0=27 rec1=00 rec2=0c rec3=010
        [0x16] rec0=00 rec1=00 rec2=0f rec3=00a
        [0x17] rec0=29 rec1=00 rec2=0d rec3=00c
        [0x18] rec0=0e rec1=00 rec2=0e rec3=001
    tail 0x2170b5d0c8225429ec9bf 0x42a00088462063203