-- --
-- B o d y --
-- --
--- Copyright (C) 1992-2005 Free Software Foundation, Inc. --
+-- Copyright (C) 1992-2006, Free Software Foundation, Inc. --
-- --
-- GNAT is free software; you can redistribute it and/or modify it under --
-- terms of the GNU General Public License as published by the Free Soft- --
-- value, since the issue is host representation of integer values.
Uint_Int_Last : Uint;
- -- Uint value containing Int'Last value set by Initialize.
+ -- Uint value containing Int'Last value set by Initialize
UI_Power_2 : array (Int range 0 .. 64) of Uint;
-- This table is used to memoize exponentiations by powers of 2. The Nth
-- digit of Vec contains the sign, all other digits are always non-
-- negative. Note that the input may be directly represented, and in
-- this case Vec will contain the corresponding one or two digit value.
+ -- The low bound of Vec is always 1.
function Least_Sig_Digit (Arg : Uint) return Int;
pragma Inline (Least_Sig_Digit);
function Sum_Double_Digits (Left : Uint; Sign : Int) return Int;
-- Same as above but work in New_Base = Base * Base
+ procedure UI_Div_Rem
+ (Left, Right : Uint;
+ Quotient : out Uint;
+ Remainder : out Uint;
+ Discard_Quotient : Boolean;
+ Discard_Remainder : Boolean);
+ -- Compute euclidian division of Left by Right, and return Quotient and
+ -- signed Remainder (Left rem Right).
+ -- If Discard_Quotient is True, Quotient is left unchanged.
+ -- If Discard_Remainder is True, Remainder is left unchanged.
+
function Vector_To_Uint
(In_Vec : UI_Vector;
- Negative : Boolean)
- return Uint;
+ Negative : Boolean) return Uint;
-- Functions that calculate values in UI_Vectors, call this function
-- to create and return the Uint value. In_Vec contains the multiple
-- precision (Base) representation of a non-negative value. Leading
procedure Init_Operand (UI : Uint; Vec : out UI_Vector) is
Loc : Int;
+ pragma Assert (Vec'First = Int'(1));
+
begin
if Direct (UI) then
Vec (1) := Direct_Val (UI);
Num : Nat;
begin
- if UI_Is_In_Int_Range (Input) then
+ -- Largest negative number has to be handled specially, since it is in
+ -- Int_Range, but we cannot take the absolute value.
+
+ if Input = Uint_Int_First then
+ return Int'Size;
+
+ -- For any other number in Int_Range, get absolute value of number
+
+ elsif UI_Is_In_Int_Range (Input) then
Num := abs (UI_To_Int (Input));
Bits := 0;
+ -- If not in Int_Range then initialize bit count for all low order
+ -- words, and set number to high order digit.
+
else
Bits := Base_Bits * (Uints.Table (Input).Length - 1);
Num := abs (Udigits.Table (Uints.Table (Input).Loc));
end if;
+ -- Increase bit count for remaining value in Num
+
while Types.">" (Num, 0) loop
Num := Num / 2;
Bits := Bits + 1;
-- Mathematically: assume base congruent to 1 and compute an equivelent
-- integer to Left.
- -- If Sign = -1 return the alternating sum of the "digits".
+ -- If Sign = -1 return the alternating sum of the "digits"
- -- D1 - D2 + D3 - D4 + D5 . . .
+ -- D1 - D2 + D3 - D4 + D5 ...
-- (where D1 is Least Significant Digit)
if Tmp_Int >= Base then
- -- Sign must be 1.
+ -- Sign must be 1
Tmp_Int := (Tmp_Int / Base) + 1;
end UI_Div;
function UI_Div (Left, Right : Uint) return Uint is
+ Quotient : Uint;
+ Remainder : Uint;
+ begin
+ UI_Div_Rem
+ (Left, Right,
+ Quotient, Remainder,
+ Discard_Quotient => False,
+ Discard_Remainder => True);
+ return Quotient;
+ end UI_Div;
+
+ ----------------
+ -- UI_Div_Rem --
+ ----------------
+
+ procedure UI_Div_Rem
+ (Left, Right : Uint;
+ Quotient : out Uint;
+ Remainder : out Uint;
+ Discard_Quotient : Boolean;
+ Discard_Remainder : Boolean)
+ is
begin
pragma Assert (Right /= Uint_0);
-- Cases where both operands are represented directly
if Direct (Left) and then Direct (Right) then
- return UI_From_Int (Direct_Val (Left) / Direct_Val (Right));
+ declare
+ DV_Left : constant Int := Direct_Val (Left);
+ DV_Right : constant Int := Direct_Val (Right);
+
+ begin
+ if not Discard_Quotient then
+ Quotient := UI_From_Int (DV_Left / DV_Right);
+ end if;
+
+ if not Discard_Remainder then
+ Remainder := UI_From_Int (DV_Left rem DV_Right);
+ end if;
+
+ return;
+ end;
end if;
declare
L_Vec : UI_Vector (1 .. L_Length);
R_Vec : UI_Vector (1 .. R_Length);
D : Int;
- Remainder : Int;
+ Remainder_I : Int;
Tmp_Divisor : Int;
Carry : Int;
Tmp_Int : Int;
Tmp_Dig : Int;
+ procedure UI_Div_Vector
+ (L_Vec : UI_Vector;
+ R_Int : Int;
+ Quotient : out UI_Vector;
+ Remainder : out Int);
+ pragma Inline (UI_Div_Vector);
+ -- Specialised variant for case where the divisor is a single digit
+
+ procedure UI_Div_Vector
+ (L_Vec : UI_Vector;
+ R_Int : Int;
+ Quotient : out UI_Vector;
+ Remainder : out Int)
+ is
+ Tmp_Int : Int;
+
+ begin
+ Remainder := 0;
+ for J in L_Vec'Range loop
+ Tmp_Int := Remainder * Base + abs L_Vec (J);
+ Quotient (Quotient'First + J - L_Vec'First) := Tmp_Int / R_Int;
+ Remainder := Tmp_Int rem R_Int;
+ end loop;
+
+ if L_Vec (L_Vec'First) < Int_0 then
+ Remainder := -Remainder;
+ end if;
+ end UI_Div_Vector;
+
+ -- Start of processing for UI_Div_Rem
+
begin
-- Result is zero if left operand is shorter than right
if L_Length < R_Length then
- return Uint_0;
+ if not Discard_Quotient then
+ Quotient := Uint_0;
+ end if;
+ if not Discard_Remainder then
+ Remainder := Left;
+ end if;
+ return;
end if;
Init_Operand (Left, L_Vec);
-- ordinary long division by hand).
if R_Length = Int_1 then
- Remainder := 0;
Tmp_Divisor := abs R_Vec (1);
declare
- Quotient : UI_Vector (1 .. L_Length);
+ Quotient_V : UI_Vector (1 .. L_Length);
begin
- for J in L_Vec'Range loop
- Tmp_Int := Remainder * Base + abs L_Vec (J);
- Quotient (J) := Tmp_Int / Tmp_Divisor;
- Remainder := Tmp_Int rem Tmp_Divisor;
- end loop;
+ UI_Div_Vector (L_Vec, Tmp_Divisor, Quotient_V, Remainder_I);
- return
- Vector_To_Uint
- (Quotient, (L_Vec (1) < Int_0 xor R_Vec (1) < Int_0));
+ if not Discard_Quotient then
+ Quotient :=
+ Vector_To_Uint
+ (Quotient_V, (L_Vec (1) < Int_0 xor R_Vec (1) < Int_0));
+ end if;
+
+ if not Discard_Remainder then
+ Remainder := UI_From_Int (Remainder_I);
+ end if;
+ return;
end;
end if;
Algorithm_D : declare
Dividend : UI_Vector (1 .. L_Length + 1);
Divisor : UI_Vector (1 .. R_Length);
- Quotient : UI_Vector (1 .. Q_Length);
+ Quotient_V : UI_Vector (1 .. Q_Length);
Divisor_Dig1 : Int;
Divisor_Dig2 : Int;
Q_Guess : Int;
Carry := Tmp_Int / Base;
end loop;
- -- Multiply Divisor by d.
+ -- Multiply Divisor by d
Carry := 0;
for J in reverse Divisor'Range loop
end loop;
end if;
- -- Main loop of long division algorithm.
+ -- Main loop of long division algorithm
Divisor_Dig1 := Divisor (1);
Divisor_Dig2 := Divisor (2);
- for J in Quotient'Range loop
+ for J in Quotient_V'Range loop
- -- [ CALCULATE Q (hat) ] (step D3 in the algorithm).
+ -- [ CALCULATE Q (hat) ] (step D3 in the algorithm)
Tmp_Int := Dividend (J) * Base + Dividend (J + 1);
Q_Guess := Q_Guess - 1;
end loop;
- -- [ MULTIPLY & SUBTRACT] (step D4). Q_Guess * Divisor is
+ -- [ MULTIPLY & SUBTRACT ] (step D4). Q_Guess * Divisor is
-- subtracted from the remaining dividend.
Carry := 0;
-- Finally we can get the next quotient digit
- Quotient (J) := Q_Guess;
+ Quotient_V (J) := Q_Guess;
end loop;
- return Vector_To_Uint
- (Quotient, (L_Vec (1) < Int_0 xor R_Vec (1) < Int_0));
+ -- [ UNNORMALIZE ] (step D8)
+
+ if not Discard_Quotient then
+ Quotient := Vector_To_Uint
+ (Quotient_V, (L_Vec (1) < Int_0 xor R_Vec (1) < Int_0));
+ end if;
+ if not Discard_Remainder then
+ declare
+ Remainder_V : UI_Vector (1 .. R_Length);
+ Discard_Int : Int;
+ begin
+ UI_Div_Vector
+ (Dividend (Dividend'Last - R_Length + 1 .. Dividend'Last),
+ D,
+ Remainder_V, Discard_Int);
+ Remainder := Vector_To_Uint (Remainder_V, L_Vec (1) < Int_0);
+ end;
+ end if;
end Algorithm_D;
end;
- end UI_Div;
+ end UI_Div_Rem;
------------
-- UI_Eq --
if Right = Uint_0 then
return Uint_1;
- -- 0 to any positive power is 0.
+ -- 0 to any positive power is 0
elsif Left = Uint_0 then
return Uint_0;
-- UI_GCD --
------------
- -- Lehmer's algorithm for GCD.
+ -- Lehmer's algorithm for GCD
-- The idea is to avoid using multiple precision arithmetic wherever
-- possible, substituting Int arithmetic instead. See Knuth volume II,
loop
-- We might overflow and get division by zero here. This just
- -- means we can not take the single precision step
+ -- means we cannot take the single precision step
Den1 := V_Hat + C;
Den2 := V_Hat + D;
if B = Int_0 then
- -- No single precision steps take a regular Euclid step.
+ -- No single precision steps take a regular Euclid step
Tmp_UI := U rem V;
U := V;
V := Tmp_UI;
else
- -- Use prior single precision steps to compute this Euclid step.
+ -- Use prior single precision steps to compute this Euclid step
-- Fixed bug 1415-008 spends 80% of its time working on this
-- step. Perhaps we need a special case Int / Uint dot
end if;
end UI_Mod;
+ -------------------------------
+ -- UI_Modular_Exponentiation --
+ -------------------------------
+
+ function UI_Modular_Exponentiation
+ (B : Uint;
+ E : Uint;
+ Modulo : Uint) return Uint
+ is
+ M : constant Save_Mark := Mark;
+
+ Result : Uint := Uint_1;
+ Base : Uint := B;
+ Exponent : Uint := E;
+
+ begin
+ while Exponent /= Uint_0 loop
+ if Least_Sig_Digit (Exponent) rem Int'(2) = Int'(1) then
+ Result := (Result * Base) rem Modulo;
+ end if;
+
+ Exponent := Exponent / Uint_2;
+ Base := (Base * Base) rem Modulo;
+ end loop;
+
+ Release_And_Save (M, Result);
+ return Result;
+ end UI_Modular_Exponentiation;
+
+ ------------------------
+ -- UI_Modular_Inverse --
+ ------------------------
+
+ function UI_Modular_Inverse (N : Uint; Modulo : Uint) return Uint is
+ M : constant Save_Mark := Mark;
+ U : Uint;
+ V : Uint;
+ Q : Uint;
+ R : Uint;
+ X : Uint;
+ Y : Uint;
+ T : Uint;
+ S : Int := 1;
+
+ begin
+ U := Modulo;
+ V := N;
+
+ X := Uint_1;
+ Y := Uint_0;
+
+ loop
+ UI_Div_Rem
+ (U, V,
+ Quotient => Q, Remainder => R,
+ Discard_Quotient => False,
+ Discard_Remainder => False);
+
+ U := V;
+ V := R;
+
+ T := X;
+ X := Y + Q * X;
+ Y := T;
+ S := -S;
+
+ exit when R = Uint_1;
+ end loop;
+
+ if S = Int'(-1) then
+ X := Modulo - X;
+ end if;
+
+ Release_And_Save (M, X);
+ return X;
+ end UI_Modular_Inverse;
+
------------
-- UI_Mul --
------------
return UI_From_Int (Direct_Val (Left) rem Direct_Val (Right));
else
+
-- Special cases when Right is less than 13 and Left is larger
-- larger than one digit. All of these algorithms depend on the
-- base being 2 ** 15 We work with Abs (Left) and Abs(Right)
-- and replace the rem with simpler operations where
-- possible.
- -- Least_Sig_Digit might return Negative numbers.
+ -- Least_Sig_Digit might return Negative numbers
when 2 =>
return UI_From_Int (
end if;
- -- Else fall through to general case.
+ -- Else fall through to general case
- -- ???This needs to be improved. We have the Rem when we do the
- -- Div. Div throws it away!
-
- -- The special case Length (Left) = Length(right) = 1 in Div
+ -- The special case Length (Left) = Length (Right) = 1 in Div
-- looks slow. It uses UI_To_Int when Int should suffice. ???
end if;
end if;
- return Left - (Left / Right) * Right;
+ declare
+ Quotient, Remainder : Uint;
+ begin
+ UI_Div_Rem
+ (Left, Right, Quotient, Remainder,
+ Discard_Quotient => True,
+ Discard_Remainder => False);
+ return Remainder;
+ end;
end UI_Rem;
------------
OSDN Git Service
