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;