1with Interfaces; use Interfaces;
2
3with Ada.Unchecked_Conversion;
4
5package body Opt61_Pkg is
6
7   pragma Suppress (Overflow_Check);
8   pragma Suppress (Range_Check);
9
10   subtype Uns64 is Unsigned_64;
11
12   function To_Int is new Ada.Unchecked_Conversion (Uns64, Int64);
13
14   subtype Uns32 is Unsigned_32;
15
16   -----------------------
17   -- Local Subprograms --
18   -----------------------
19
20   function "+" (A : Uns64; B : Uns32) return Uns64 is (A + Uns64 (B));
21   --  Length doubling additions
22
23   function "*" (A, B : Uns32) return Uns64 is (Uns64 (A) * Uns64 (B));
24   --  Length doubling multiplication
25
26   function "&" (Hi, Lo : Uns32) return Uns64 is
27     (Shift_Left (Uns64 (Hi), 32) or Uns64 (Lo));
28   --  Concatenate hi, lo values to form 64-bit result
29
30   function "abs" (X : Int64) return Uns64 is
31     (if X = Int64'First then 2**63 else Uns64 (Int64'(abs X)));
32   --  Convert absolute value of X to unsigned. Note that we can't just use
33   --  the expression of the Else, because it overflows for X = Int64'First.
34
35   function Lo (A : Uns64) return Uns32 is (Uns32 (A and 16#FFFF_FFFF#));
36   --  Low order half of 64-bit value
37
38   function Hi (A : Uns64) return Uns32 is (Uns32 (Shift_Right (A, 32)));
39   --  High order half of 64-bit value
40
41   -------------------
42   -- Double_Divide --
43   -------------------
44
45   procedure Double_Divide
46     (X, Y, Z : Int64;
47      Q, R    : out Int64;
48      Round   : Boolean)
49   is
50      Xu  : constant Uns64 := abs X;
51      Yu  : constant Uns64 := abs Y;
52
53      Yhi : constant Uns32 := Hi (Yu);
54      Ylo : constant Uns32 := Lo (Yu);
55
56      Zu  : constant Uns64 := abs Z;
57      Zhi : constant Uns32 := Hi (Zu);
58      Zlo : constant Uns32 := Lo (Zu);
59
60      T1, T2     : Uns64;
61      Du, Qu, Ru : Uns64;
62      Den_Pos    : Boolean;
63
64   begin
65      if Yu = 0 or else Zu = 0 then
66         raise Constraint_Error;
67      end if;
68
69      --  Compute Y * Z. Note that if the result overflows 64 bits unsigned,
70      --  then the rounded result is clearly zero (since the dividend is at
71      --  most 2**63 - 1, the extra bit of precision is nice here).
72
73      if Yhi /= 0 then
74         if Zhi /= 0 then
75            Q := 0;
76            R := X;
77            return;
78         else
79            T2 := Yhi * Zlo;
80         end if;
81
82      else
83         T2 := (if Zhi /= 0 then Ylo * Zhi else 0);
84      end if;
85
86      T1 := Ylo * Zlo;
87      T2 := T2 + Hi (T1);
88
89      if Hi (T2) /= 0 then
90         Q := 0;
91         R := X;
92         return;
93      end if;
94
95      Du := Lo (T2) & Lo (T1);
96
97      --  Set final signs (RM 4.5.5(27-30))
98
99      Den_Pos := (Y < 0) = (Z < 0);
100
101      --  Check overflow case of largest negative number divided by 1
102
103      if X = Int64'First and then Du = 1 and then not Den_Pos then
104         raise Constraint_Error;
105      end if;
106
107      --  Perform the actual division
108
109      Qu := Xu / Du;
110      Ru := Xu rem Du;
111
112      --  Deal with rounding case
113
114      if Round and then Ru > (Du - Uns64'(1)) / Uns64'(2) then
115         Qu := Qu + Uns64'(1);
116      end if;
117
118      --  Case of dividend (X) sign positive
119
120      if X >= 0 then
121         R := To_Int (Ru);
122         Q := (if Den_Pos then To_Int (Qu) else -To_Int (Qu));
123
124      --  Case of dividend (X) sign negative
125
126      else
127         R := -To_Int (Ru);
128         Q := (if Den_Pos then -To_Int (Qu) else To_Int (Qu));
129      end if;
130   end Double_Divide;
131
132end Opt61_Pkg;
133