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

⟦e13771806⟧ TextFile

    Length: 12287 (0x2fff)
    Types: TextFile
    Names: »B«

Derivation

└─⟦5f3412b64⟧ Bits:30000745 8mm tape, Rational 1000, ENVIRONMENT 12_6_5 TOOLS 
    └─ ⟦91c658230⟧ »DATA« 
        └─⟦458657fb6⟧ 
            └─⟦1472c4407⟧ 
                └─⟦this⟧ 
└─⟦d10a02448⟧ Bits:30000409 8mm tape, Rational 1000, ENVIRONMENT, D_12_7_3
    └─ ⟦fc9b38f02⟧ »DATA« 
        └─⟦9b46a407a⟧ 
            └─⟦2e03b931c⟧ 
                └─⟦this⟧ 

TextFile

-- This package body contains an implementation of the proposed standard for
-- primitive real arithmetic functions, as defined by the
-- ISO-IEC/JTC1/SC22/WG9 (Ada) committee (draft 1.0, December 1990)
--
-- Users recognize that the Math Support packages provided as part of the
-- Rational Environment are delivered "AS IS". No warranties are given,
-- whether express, implied, or statutory, including implied warranties of
-- fitness for a particular purpose or mechantability. In no event will
-- Rational be liable in tort, negligence, or other liability incurred as a
-- result of the use of this Math Support packages.
--
-- In addition, Rational does not guarantee conformance to the proposed
-- standards, especially in the accuracy of the results of the various math
-- functions.
--
-- This code can be distributed and modified freely, although Rational
-- would be happy to receive comments, bug reports, and suggestions for
-- improvements at:
--                    Rational
--                    3320 Scott Boulevard
--                    Santa Clara, CA 95054 USA
--                 or by e-mail at:
--                    Support@Rational.COM
--

with Unchecked_Conversion;
package body Generic_Primitive_Functions is

-- This is a non-portable implementation for the R1000 architecture.
--
-- The R1000 supports only one floating-point type: FLOAT.
-- Values of the type FLOAT are represented using the double (ie. 64 bits)
-- float format specified by the IEEE Standard for Binary Floating Point
-- Arithmetic (ANSI/IEEE Std.754-1985).
-- The R1000 compiler does not support "not a number" values and the value
-- for negative zero.

   type Sign is new Boolean;
   for Sign'Size use 1;

   type Biased_Exponent is range 0 .. 2 ** 11 - 1;
   for Biased_Exponent'Size use 11;
   Exponent_Bias         : constant := -1022;
   Not_A_Number_Exp      : constant := 1025;
   Gradual_Underflow_Exp : constant := Exponent_Bias;

   type Mantissa is range 0 .. 2 ** 52 - 1;
   for Mantissa'Size use 52;

   type Float_Representation is
      record
         Sign_Bit      : Sign;
         Exponent_Part : Biased_Exponent;
         Mantissa_Part : Mantissa;
      end record;

   -- for Float_Representation use
   --     record
   --         Sign_Bit      at 0 range 0 .. 0;
   --         Exponent_Part at 0 range 1 .. 11;
   --         Mantissa_Part at 0 range 12 .. 63;
   --     end record;

   function Conv is new Unchecked_Conversion (Float_Type, Float_Representation);
   function Conv is new Unchecked_Conversion (Float_Representation, Float_Type);

   Largest_Representable_Positive : constant   := 2#1.0#E+53;
   Smallest_Positive              : Float_Type :=
      Conv (Float_Representation'
            (Sign_Bit => False, Exponent_Part => 0, Mantissa_Part => 1));
   Smallest_Negative              : Float_Type :=
      Conv (Float_Representation'
            (Sign_Bit => True, Exponent_Part => 0, Mantissa_Part => 1));


   type Complete_Mantissa is range 0 .. 2 ** 53 - 1;
   Mantissa_Complement : constant := 2 ** 52;

   type Mantissa_Bit_String is array (1 .. 52) of Boolean;
   pragma Pack (Mantissa_Bit_String);

   function Conv is new Unchecked_Conversion (Mantissa_Bit_String, Mantissa);
   function Conv is new Unchecked_Conversion (Mantissa, Mantissa_Bit_String);

