|
DataMuseum.dkPresents historical artifacts from the history of: Rational R1000/400 Tapes |
This is an automatic "excavation" of a thematic subset of
See our Wiki for more about Rational R1000/400 Tapes Excavated with: AutoArchaeologist - Free & Open Source Software. |
top - downloadIndex: ┃ B T ┃
Length: 12287 (0x2fff)
Types: TextFile
Names: »B«
└─⟦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⟧
-- 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;