--- /dev/null
+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;