File : s-arit64.adb
1 ------------------------------------------------------------------------------
2 -- --
3 -- GNAT RUN-TIME COMPONENTS --
4 -- --
5 -- S Y S T E M . A R I T H _ 6 4 --
6 -- --
7 -- B o d y --
8 -- --
9 -- Copyright (C) 1992-2015, Free Software Foundation, Inc. --
10 -- --
11 -- GNAT is free software; you can redistribute it and/or modify it under --
12 -- terms of the GNU General Public License as published by the Free Soft- --
13 -- ware Foundation; either version 3, or (at your option) any later ver- --
14 -- sion. GNAT is distributed in the hope that it will be useful, but WITH- --
15 -- OUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY --
16 -- or FITNESS FOR A PARTICULAR PURPOSE. --
17 -- --
18 -- --
19 -- --
20 -- --
21 -- --
22 -- You should have received a copy of the GNU General Public License and --
23 -- a copy of the GCC Runtime Library Exception along with this program; --
24 -- see the files COPYING3 and COPYING.RUNTIME respectively. If not, see --
25 -- <http://www.gnu.org/licenses/>. --
26 -- --
27 -- GNAT was originally developed by the GNAT team at New York University. --
28 -- Extensive contributions were provided by Ada Core Technologies Inc. --
29 -- --
30 ------------------------------------------------------------------------------
31
32 with Interfaces; use Interfaces;
33
34 with Ada.Unchecked_Conversion;
35
36 package body System.Arith_64 is
37
38 pragma Suppress (Overflow_Check);
39 pragma Suppress (Range_Check);
40
41 subtype Uns64 is Unsigned_64;
42 function To_Uns is new Ada.Unchecked_Conversion (Int64, Uns64);
43 function To_Int is new Ada.Unchecked_Conversion (Uns64, Int64);
44
45 subtype Uns32 is Unsigned_32;
46
47 -----------------------
48 -- Local Subprograms --
49 -----------------------
50
51 function "+" (A, B : Uns32) return Uns64 is (Uns64 (A) + Uns64 (B));
52 function "+" (A : Uns64; B : Uns32) return Uns64 is (A + Uns64 (B));
53 -- Length doubling additions
54
55 function "*" (A, B : Uns32) return Uns64 is (Uns64 (A) * Uns64 (B));
56 -- Length doubling multiplication
57
58 function "/" (A : Uns64; B : Uns32) return Uns64 is (A / Uns64 (B));
59 -- Length doubling division
60
61 function "&" (Hi, Lo : Uns32) return Uns64 is
62 (Shift_Left (Uns64 (Hi), 32) or Uns64 (Lo));
63 -- Concatenate hi, lo values to form 64-bit result
64
65 function "abs" (X : Int64) return Uns64 is
66 (if X = Int64'First then 2**63 else Uns64 (Int64'(abs X)));
67 -- Convert absolute value of X to unsigned. Note that we can't just use
68 -- the expression of the Else, because it overflows for X = Int64'First.
69
70 function "rem" (A : Uns64; B : Uns32) return Uns64 is (A rem Uns64 (B));
71 -- Length doubling remainder
72
73 function Le3 (X1, X2, X3 : Uns32; Y1, Y2, Y3 : Uns32) return Boolean;
74 -- Determines if 96 bit value X1&X2&X3 <= Y1&Y2&Y3
75
76 function Lo (A : Uns64) return Uns32 is (Uns32 (A and 16#FFFF_FFFF#));
77 -- Low order half of 64-bit value
78
79 function Hi (A : Uns64) return Uns32 is (Uns32 (Shift_Right (A, 32)));
80 -- High order half of 64-bit value
81
82 procedure Sub3 (X1, X2, X3 : in out Uns32; Y1, Y2, Y3 : Uns32);
83 -- Computes X1&X2&X3 := X1&X2&X3 - Y1&Y1&Y3 with mod 2**96 wrap
84
85 function To_Neg_Int (A : Uns64) return Int64 with Inline;
86 -- Convert to negative integer equivalent. If the input is in the range
87 -- 0 .. 2 ** 63, then the corresponding negative signed integer (obtained
88 -- by negating the given value) is returned, otherwise constraint error
89 -- is raised.
90
91 function To_Pos_Int (A : Uns64) return Int64 with Inline;
92 -- Convert to positive integer equivalent. If the input is in the range
93 -- 0 .. 2 ** 63-1, then the corresponding non-negative signed integer is
94 -- returned, otherwise constraint error is raised.
95
96 procedure Raise_Error with Inline;
97 pragma No_Return (Raise_Error);
98 -- Raise constraint error with appropriate message
99
100 --------------------------
101 -- Add_With_Ovflo_Check --
102 --------------------------
103
104 function Add_With_Ovflo_Check (X, Y : Int64) return Int64 is
105 R : constant Int64 := To_Int (To_Uns (X) + To_Uns (Y));
106
107 begin
108 if X >= 0 then
109 if Y < 0 or else R >= 0 then
110 return R;
111 end if;
112
113 else -- X < 0
114 if Y > 0 or else R < 0 then
115 return R;
116 end if;
117 end if;
118
119 Raise_Error;
120 end Add_With_Ovflo_Check;
121
122 -------------------
123 -- Double_Divide --
124 -------------------
125
126 procedure Double_Divide
127 (X, Y, Z : Int64;
128 Q, R : out Int64;
129 Round : Boolean)
130 is
131 Xu : constant Uns64 := abs X;
132 Yu : constant Uns64 := abs Y;
133
134 Yhi : constant Uns32 := Hi (Yu);
135 Ylo : constant Uns32 := Lo (Yu);
136
137 Zu : constant Uns64 := abs Z;
138 Zhi : constant Uns32 := Hi (Zu);
139 Zlo : constant Uns32 := Lo (Zu);
140
141 T1, T2 : Uns64;
142 Du, Qu, Ru : Uns64;
143 Den_Pos : Boolean;
144
145 begin
146 if Yu = 0 or else Zu = 0 then
147 Raise_Error;
148 end if;
149
150 -- Compute Y * Z. Note that if the result overflows 64 bits unsigned,
151 -- then the rounded result is clearly zero (since the dividend is at
152 -- most 2**63 - 1, the extra bit of precision is nice here).
153
154 if Yhi /= 0 then
155 if Zhi /= 0 then
156 Q := 0;
157 R := X;
158 return;
159 else
160 T2 := Yhi * Zlo;
161 end if;
162
163 else
164 T2 := (if Zhi /= 0 then Ylo * Zhi else 0);
165 end if;
166
167 T1 := Ylo * Zlo;
168 T2 := T2 + Hi (T1);
169
170 if Hi (T2) /= 0 then
171 Q := 0;
172 R := X;
173 return;
174 end if;
175
176 Du := Lo (T2) & Lo (T1);
177
178 -- Set final signs (RM 4.5.5(27-30))
179
180 Den_Pos := (Y < 0) = (Z < 0);
181
182 -- Check overflow case of largest negative number divided by 1
183
184 if X = Int64'First and then Du = 1 and then not Den_Pos then
185 Raise_Error;
186 end if;
187
188 -- Perform the actual division
189
190 Qu := Xu / Du;
191 Ru := Xu rem Du;
192
193 -- Deal with rounding case
194
195 if Round and then Ru > (Du - Uns64'(1)) / Uns64'(2) then
196 Qu := Qu + Uns64'(1);
197 end if;
198
199 -- Case of dividend (X) sign positive
200
201 if X >= 0 then
202 R := To_Int (Ru);
203 Q := (if Den_Pos then To_Int (Qu) else -To_Int (Qu));
204
205 -- Case of dividend (X) sign negative
206
207 else
208 R := -To_Int (Ru);
209 Q := (if Den_Pos then -To_Int (Qu) else To_Int (Qu));
210 end if;
211 end Double_Divide;
212
213 ---------
214 -- Le3 --
215 ---------
216
217 function Le3 (X1, X2, X3 : Uns32; Y1, Y2, Y3 : Uns32) return Boolean is
218 begin
219 if X1 < Y1 then
220 return True;
221 elsif X1 > Y1 then
222 return False;
223 elsif X2 < Y2 then
224 return True;
225 elsif X2 > Y2 then
226 return False;
227 else
228 return X3 <= Y3;
229 end if;
230 end Le3;
231
232 -------------------------------
233 -- Multiply_With_Ovflo_Check --
234 -------------------------------
235
236 function Multiply_With_Ovflo_Check (X, Y : Int64) return Int64 is
237 Xu : constant Uns64 := abs X;
238 Xhi : constant Uns32 := Hi (Xu);
239 Xlo : constant Uns32 := Lo (Xu);
240
241 Yu : constant Uns64 := abs Y;
242 Yhi : constant Uns32 := Hi (Yu);
243 Ylo : constant Uns32 := Lo (Yu);
244
245 T1, T2 : Uns64;
246
247 begin
248 if Xhi /= 0 then
249 if Yhi /= 0 then
250 Raise_Error;
251 else
252 T2 := Xhi * Ylo;
253 end if;
254
255 elsif Yhi /= 0 then
256 T2 := Xlo * Yhi;
257
258 else -- Yhi = Xhi = 0
259 T2 := 0;
260 end if;
261
262 -- Here we have T2 set to the contribution to the upper half of the
263 -- result from the upper halves of the input values.
264
265 T1 := Xlo * Ylo;
266 T2 := T2 + Hi (T1);
267
268 if Hi (T2) /= 0 then
269 Raise_Error;
270 end if;
271
272 T2 := Lo (T2) & Lo (T1);
273
274 if X >= 0 then
275 if Y >= 0 then
276 return To_Pos_Int (T2);
277 else
278 return To_Neg_Int (T2);
279 end if;
280 else -- X < 0
281 if Y < 0 then
282 return To_Pos_Int (T2);
283 else
284 return To_Neg_Int (T2);
285 end if;
286 end if;
287
288 end Multiply_With_Ovflo_Check;
289
290 -----------------
291 -- Raise_Error --
292 -----------------
293
294 procedure Raise_Error is
295 begin
296 raise Constraint_Error with "64-bit arithmetic overflow";
297 end Raise_Error;
298
299 -------------------
300 -- Scaled_Divide --
301 -------------------
302
303 procedure Scaled_Divide
304 (X, Y, Z : Int64;
305 Q, R : out Int64;
306 Round : Boolean)
307 is
308 Xu : constant Uns64 := abs X;
309 Xhi : constant Uns32 := Hi (Xu);
310 Xlo : constant Uns32 := Lo (Xu);
311
312 Yu : constant Uns64 := abs Y;
313 Yhi : constant Uns32 := Hi (Yu);
314 Ylo : constant Uns32 := Lo (Yu);
315
316 Zu : Uns64 := abs Z;
317 Zhi : Uns32 := Hi (Zu);
318 Zlo : Uns32 := Lo (Zu);
319
320 D : array (1 .. 4) of Uns32;
321 -- The dividend, four digits (D(1) is high order)
322
323 Qd : array (1 .. 2) of Uns32;
324 -- The quotient digits, two digits (Qd(1) is high order)
325
326 S1, S2, S3 : Uns32;
327 -- Value to subtract, three digits (S1 is high order)
328
329 Qu : Uns64;
330 Ru : Uns64;
331 -- Unsigned quotient and remainder
332
333 Scale : Natural;
334 -- Scaling factor used for multiple-precision divide. Dividend and
335 -- Divisor are multiplied by 2 ** Scale, and the final remainder is
336 -- divided by the scaling factor. The reason for this scaling is to
337 -- allow more accurate estimation of quotient digits.
338
339 T1, T2, T3 : Uns64;
340 -- Temporary values
341
342 begin
343 -- First do the multiplication, giving the four digit dividend
344
345 T1 := Xlo * Ylo;
346 D (4) := Lo (T1);
347 D (3) := Hi (T1);
348
349 if Yhi /= 0 then
350 T1 := Xlo * Yhi;
351 T2 := D (3) + Lo (T1);
352 D (3) := Lo (T2);
353 D (2) := Hi (T1) + Hi (T2);
354
355 if Xhi /= 0 then
356 T1 := Xhi * Ylo;
357 T2 := D (3) + Lo (T1);
358 D (3) := Lo (T2);
359 T3 := D (2) + Hi (T1);
360 T3 := T3 + Hi (T2);
361 D (2) := Lo (T3);
362 D (1) := Hi (T3);
363
364 T1 := (D (1) & D (2)) + Uns64'(Xhi * Yhi);
365 D (1) := Hi (T1);
366 D (2) := Lo (T1);
367
368 else
369 D (1) := 0;
370 end if;
371
372 else
373 if Xhi /= 0 then
374 T1 := Xhi * Ylo;
375 T2 := D (3) + Lo (T1);
376 D (3) := Lo (T2);
377 D (2) := Hi (T1) + Hi (T2);
378
379 else
380 D (2) := 0;
381 end if;
382
383 D (1) := 0;
384 end if;
385
386 -- Now it is time for the dreaded multiple precision division. First an
387 -- easy case, check for the simple case of a one digit divisor.
388
389 if Zhi = 0 then
390 if D (1) /= 0 or else D (2) >= Zlo then
391 Raise_Error;
392
393 -- Here we are dividing at most three digits by one digit
394
395 else
396 T1 := D (2) & D (3);
397 T2 := Lo (T1 rem Zlo) & D (4);
398
399 Qu := Lo (T1 / Zlo) & Lo (T2 / Zlo);
400 Ru := T2 rem Zlo;
401 end if;
402
403 -- If divisor is double digit and too large, raise error
404
405 elsif (D (1) & D (2)) >= Zu then
406 Raise_Error;
407
408 -- This is the complex case where we definitely have a double digit
409 -- divisor and a dividend of at least three digits. We use the classical
410 -- multiple division algorithm (see section (4.3.1) of Knuth's "The Art
411 -- of Computer Programming", Vol. 2 for a description (algorithm D).
412
413 else
414 -- First normalize the divisor so that it has the leading bit on.
415 -- We do this by finding the appropriate left shift amount.
416
417 Scale := 0;
418
419 if (Zhi and 16#FFFF0000#) = 0 then
420 Scale := 16;
421 Zu := Shift_Left (Zu, 16);
422 end if;
423
424 if (Hi (Zu) and 16#FF00_0000#) = 0 then
425 Scale := Scale + 8;
426 Zu := Shift_Left (Zu, 8);
427 end if;
428
429 if (Hi (Zu) and 16#F000_0000#) = 0 then
430 Scale := Scale + 4;
431 Zu := Shift_Left (Zu, 4);
432 end if;
433
434 if (Hi (Zu) and 16#C000_0000#) = 0 then
435 Scale := Scale + 2;
436 Zu := Shift_Left (Zu, 2);
437 end if;
438
439 if (Hi (Zu) and 16#8000_0000#) = 0 then
440 Scale := Scale + 1;
441 Zu := Shift_Left (Zu, 1);
442 end if;
443
444 Zhi := Hi (Zu);
445 Zlo := Lo (Zu);
446
447 -- Note that when we scale up the dividend, it still fits in four
448 -- digits, since we already tested for overflow, and scaling does
449 -- not change the invariant that (D (1) & D (2)) >= Zu.
450
451 T1 := Shift_Left (D (1) & D (2), Scale);
452 D (1) := Hi (T1);
453 T2 := Shift_Left (0 & D (3), Scale);
454 D (2) := Lo (T1) or Hi (T2);
455 T3 := Shift_Left (0 & D (4), Scale);
456 D (3) := Lo (T2) or Hi (T3);
457 D (4) := Lo (T3);
458
459 -- Loop to compute quotient digits, runs twice for Qd(1) and Qd(2)
460
461 for J in 0 .. 1 loop
462
463 -- Compute next quotient digit. We have to divide three digits by
464 -- two digits. We estimate the quotient by dividing the leading
465 -- two digits by the leading digit. Given the scaling we did above
466 -- which ensured the first bit of the divisor is set, this gives
467 -- an estimate of the quotient that is at most two too high.
468
469 Qd (J + 1) := (if D (J + 1) = Zhi
470 then 2 ** 32 - 1
471 else Lo ((D (J + 1) & D (J + 2)) / Zhi));
472
473 -- Compute amount to subtract
474
475 T1 := Qd (J + 1) * Zlo;
476 T2 := Qd (J + 1) * Zhi;
477 S3 := Lo (T1);
478 T1 := Hi (T1) + Lo (T2);
479 S2 := Lo (T1);
480 S1 := Hi (T1) + Hi (T2);
481
482 -- Adjust quotient digit if it was too high
483
484 loop
485 exit when Le3 (S1, S2, S3, D (J + 1), D (J + 2), D (J + 3));
486 Qd (J + 1) := Qd (J + 1) - 1;
487 Sub3 (S1, S2, S3, 0, Zhi, Zlo);
488 end loop;
489
490 -- Now subtract S1&S2&S3 from D1&D2&D3 ready for next step
491
492 Sub3 (D (J + 1), D (J + 2), D (J + 3), S1, S2, S3);
493 end loop;
494
495 -- The two quotient digits are now set, and the remainder of the
496 -- scaled division is in D3&D4. To get the remainder for the
497 -- original unscaled division, we rescale this dividend.
498
499 -- We rescale the divisor as well, to make the proper comparison
500 -- for rounding below.
501
502 Qu := Qd (1) & Qd (2);
503 Ru := Shift_Right (D (3) & D (4), Scale);
504 Zu := Shift_Right (Zu, Scale);
505 end if;
506
507 -- Deal with rounding case
508
509 if Round and then Ru > (Zu - Uns64'(1)) / Uns64'(2) then
510 Qu := Qu + Uns64 (1);
511 end if;
512
513 -- Set final signs (RM 4.5.5(27-30))
514
515 -- Case of dividend (X * Y) sign positive
516
517 if (X >= 0 and then Y >= 0) or else (X < 0 and then Y < 0) then
518 R := To_Pos_Int (Ru);
519 Q := (if Z > 0 then To_Pos_Int (Qu) else To_Neg_Int (Qu));
520
521 -- Case of dividend (X * Y) sign negative
522
523 else
524 R := To_Neg_Int (Ru);
525 Q := (if Z > 0 then To_Neg_Int (Qu) else To_Pos_Int (Qu));
526 end if;
527 end Scaled_Divide;
528
529 ----------
530 -- Sub3 --
531 ----------
532
533 procedure Sub3 (X1, X2, X3 : in out Uns32; Y1, Y2, Y3 : Uns32) is
534 begin
535 if Y3 > X3 then
536 if X2 = 0 then
537 X1 := X1 - 1;
538 end if;
539
540 X2 := X2 - 1;
541 end if;
542
543 X3 := X3 - Y3;
544
545 if Y2 > X2 then
546 X1 := X1 - 1;
547 end if;
548
549 X2 := X2 - Y2;
550 X1 := X1 - Y1;
551 end Sub3;
552
553 -------------------------------
554 -- Subtract_With_Ovflo_Check --
555 -------------------------------
556
557 function Subtract_With_Ovflo_Check (X, Y : Int64) return Int64 is
558 R : constant Int64 := To_Int (To_Uns (X) - To_Uns (Y));
559
560 begin
561 if X >= 0 then
562 if Y > 0 or else R >= 0 then
563 return R;
564 end if;
565
566 else -- X < 0
567 if Y <= 0 or else R < 0 then
568 return R;
569 end if;
570 end if;
571
572 Raise_Error;
573 end Subtract_With_Ovflo_Check;
574
575 ----------------
576 -- To_Neg_Int --
577 ----------------
578
579 function To_Neg_Int (A : Uns64) return Int64 is
580 R : constant Int64 := (if A = 2**63 then Int64'First else -To_Int (A));
581 -- Note that we can't just use the expression of the Else, because it
582 -- overflows for A = 2**63.
583 begin
584 if R <= 0 then
585 return R;
586 else
587 Raise_Error;
588 end if;
589 end To_Neg_Int;
590
591 ----------------
592 -- To_Pos_Int --
593 ----------------
594
595 function To_Pos_Int (A : Uns64) return Int64 is
596 R : constant Int64 := To_Int (A);
597 begin
598 if R >= 0 then
599 return R;
600 else
601 Raise_Error;
602 end if;
603 end To_Pos_Int;
604
605 end System.Arith_64;