133 lines
3.3 KiB
Ada
133 lines
3.3 KiB
Ada
|
with Interfaces; use Interfaces;
|
||
|
|
||
|
with Ada.Unchecked_Conversion;
|
||
|
|
||
|
package body Opt61_Pkg is
|
||
|
|
||
|
pragma Suppress (Overflow_Check);
|
||
|
pragma Suppress (Range_Check);
|
||
|
|
||
|
subtype Uns64 is Unsigned_64;
|
||
|
|
||
|
function To_Int is new Ada.Unchecked_Conversion (Uns64, Int64);
|
||
|
|
||
|
subtype Uns32 is Unsigned_32;
|
||
|
|
||
|
-----------------------
|
||
|
-- Local Subprograms --
|
||
|
-----------------------
|
||
|
|
||
|
function "+" (A : Uns64; B : Uns32) return Uns64 is (A + Uns64 (B));
|
||
|
-- Length doubling additions
|
||
|
|
||
|
function "*" (A, B : Uns32) return Uns64 is (Uns64 (A) * Uns64 (B));
|
||
|
-- Length doubling multiplication
|
||
|
|
||
|
function "&" (Hi, Lo : Uns32) return Uns64 is
|
||
|
(Shift_Left (Uns64 (Hi), 32) or Uns64 (Lo));
|
||
|
-- Concatenate hi, lo values to form 64-bit result
|
||
|
|
||
|
function "abs" (X : Int64) return Uns64 is
|
||
|
(if X = Int64'First then 2**63 else Uns64 (Int64'(abs X)));
|
||
|
-- Convert absolute value of X to unsigned. Note that we can't just use
|
||
|
-- the expression of the Else, because it overflows for X = Int64'First.
|
||
|
|
||
|
function Lo (A : Uns64) return Uns32 is (Uns32 (A and 16#FFFF_FFFF#));
|
||
|
-- Low order half of 64-bit value
|
||
|
|
||
|
function Hi (A : Uns64) return Uns32 is (Uns32 (Shift_Right (A, 32)));
|
||
|
-- High order half of 64-bit value
|
||
|
|
||
|
-------------------
|
||
|
-- Double_Divide --
|
||
|
-------------------
|
||
|
|
||
|
procedure Double_Divide
|
||
|
(X, Y, Z : Int64;
|
||
|
Q, R : out Int64;
|
||
|
Round : Boolean)
|
||
|
is
|
||
|
Xu : constant Uns64 := abs X;
|
||
|
Yu : constant Uns64 := abs Y;
|
||
|
|
||
|
Yhi : constant Uns32 := Hi (Yu);
|
||
|
Ylo : constant Uns32 := Lo (Yu);
|
||
|
|
||
|
Zu : constant Uns64 := abs Z;
|
||
|
Zhi : constant Uns32 := Hi (Zu);
|
||
|
Zlo : constant Uns32 := Lo (Zu);
|
||
|
|
||
|
T1, T2 : Uns64;
|
||
|
Du, Qu, Ru : Uns64;
|
||
|
Den_Pos : Boolean;
|
||
|
|
||
|
begin
|
||
|
if Yu = 0 or else Zu = 0 then
|
||
|
raise Constraint_Error;
|
||
|
end if;
|
||
|
|
||
|
-- Compute Y * Z. Note that if the result overflows 64 bits unsigned,
|
||
|
-- then the rounded result is clearly zero (since the dividend is at
|
||
|
-- most 2**63 - 1, the extra bit of precision is nice here).
|
||
|
|
||
|
if Yhi /= 0 then
|
||
|
if Zhi /= 0 then
|
||
|
Q := 0;
|
||
|
R := X;
|
||
|
return;
|
||
|
else
|
||
|
T2 := Yhi * Zlo;
|
||
|
end if;
|
||
|
|
||
|
else
|
||
|
T2 := (if Zhi /= 0 then Ylo * Zhi else 0);
|
||
|
end if;
|
||
|
|
||
|
T1 := Ylo * Zlo;
|
||
|
T2 := T2 + Hi (T1);
|
||
|
|
||
|
if Hi (T2) /= 0 then
|
||
|
Q := 0;
|
||
|
R := X;
|
||
|
return;
|
||
|
end if;
|
||
|
|
||
|
Du := Lo (T2) & Lo (T1);
|
||
|
|
||
|
-- Set final signs (RM 4.5.5(27-30))
|
||
|
|
||
|
Den_Pos := (Y < 0) = (Z < 0);
|
||
|
|
||
|
-- Check overflow case of largest negative number divided by 1
|
||
|
|
||
|
if X = Int64'First and then Du = 1 and then not Den_Pos then
|
||
|
raise Constraint_Error;
|
||
|
end if;
|
||
|
|
||
|
-- Perform the actual division
|
||
|
|
||
|
Qu := Xu / Du;
|
||
|
Ru := Xu rem Du;
|
||
|
|
||
|
-- Deal with rounding case
|
||
|
|
||
|
if Round and then Ru > (Du - Uns64'(1)) / Uns64'(2) then
|
||
|
Qu := Qu + Uns64'(1);
|
||
|
end if;
|
||
|
|
||
|
-- Case of dividend (X) sign positive
|
||
|
|
||
|
if X >= 0 then
|
||
|
R := To_Int (Ru);
|
||
|
Q := (if Den_Pos then To_Int (Qu) else -To_Int (Qu));
|
||
|
|
||
|
-- Case of dividend (X) sign negative
|
||
|
|
||
|
else
|
||
|
R := -To_Int (Ru);
|
||
|
Q := (if Den_Pos then -To_Int (Qu) else To_Int (Qu));
|
||
|
end if;
|
||
|
end Double_Divide;
|
||
|
|
||
|
end Opt61_Pkg;
|