|
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 - metrics - downloadIndex: B T
Length: 24097 (0x5e21)
Types: TextFile
Names: »B«
└─⟦180fe333a⟧ Bits:30000405 8mm tape, Rational 1000, SW CATALOG, 10_20_0
└─⟦180fe333a⟧ Bits:30000537 8mm tape, Rational 1000, SW Catalog 10_20_0
└─⟦5cb1d1d7f⟧ »DATA«
└─⟦3b1ee7bd8⟧
└─⟦this⟧
with System;
with Unchecked_Conversion;
package body Universal_Integer_Arithmetic is
-- A universal integer consists of a sign and a magnitude. The
-- magnitude is a vector of non-negative integers giving from
-- most significant to least significant the "digits" of the
-- number in some convenient base. There are no leading zero digits,
-- unless the value is zero. Universal integers are always normalized.
-- The lower bound of the universal integer vector is always one.
-- Thus, the magnitude for the vector V(1 .. k) is given by :
--
-- V(1) * BASE**(k - 1) + V(2) * BASE**(k - 2) + ... + V(k)
--
-- The maximum number of digits in a universal integer is limited
-- in this implementation only by the amount of available memory.
--
-- For this implementation the BASE is 2**31. The universal digits are
-- integers in the range 0 .. BASE - 1. This choice of BASE means that
-- slightly less than half of the integer range is used. However, the
-- choice does ensure that the product of two universal digits is a
-- long integer.
-- Also, the number of universal digits required to represent an integer value
-- as a universal integer is at most one.
--
-- To complete the representation the high order universal digit has the sign
-- of the universal integer.
Base_B : constant := 31;
Base : constant Long_Integer := 2 ** Base_B;
Int_D : constant := 3;
subtype Universal_Digit is Long_Integer range -(Base - 1) .. (Base - 1);
type Vector is array (Positive range <>) of Universal_Digit;
-- I_Zero : constant Universal_Integer := new Vector'(1 => 0);
-- I_One : constant Universal_Integer := new Vector'(1 => 1);
-- I_Two : constant Universal_Integer := new Vector'(1 => 2);
-- I_Ten : constant Universal_Integer := new Vector'(1 => 10);
function I_Zero return Universal_Integer is
begin
return new Vector'(1 => 0);
end I_Zero;
function I_One return Universal_Integer is
begin
return new Vector'(1 => 1);
end I_One;
function I_Two return Universal_Integer is
begin
return new Vector'(1 => 2);
end I_Two;
function I_Ten return Universal_Integer is
begin
return new Vector'(1 => 10);
end I_Ten;
function Ui (V : Vector; S : Boolean := False) return Universal_Integer is
-- Constructs a universal integer from a vector and a sign; the vector
-- need not be normalized. The boolean s is true if the number is negative.
T : Universal_Integer;
begin
-- The representation used in this package requires that all
-- Universal_integer values be normalized. The first digit of any
-- value, except zero, must be non-zero.
for J in V'Range loop
if V (J) /= 0 then
T := new Vector (1 .. V'Last - J + 1);
-- ensure lower bound of one
T.all := V (J .. V'Last);
if S then
T (1) := -T (1);
end if;
return T;
end if;
end loop;
return I_Zero;
end Ui;
function Uv (V : Vector; S : Boolean := False) return Vector is
-- Constructs a universal integer from a vector and a sign; the vector
-- need not be normalized. The boolean s is true if the number is negative.
T : Universal_Integer;
R : Vector (1 .. V'Length);
begin
-- The representation used in this package requires that all
-- Universal_integer values be normalized. The first digit of any
-- value, except zero, must be non-zero.
for J in V'Range loop
if V (J) /= 0 then
-- ensure lower bound of one
R (1 .. R'Last - J + V'First) := V (J .. V'Last);
if S then
R (1) := -R (1);
end if;
return R (1 .. R'Last - J + V'First);
end if;
end loop;
return (1 => 0);
end Uv;
function Ui (I : Long_Integer) return Universal_Integer is
Y : Vector (1 .. Int_D) := (1 .. Int_D => 0);
Z : Long_Integer := Long_Integer (I);
begin
if I > -Base and then I < Base then
return new Vector'(1 => Z);
end if;
for J in reverse Y'Range loop
Y (J) := abs (Z rem Base);
Z := Z / Base;
end loop;
return Ui (Y, I < 0);
end Ui;
function Ui (I : Integer) return Universal_Integer is
begin
return Ui (Long_Integer (I));
end Ui;
function Int (X : Universal_Integer) return Long_Integer is
Y : Long_Integer;
begin
if X'Length = 1 then
return X (1);
end if;
Y := 0;
for I in X'Range loop
-- convert as a negative integer
Y := Y * Base - abs X (I);
-- this may raise NUMERIC_ERROR, but
-- only if the magnitude of x is too large.
end loop;
if X (1) < 0 then
return Y;
else
return -Y;
end if;
end Int;
function Int (X : Universal_Integer) return Integer is
begin
return Integer (Long_Integer'(Int (X)));
end Int;
function Image (X : Vector) return System.Byte_String is
subtype Int is Vector (1 .. X'Length);
Y : Int := X;
subtype Str is System.Byte_String
(1 .. X'Length * Universal_Digit'Size /
System.Byte'Size);
function Xbytes is new Unchecked_Conversion (Int, Str);
begin
return Str'(Xbytes (Int'(Y)));
end Image;
function Value (S : System.Byte_String) return Vector is
subtype Int is Vector (1 .. (S'Length * System.Byte'Size) /
Universal_Digit'Size);
subtype Str is System.Byte_String (1 .. S'Length);
function Xdigits is new Unchecked_Conversion (Str, Int);
begin
return Int'(Xdigits (Str (S)));
end Value;
function Value (S : System.Byte_String) return Universal_Integer is
begin
return new Vector'(Value (S));
end Value;
function Image (I : Universal_Integer) return System.Byte_String is
begin
return Image (I.all);
end Image;
function Image (X : Vector) return String is
Q : Vector (1 .. X'Length) := X; -- normalizes the vector at 1-origin
R : Long_Integer := 0;
T : Long_Integer;
begin
for N in Q'Range loop
T := R * Base + abs Q (N);
Q (N) := T / 10;
R := T rem 10;
end loop;
for N in Q'Range loop
if Q (N) /= 0 then
return Image (Q (N .. Q'Last)) &
Character'Val (Character'Pos ('0') + R);
end if;
end loop;
return String'(1 => Character'Val (Character'Pos ('0') + R));
end Image;
function Image (X : Universal_Integer) return String is
begin
if X (1) = 0 then
return "0";
elsif X (1) < 0 then
return '-' & Image (X.all);
else
return ' ' & Image (X.all);
end if;
end Image;
function Value (S : String) return Universal_Integer is
Num : Universal_Integer := I_Zero;
Exp : Integer := 0;
Signed : Boolean := False;
Has_Exp : Boolean := False;
C : Character;
J : Integer;
begin
if S'Length = 0 then
raise Constraint_Error;
end if;
J := S'First;
C := S (J);
if C = '-' or else C = '+' then
J := J + 1;
if S (J) not in '0' .. '9' then
-- index out of range may also raise
raise Constraint_Error; -- constraint_error here
end if;
Signed := C = '-';
end if;
while J <= S'Last loop
C := S (J);
case C is
when '0' .. '9' =>
if Has_Exp then
Exp := Exp * 10 + (Character'Pos (C) -
Character'Pos ('0'));
else
Num := Num * I_Ten + Ui (Integer (Character'Pos (C) -
Character'Pos ('0')));
end if;
when '_' =>
if S (J - 1) not in '0' .. '9' or else
S (J + 1) not in '0' .. '9' then
raise Constraint_Error;
end if;
when 'E' | 'e' =>
if Has_Exp or else S (J - 1) not in '0' .. '9' then
raise Constraint_Error;
end if;
Has_Exp := True;
if S (J + 1) = '+' then
J := J + 1;
end if;
if S (J + 1) not in '0' .. '9' then
raise Constraint_Error;
end if;
when others =>
raise Constraint_Error;
end case;
J := J + 1;
end loop;
if Has_Exp then
Num := Num * I_Ten ** Exp;
end if;
if Signed then
Num := -Num;
end if;
return Num;
end Value;
function "/" (X, Y : Vector) return Vector is
M : Integer;
X1, Y1 : Integer;
E : Long_Integer;
D, R, T : Long_Integer;
Qe : Long_Integer; -- quotient digit estimate
V1, V2 : Long_Integer;
begin
X1 := X'Length;
Y1 := Y'Length;
if X1 = 1 and then Y1 = 1 then
-- can use simple integer division
return (1 => X (1) / Y (1));
-- integer divide catches zero divisor
elsif Y1 = 1 then
-- divisor has single digit
-- dividend has more than one digit,
-- important special case for which
-- an efficient algorithm is used
R := 0;
V1 := abs Y (1);
if V1 = 0 then
-- divisor is zero
raise Numeric_Error;
end if;
declare
Q : Vector (1 .. X1);
begin
for J in X'Range loop
T := R * Base + abs X (J);
Q (J) := T / V1;
R := T rem V1;
end loop;
return Uv (Q, (X (1) < 0) xor (Y (1) < 0));
end;
end if;
-- At this point the length of the dividend is at least two and
-- at least as much as the length of the divisor. We must do a
-- full long division. The algorithm used here is from Knuth,
-- "The Art of Programming", Volume 2, Section 4.3.1, Algorithm D.
-- The first step is to multiply both the divisor and dividend
-- by a scale factor to ensure that the first digit of the divisor
-- is at least BASE / 2. This condition is required by the
-- quotient digit estimation algorithm used in the division loop.
-- Note that this may increase the size of the dividend by one digit
-- and thus the scaled dividend is placed in u.
M := X1 - Y1 + 1;
declare
U : Vector (1 .. X1 + 1); -- the dividend
V : Vector (1 .. Y1); -- the divisor
Q : Vector (1 .. M); -- the quotient
begin
U := 0 & abs X (1) & X (2 .. X1);
V := abs Y (1) & Y (2 .. Y1);
V1 := V (1);
D := Base / (V1 + 1); -- scale factor
if D > 1 then
-- scale dividend and divisor
R := 0;
for J in reverse U'Range loop
T := U (J) * D + R;
U (J) := T rem Base;
R := T / Base;
end loop;
R := 0;
for J in reverse V'Range loop
T := V (J) * D + R;
V (J) := T rem Base;
R := T / Base;
end loop;
end if;
-- This is the major loop, corresponding to long division steps.
V1 := V (1);
V2 := V (2);
for J in Q'Range loop
-- Guess the next quotient digit, qe, by dividing the first two
-- remaining dividend digits by the high order divisor digit.
-- This estimate is never low and is at most 2 high.
T := U (J) * Base + U (J + 1);
if U (J) /= V1 then
Qe := T / V1;
else
Qe := Base - 1;
end if;
-- Now refine this guess so that it is almost always correct and
-- is at worst one too high.
while V2 * Qe > (T - Qe * V1) * Base + U (J + 2) loop
Qe := Qe - 1;
end loop;
-- Using qe as the quotient digit, we multiply the divisor by
-- qe and subtract from the remaining dividend.
R := 0;
for K in reverse V'Range loop
T := U (J + K) - Qe * V (K) + R;
E := T rem Base;
R := T / Base;
if E < 0 then
E := E + Base;
R := R - 1;
end if;
U (J + K) := E;
end loop;
U (J) := U (J) + R;
-- If qe was off by one, then u(j) went negative when the last
-- carry was added. So we correct the error by subtracting one
-- from the quotient digit and adding back the divisor to the
-- relevant portion of the dividend.
if U (J) < 0 then
Qe := Qe - 1;
R := 0;
for K in reverse V'Range loop
T := U (J + K) + V (K) + R;
if T >= Base then
T := T - Base;
R := 1;
else
R := 0;
end if;
U (J + K) := T;
end loop;
U (J) := U (J) + R;
end if;
-- Store the next quotient digit.
Q (J) := Qe;
end loop;
return Uv (Q, (X (1) < 0) xor (Y (1) < 0));
end;
end "/";
procedure Scaled_Value (X : Universal_Integer;
Mantissa : out Vector;
Exponent : out Integer) is
-- normalize X to the length of Mantissa and set Exponent to the
-- correct scale factor to preserve the value of X
Q : Vector (Mantissa'First .. Mantissa'Last + 1);
Exp : Integer := (X'Length - Mantissa'Length) * Base_B;
B : Long_Integer := Base / 2;
Bb : Long_Integer := 1;
T : Long_Integer := 0;
begin
for I in Q'Range loop
if I > X'Last then
Q (I) := 0;
else
Q (I) := abs X (I);
end if;
end loop;
if Q (1) = 0 then
Exp := 0;
else
-- normalize the mantissa
while Q (1) < B loop
Exp := Exp - 1;
B := B / 2;
Bb := 2 * Bb;
end loop;
if Bb > 1 then
for I in reverse Q'Range loop
T := (T / Base) + Q (I) * Bb;
Q (I) := T mod Base;
end loop;
end if;
end if;
if X (1) < 0 then
Q (1) := -Q (1);
end if;
Mantissa := Q (Mantissa'Range);
Exponent := Exp;
end Scaled_Value;
function Scaled_Value (V : Vector; Exp : Integer) return Float is
R : Float := 0.0;
S : Float := 0.0;
B : Float := Float (Base);
X : Integer := Exp;
begin
for I in reverse V'Range loop
if S = 0.0 then
-- First time through or exponentiation has been
-- underflowing up to now. X keeps track of correct
-- exponent
S := 2.0 ** X;
X := X + Base_B;
else
S := S * B;
end if;
R := R + Float (abs V (I)) * S;
end loop;
if V (V'First) < 0 then
return -R;
else
return R;
end if;
end Scaled_Value;
function Scaled_Value
(X : Universal_Integer; Y : Universal_Integer) return Float is
Num : Vector (1 .. 4);
Den : Vector (1 .. 2);
Nx, Dx : Integer;
begin
Scaled_Value (X, Num, Nx);
Scaled_Value (Y, Den, Dx);
return Scaled_Value (Num / Den, Nx - Dx);
end Scaled_Value;
function "-" (X : Universal_Integer) return Universal_Integer is
begin
return new Vector'(-X (1) & X (2 .. X'Last));
end "-";
function "abs" (X : Universal_Integer) return Universal_Integer is
begin
return new Vector'(abs X (1) & X (2 .. X'Last));
end "abs";
function "+" (X, Y : Universal_Integer) return Universal_Integer is
M : Integer;
K, R : Long_Integer;
Xl, Yl : Integer;
Xs, Ys : Boolean;
begin
Xl := X'Length;
Yl := Y'Length;
if Xl = 1 and then Yl = 1 then
-- each one has a digit
return Ui (X (1) + Y (1));
else
-- either or both operands have > 1 digits
if Xl < Yl then
M := Yl + 1;
else
M := Xl + 1;
end if;
declare
U, V : Vector (1 .. M);
begin
Xs := X (1) < 0;
Ys := Y (1) < 0;
U := (1 .. M - Xl => 0) & abs X (1) & X (2 .. Xl);
V := (1 .. M - Yl => 0) & abs Y (1) & Y (2 .. Yl);
if Xs = Ys then
-- signs agree so add
K := 0;
for I in reverse 1 .. M loop
R := U (I) + V (I) + K;
if R >= Base then
R := R - Base;
K := 1;
else
K := 0;
end if;
U (I) := R;
end loop;
return Ui (U, Xs);
else
-- signs different, subtract smaller from larger
K := 0;
for I in reverse 1 .. M loop
R := U (I) - V (I) + K;
if R < 0 then
R := R + Base;
K := -1;
else
K := 0;
end if;
U (I) := R;
end loop;
if K = 0 then
-- x has the larger magnitude
return Ui (U, Xs);
else
-- y has the larger magnitude, so recomplement
K := 1;
for I in reverse 1 .. M loop
R := Base - 1 - U (I) + K;
if R = Base then
R := 0;
K := 1;
else
K := 0;
end if;
U (I) := R;
end loop;
return Ui (U, Ys);
end if;
end if;
end;
end if;
end "+";
function "-" (X, Y : Universal_Integer) return Universal_Integer is
begin
return X + (-Y);
end "-";
function "*" (X, Y : Universal_Integer) return Universal_Integer is
-- This function returns the product of the universal integers x
-- and y using essentially the familiar hand algorithm.
X1, Y1 : Integer;
begin
X1 := X'Length;
Y1 := Y'Length;
if X1 = 1 and Y1 = 1 then
-- both have a single digit
return Ui (X (1) * Y (1));
end if;
declare
W : Vector (1 .. X1 + Y1) := (1 .. X1 + Y1 => 0);
K, R : Long_Integer;
begin
for J in reverse Y'Range loop
-- outer loop through digits of the multiplier, inner loop
-- through digits of multiplicand
K := 0;
for I in reverse X'Range loop
R := abs (X (I) * Y (J)) + W (I + J) + K;
W (I + J) := R rem Base;
K := R / Base;
end loop;
W (J) := K;
end loop;
return Ui (W, (X (1) < 0) xor (Y (1) < 0));
end;
end "*";
function "/" (X, Y : Universal_Integer) return Universal_Integer is
begin
return new Vector'(X.all / Y.all);
end "/";
function "rem" (X, Y : Universal_Integer) return Universal_Integer is
begin
if X'Length = 1 and then Y'Length = 1 then
return Ui (X (1) rem Y (1));
else
return X - (X / Y) * Y;
end if;
end "rem";
function "mod" (X, Y : Universal_Integer) return Universal_Integer is
R : constant Universal_Integer := X rem Y;
begin
if (X (1) < 0) = (Y (1) < 0) or else R (1) = 0 then
return R;
else
return Y + R;
end if;
end "mod";
function "**" (X : Universal_Integer; Y : Integer)
return Universal_Integer is
-- Raise a universal integer to an integer power using the binary
-- representation of the exponent.
R : Universal_Integer := I_One;
V : Integer := Y;
T : Universal_Integer := abs X;
begin
if Y < 0 then
raise Constraint_Error;
elsif Y = 0 then
return I_One;
elsif X (1) = 0 then
return I_Zero;
elsif Y > 5280 then
raise Numeric_Error;
end if;
-- Starting the variable r at 1 and t at x loop through the binary
-- digits of v, squaring t each time, and multiplying the result r
-- by the current value of t each time a 1-bit is found.
while V /= 0 loop
-- v is odd
if V rem 2 = 1 then
R := R * T;
end if;
T := T * T; -- halve v
V := V / 2;
end loop;
-- compute the sign of the result: positive if y is even, the sign of
-- x if y is odd.
if X (1) < 0 and then Y rem 2 = 1 then
R (1) := -R (1);
end if;
return R;
end "**";
function ">=" (X, Y : Universal_Integer) return Boolean is
Z : Universal_Integer := X - Y;
begin
return Z (1) >= 0;
end ">=";
function "<=" (X, Y : Universal_Integer) return Boolean is
Z : Universal_Integer := X - Y;
begin
return Z (1) <= 0;
end "<=";
function "<" (X, Y : Universal_Integer) return Boolean is
Z : Universal_Integer := X - Y;
begin
return Z (1) < 0;
end "<";
function ">" (X, Y : Universal_Integer) return Boolean is
Z : Universal_Integer := X - Y;
begin
return X (1) > 0;
end ">";
function Eql (X, Y : Universal_Integer) return Boolean is
begin
return X.all = Y.all;
end Eql;
end Universal_Integer_Arithmetic;