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