-------------------------------------------------------------------------------

   function Exponent (X : Float_Type) return Exponent_Type is
      F : Float_Type;
      E : Exponent_Type;
   begin
      Decompose (X, Fraction => F, Exponent => E);
      return E;
   end Exponent;


   function Fraction (X : Float_Type) return Float_Type is
      F : Float_Type;
      E : Exponent_Type;
   begin
      Decompose (X, Fraction => F, Exponent => E);
      return F;
   end Fraction;


   procedure Decompose (X        : in  Float_Type;
                        Fraction : out Float_Type;
                        Exponent : out Exponent_Type) is

      X_Repr : Float_Representation := Conv (X);
      Exp    : Exponent_Type := Exponent_Type
                                   (X_Repr.Exponent_Part + Exponent_Bias);
      M      : Complete_Mantissa;
   begin
      if X = 0.0 then
         Fraction := 0.0;
         Exponent := 0;

      elsif Exp = Gradual_Underflow_Exp then
         -- renormalized number
         M := Complete_Mantissa (X_Repr.Mantissa_Part);
         while M < Mantissa_Complement loop
            M   := M * 2;
            Exp := Exp - 1;
         end loop;
         X_Repr.Mantissa_Part := Mantissa (M - Mantissa_Complement);
         X_Repr.Exponent_Part := Biased_Exponent (0 - Exponent_Bias);
         Fraction             := Conv (X_Repr);
         Exponent             := Exp;

      elsif Exp >= Not_A_Number_Exp then
         raise Constraint_Error;

      else -- normal case
         X_Repr.Exponent_Part := Biased_Exponent (0 - Exponent_Bias);
         Fraction             := Conv (X_Repr);
         Exponent             := Exp;
      end if;

   end Decompose;


   function Compose (Fraction : Float_Type; Exponent : Exponent_Type)
                    return Float_Type is
      X_Repr : Float_Representation := Conv (Fraction);
   begin
      if Fraction = 0.0 then
         return 0.0;

      elsif Exponent >= Not_A_Number_Exp then
         raise Constraint_Error;

      elsif Exponent <= Gradual_Underflow_Exp then
         return Copy_Sign (Float_Type'Small, Fraction); -- can be improved

      else
         X_Repr.Exponent_Part := Biased_Exponent (Exponent - Exponent_Bias);
         return Conv (X_Repr);
      end if;
   end Compose;

   function Scale (X : Float_Type; Exponent : Exponent_Type)
                  return Float_Type is
      X_Repr   : Float_Representation := Conv (X);
      Temp_Exp : Exponent_Type;
      X_Temp   : Float_Type;
   begin
      if X = 0.0 then
         return 0.0;
      end if;

      Temp_Exp := Exponent_Type (X_Repr.Exponent_Part) +
                     Exponent_Bias + Exponent;

      if Temp_Exp >= Not_A_Number_Exp then
         raise Constraint_Error;

      elsif Temp_Exp <= Gradual_Underflow_Exp then
         X_Repr.Exponent_Part := 1;
         X_Temp               := Conv (X_Repr);
         while Temp_Exp <= Gradual_Underflow_Exp loop
            X_Temp   := X_Temp * 0.5;
            Temp_Exp := Temp_Exp + 1;
         end loop;
         return X_Temp;
      else
         X_Repr.Exponent_Part := Biased_Exponent (Temp_Exp - Exponent_Bias);
         return Conv (X_Repr);
      end if;

   end Scale;


   function Floor (X : Float_Type) return Float_Type is
      Trunc : Float_Type;
   begin
      if abs X > Largest_Representable_Positive then
         return X;
      end if;
      Trunc := Truncate (X);
      if X >= 0.0 or else Trunc = X then
         return Trunc;
      else
         return Trunc - 1.0;
      end if;
   end Floor;


   function Ceiling (X : Float_Type) return Float_Type is
      Trunc : Float_Type;
   begin
      if abs X > Largest_Representable_Positive then
         return X;
      end if;
      Trunc := Truncate (X);
      if X <= 0.0 or else Trunc = X then
         return Trunc;
      else
         return Trunc + 1.0;
      end if;
   end Ceiling;

   function Is_Odd (F : Float_Type) return Boolean is
      F_Over_2 : Float_Type := F / 2.0;
   begin
      return (Truncate (F_Over_2) /= F_Over_2);
   end Is_Odd;

   function Round (X : Float_Type) return Float_Type is
      Abs_X       : Float_Type := abs X;
      Trunc, Frac : Float_Type;
   begin
      if Abs_X > Largest_Representable_Positive then
         return X;
      end if;
      Trunc := Truncate (Abs_X);
      Frac  := Abs_X - Trunc;
      if Frac > 0.5 then
         return Copy_Sign (Trunc + 1.0, X);
      elsif Frac < 0.5 then
         return Copy_Sign (Trunc, X);

         -- Frac = 0.5
         -- we have a tie: need to go to the even integer
      elsif Is_Odd (Trunc) then
         return Copy_Sign (Trunc + 1.0, X);
      else
         return Copy_Sign (Trunc, X);
      end if;

   end Round;


   function Truncate (X : Float_Type) return Float_Type is
   begin
      if abs X > Largest_Representable_Positive then
         return X;
      elsif abs X < 1.0 then
         return 0.0;
      else
         return Leading_Part (X, Integer (Exponent (X)));
      end if;
   end Truncate;


   function Remainder (X, Y : Float_Type) return Float_Type is

      -- Written by: P. T. P. Tang, Argonne National Laboratory
      -- Last revised: August 2, 1990
      A, B, Arg, P, Arg_Frac, P_Frac, Sign_X, Ieee_Rem : Float_Type;
      Arg_Exp, P_Exp : Exponent_Type;
      K : Integer;
      type Parity is (Even, Odd);
      P_Flag : Parity;

   begin

      if X > 0.0 then
         Sign_X := 1.0;
      else
         Sign_X := -1.0;
      end if;


      Arg := abs (X);
      P   := abs (Y);
      if Arg < P then
         P_Flag   := Even;
         Ieee_Rem := Arg;
      else
         Decompose (Arg, Arg_Frac, Arg_Exp);
         Decompose (P, P_Frac, P_Exp);
         P        := Compose (P_Frac, Arg_Exp);
         K        := Integer (Arg_Exp) - Integer (P_Exp);
         P_Flag   := Even;
         Ieee_Rem := Arg;
         for Cnt in reverse 0 .. K loop
            if Ieee_Rem >= P then
               P_Flag   := Odd;
               Ieee_Rem := Ieee_Rem - P;
            else
               P_Flag := Even;
            end if;
            P := P * 0.5;
         end loop;
      end if;
      --end of modulus remainder

      --final step to get IEEE remainder
      --here we need to compare Rem with abs(Y)/2; must be careful
      --of unrepresentable Y/2 possibly caused by subnormal numbers
      if P_Exp >= 0 then
         A := Ieee_Rem;
         B := abs (Y) * 0.5;
      else
         A := Ieee_Rem * 0.5;
         B := Y;
      end if;

      if ((A > B) or (A = B and P_Flag = Odd)) then
         Ieee_Rem := Ieee_Rem - abs (Y);
      end if;

      return (Sign_X * Ieee_Rem);

   end Remainder;


   function Adjacent (X, Towards : Float_Type) return Float_Type is
   begin
      if X = Towards then
         return X;
      elsif X < Towards then
         return Successor (X);
      else
         return Predecessor (X);
      end if;
   end Adjacent;


   function Successor (X : Float_Type) return Float_Type is
      F      : Float_Type;
      E      : Exponent_Type;
      X_Repr : Float_Representation := Conv (X);
      M      : Complete_Mantissa;
   begin
      if X >= Float_Type'Base'Last then
         raise Constraint_Error;
      end if;
      if X = 0.0 then
         return Smallest_Positive;
      end if;
      M := Complete_Mantissa (X_Repr.Mantissa_Part) + 1;
      if M > Complete_Mantissa (Mantissa'Last) then
         X_Repr.Mantissa_Part := 0;
         X_Repr.Exponent_Part := X_Repr.Exponent_Part + 1;
      else
         X_Repr.Mantissa_Part := Mantissa (M);
      end if;
      return Conv (X_Repr);
   end Successor;


   function Predecessor (X : Float_Type) return Float_Type is
      X_Repr : Float_Representation := Conv (X);
      M      : Mantissa;
   begin
      if X <= -Float_Type'Base'Last then
         raise Constraint_Error;
      end if;
      if X = 0.0 then
         return Smallest_Negative;
      end if;
      M := X_Repr.Mantissa_Part;
      if M = 0 then
         X_Repr.Mantissa_Part := Mantissa'Last;
         X_Repr.Exponent_Part := X_Repr.Exponent_Part - 1;
      else
         X_Repr.Mantissa_Part := M - 1;
      end if;
      return Conv (X_Repr);
   end Predecessor;


   function Copy_Sign (Value, Sign : Float_Type) return Float_Type is
   begin
      if Sign >= 0.0 then
         return abs Value;
      else
         return -(abs Value);
      end if;
   end Copy_Sign;


   function Leading_Part
               (X : Float_Type; Radix_Digits : Positive) return Float_Type is
      F      : Float_Type;
      E      : Exponent_Type;
      F_Repr : Float_Representation;
      M_Bits : Mantissa_Bit_String;
   begin
      Decompose (X, F, E);
      F_Repr := Conv (F);
      M_Bits := Conv (F_Repr.Mantissa_Part);
      for Bit in Radix_Digits .. Mantissa_Bit_String'Last loop
         M_Bits (Bit) := False;
      end loop;
      F_Repr.Mantissa_Part := Mantissa (Conv (M_Bits));
      F := Conv (F_Repr);
      return Compose (F, E);
   end Leading_Part;

end Generic_Primitive_Functions;