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. -- -- The base is 10 ** ((INTEGER'WIDTH - 2) / 2). 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 an integer. -- Also, the number of universal digits required to represent an integer value -- as a universal integer is at most four. -- -- To complete the representation the high order universal digit has the sign -- of the universal integer. Base_D : constant := (Integer'Width - 2) / 2; Base : constant := 10 ** Base_D; Base_Sq : constant := Base * Base; Int_D : constant := 4; type Vector is array (Positive range <>) of Integer; 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 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 Ui (I : Integer) return Universal_Integer is Y : Vector (1 .. Int_D) := (1 .. Int_D => 0); Z : Integer; begin if I < Base and then I > -Base then return new Vector'(1 => I); end if; Z := I; 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 Int (X : Universal_Integer) return Integer is Y : 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 end loop; -- only if the magnitude of x is too large. if X (1) < 0 then return Y; else return -Y; -- this may raise NUMERIC_ERROR if x is end if; -- -(integer'first) and range is not symmetric. end Int; function Image (X : Universal_Integer) return String is M : Integer := X'Length * Base_D + 1; S : String (1 .. M); Y : Universal_Integer; J, D : Integer; begin if X (1) = 0 then return " 0"; end if; J := M; Y := abs X; while Y (1) /= 0 loop D := Int (Y rem I_Ten); Y := Y / I_Ten; S (J) := Character'Val (Character'Pos ('0') + D); J := J - 1; end loop; if X (1) < 0 then S (J) := '-'; else S (J) := ' '; end if; D := M - J + 1; S (1 .. D) := S (J .. M); return S (1 .. D); 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 (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 : 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 : Integer; Xl, Yl : Integer; Xs, Ys : Boolean; begin Xl := X'Length; Yl := Y'Length; if Xl = 1 and then Yl = 1 then -- each has one 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. Xl, Yl : Integer; begin Xl := X'Length; Yl := Y'Length; if Xl = 1 and Yl = 1 then -- both have a single digit return Ui (X (1) * Y (1)); end if; declare W : Vector (1 .. Xl + Yl) := (1 .. Xl + Yl => 0); K, R : 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 M : Integer; Xl, Yl : Integer; E : Integer; D, R, T : Integer; Qe : Integer; -- quotient digit estimate V1, V2 : Integer; begin Xl := X'Length; Yl := Y'Length; if Xl = 1 and then Yl = 1 then -- can use simple integer division return Ui (X (1) / Y (1)); -- integer divide catches zero divisor elsif Xl < Yl then -- divisor has more digits return I_Zero; elsif Yl = 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 .. Xl); begin for J in X'Range loop T := R * Base + abs X (J); Q (J) := T / V1; R := T rem V1; end loop; return Ui (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 := Xl - Yl + 1; declare U : Vector (1 .. Xl + 1); -- the dividend V : Vector (1 .. Yl); -- the divisor Q : Vector (1 .. M); -- the quotient begin U := 0 & abs X (1) & X (2 .. Xl); V := abs Y (1) & Y (2 .. Yl); 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 Ui (Q, (X (1) < 0) xor (Y (1) < 0)); end; 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; 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 if V rem 2 = 1 then -- v is odd R := R * T; end if; T := T * T; V := V / 2; -- halve v 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 Z (1) > 0; end ">"; function Eql (X, Y : Universal_Integer) return Boolean is begin return X.all = Y.all; end Eql; end Universal_Integer_Arithmetic;