1 ------------------------------------------------------------------------------
3 -- GNAT RUN-TIME COMPONENTS --
5 -- S Y S T E M . A R I T H _ 6 4 --
9 -- Copyright (C) 1992-2004 Free Software Foundation, Inc. --
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 2, 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. See the GNU General Public License --
17 -- for more details. You should have received a copy of the GNU General --
18 -- Public License distributed with GNAT; see file COPYING. If not, write --
19 -- to the Free Software Foundation, 51 Franklin Street, Fifth Floor, --
20 -- Boston, MA 02110-1301, USA. --
22 -- As a special exception, if other files instantiate generics from this --
23 -- unit, or you link this unit with other files to produce an executable, --
24 -- this unit does not by itself cause the resulting executable to be --
25 -- covered by the GNU General Public License. This exception does not --
26 -- however invalidate any other reasons why the executable file might be --
27 -- covered by the GNU Public License. --
29 -- GNAT was originally developed by the GNAT team at New York University. --
30 -- Extensive contributions were provided by Ada Core Technologies Inc. --
32 ------------------------------------------------------------------------------
34 with System.Pure_Exceptions; use System.Pure_Exceptions;
36 with Interfaces; use Interfaces;
37 with Unchecked_Conversion;
39 package body System.Arith_64 is
41 pragma Suppress (Overflow_Check);
42 pragma Suppress (Range_Check);
44 subtype Uns64 is Unsigned_64;
45 function To_Uns is new Unchecked_Conversion (Int64, Uns64);
46 function To_Int is new Unchecked_Conversion (Uns64, Int64);
48 subtype Uns32 is Unsigned_32;
50 -----------------------
51 -- Local Subprograms --
52 -----------------------
54 function "+" (A, B : Uns32) return Uns64;
55 function "+" (A : Uns64; B : Uns32) return Uns64;
57 -- Length doubling additions
59 function "*" (A, B : Uns32) return Uns64;
61 -- Length doubling multiplication
63 function "/" (A : Uns64; B : Uns32) return Uns64;
65 -- Length doubling division
67 function "rem" (A : Uns64; B : Uns32) return Uns64;
68 pragma Inline ("rem");
69 -- Length doubling remainder
71 function "&" (Hi, Lo : Uns32) return Uns64;
73 -- Concatenate hi, lo values to form 64-bit result
75 function Le3 (X1, X2, X3 : Uns32; Y1, Y2, Y3 : Uns32) return Boolean;
76 -- Determines if 96 bit value X1&X2&X3 <= Y1&Y2&Y3
78 function Lo (A : Uns64) return Uns32;
80 -- Low order half of 64-bit value
82 function Hi (A : Uns64) return Uns32;
84 -- High order half of 64-bit value
86 procedure Sub3 (X1, X2, X3 : in out Uns32; Y1, Y2, Y3 : in Uns32);
87 -- Computes X1&X2&X3 := X1&X2&X3 - Y1&Y1&Y3 with mod 2**96 wrap
89 function To_Neg_Int (A : Uns64) return Int64;
90 -- Convert to negative integer equivalent. If the input is in the range
91 -- 0 .. 2 ** 63, then the corresponding negative signed integer (obtained
92 -- by negating the given value) is returned, otherwise constraint error
95 function To_Pos_Int (A : Uns64) return Int64;
96 -- Convert to positive integer equivalent. If the input is in the range
97 -- 0 .. 2 ** 63-1, then the corresponding non-negative signed integer is
98 -- returned, otherwise constraint error is raised.
100 procedure Raise_Error;
101 pragma No_Return (Raise_Error);
102 -- Raise constraint error with appropriate message
108 function "&" (Hi, Lo : Uns32) return Uns64 is
110 return Shift_Left (Uns64 (Hi), 32) or Uns64 (Lo);
117 function "*" (A, B : Uns32) return Uns64 is
119 return Uns64 (A) * Uns64 (B);
126 function "+" (A, B : Uns32) return Uns64 is
128 return Uns64 (A) + Uns64 (B);
131 function "+" (A : Uns64; B : Uns32) return Uns64 is
133 return A + Uns64 (B);
140 function "/" (A : Uns64; B : Uns32) return Uns64 is
142 return A / Uns64 (B);
149 function "rem" (A : Uns64; B : Uns32) return Uns64 is
151 return A rem Uns64 (B);
154 --------------------------
155 -- Add_With_Ovflo_Check --
156 --------------------------
158 function Add_With_Ovflo_Check (X, Y : Int64) return Int64 is
159 R : constant Int64 := To_Int (To_Uns (X) + To_Uns (Y));
163 if Y < 0 or else R >= 0 then
168 if Y > 0 or else R < 0 then
174 end Add_With_Ovflo_Check;
180 procedure Double_Divide
185 Xu : constant Uns64 := To_Uns (abs X);
186 Yu : constant Uns64 := To_Uns (abs Y);
188 Yhi : constant Uns32 := Hi (Yu);
189 Ylo : constant Uns32 := Lo (Yu);
191 Zu : constant Uns64 := To_Uns (abs Z);
192 Zhi : constant Uns32 := Hi (Zu);
193 Zlo : constant Uns32 := Lo (Zu);
200 if Yu = 0 or else Zu = 0 then
204 -- Compute Y * Z. Note that if the result overflows 64 bits unsigned,
205 -- then the rounded result is clearly zero (since the dividend is at
206 -- most 2**63 - 1, the extra bit of precision is nice here!)
234 Du := Lo (T2) & Lo (T1);
236 -- Set final signs (RM 4.5.5(27-30))
238 Den_Pos := (Y < 0) = (Z < 0);
240 -- Check overflow case of largest negative number divided by 1
242 if X = Int64'First and then Du = 1 and then not Den_Pos then
246 -- Perform the actual division
251 -- Deal with rounding case
253 if Round and then Ru > (Du - Uns64'(1)) / Uns64'(2) then
254 Qu := Qu + Uns64'(1);
257 -- Case of dividend (X) sign positive
268 -- Case of dividend (X) sign negative
285 function Hi (A : Uns64) return Uns32 is
287 return Uns32 (Shift_Right (A, 32));
294 function Le3 (X1, X2, X3 : Uns32; Y1, Y2, Y3 : Uns32) return Boolean is
313 function Lo (A : Uns64) return Uns32 is
315 return Uns32 (A and 16#FFFF_FFFF#);
318 -------------------------------
319 -- Multiply_With_Ovflo_Check --
320 -------------------------------
322 function Multiply_With_Ovflo_Check (X, Y : Int64) return Int64 is
323 Xu : constant Uns64 := To_Uns (abs X);
324 Xhi : constant Uns32 := Hi (Xu);
325 Xlo : constant Uns32 := Lo (Xu);
327 Yu : constant Uns64 := To_Uns (abs Y);
328 Yhi : constant Uns32 := Hi (Yu);
329 Ylo : constant Uns32 := Lo (Yu);
344 else -- Yhi = Xhi = 0
348 -- Here we have T2 set to the contribution to the upper half
349 -- of the result from the upper halves of the input values.
358 T2 := Lo (T2) & Lo (T1);
362 return To_Pos_Int (T2);
364 return To_Neg_Int (T2);
368 return To_Pos_Int (T2);
370 return To_Neg_Int (T2);
374 end Multiply_With_Ovflo_Check;
380 procedure Raise_Error is
382 Raise_Exception (CE, "64-bit arithmetic overflow");
389 procedure Scaled_Divide
394 Xu : constant Uns64 := To_Uns (abs X);
395 Xhi : constant Uns32 := Hi (Xu);
396 Xlo : constant Uns32 := Lo (Xu);
398 Yu : constant Uns64 := To_Uns (abs Y);
399 Yhi : constant Uns32 := Hi (Yu);
400 Ylo : constant Uns32 := Lo (Yu);
402 Zu : Uns64 := To_Uns (abs Z);
403 Zhi : Uns32 := Hi (Zu);
404 Zlo : Uns32 := Lo (Zu);
406 D : array (1 .. 4) of Uns32;
407 -- The dividend, four digits (D(1) is high order)
409 Qd : array (1 .. 2) of Uns32;
410 -- The quotient digits, two digits (Qd(1) is high order)
413 -- Value to subtract, three digits (S1 is high order)
417 -- Unsigned quotient and remainder
420 -- Scaling factor used for multiple-precision divide. Dividend and
421 -- Divisor are multiplied by 2 ** Scale, and the final remainder
422 -- is divided by the scaling factor. The reason for this scaling
423 -- is to allow more accurate estimation of quotient digits.
429 -- First do the multiplication, giving the four digit dividend
437 T2 := D (3) + Lo (T1);
439 D (2) := Hi (T1) + Hi (T2);
443 T2 := D (3) + Lo (T1);
445 T3 := D (2) + Hi (T1);
450 T1 := (D (1) & D (2)) + Uns64'(Xhi * Yhi);
461 T2 := D (3) + Lo (T1);
463 D (2) := Hi (T1) + Hi (T2);
472 -- Now it is time for the dreaded multiple precision division. First
473 -- an easy case, check for the simple case of a one digit divisor.
476 if D (1) /= 0 or else D (2) >= Zlo then
479 -- Here we are dividing at most three digits by one digit
483 T2 := Lo (T1 rem Zlo) & D (4);
485 Qu := Lo (T1 / Zlo) & Lo (T2 / Zlo);
489 -- If divisor is double digit and too large, raise error
491 elsif (D (1) & D (2)) >= Zu then
494 -- This is the complex case where we definitely have a double digit
495 -- divisor and a dividend of at least three digits. We use the classical
496 -- multiple division algorithm (see section (4.3.1) of Knuth's "The Art
497 -- of Computer Programming", Vol. 2 for a description (algorithm D).
500 -- First normalize the divisor so that it has the leading bit on.
501 -- We do this by finding the appropriate left shift amount.
505 if (Zhi and 16#FFFF0000#) = 0 then
507 Zu := Shift_Left (Zu, 16);
510 if (Hi (Zu) and 16#FF00_0000#) = 0 then
512 Zu := Shift_Left (Zu, 8);
515 if (Hi (Zu) and 16#F000_0000#) = 0 then
517 Zu := Shift_Left (Zu, 4);
520 if (Hi (Zu) and 16#C000_0000#) = 0 then
522 Zu := Shift_Left (Zu, 2);
525 if (Hi (Zu) and 16#8000_0000#) = 0 then
527 Zu := Shift_Left (Zu, 1);
533 -- Note that when we scale up the dividend, it still fits in four
534 -- digits, since we already tested for overflow, and scaling does
535 -- not change the invariant that (D (1) & D (2)) >= Zu.
537 T1 := Shift_Left (D (1) & D (2), Scale);
539 T2 := Shift_Left (0 & D (3), Scale);
540 D (2) := Lo (T1) or Hi (T2);
541 T3 := Shift_Left (0 & D (4), Scale);
542 D (3) := Lo (T2) or Hi (T3);
545 -- Loop to compute quotient digits, runs twice for Qd(1) and Qd(2).
549 -- Compute next quotient digit. We have to divide three digits by
550 -- two digits. We estimate the quotient by dividing the leading
551 -- two digits by the leading digit. Given the scaling we did above
552 -- which ensured the first bit of the divisor is set, this gives
553 -- an estimate of the quotient that is at most two too high.
555 if D (J + 1) = Zhi then
556 Qd (J + 1) := 2 ** 32 - 1;
558 Qd (J + 1) := Lo ((D (J + 1) & D (J + 2)) / Zhi);
561 -- Compute amount to subtract
563 T1 := Qd (J + 1) * Zlo;
564 T2 := Qd (J + 1) * Zhi;
566 T1 := Hi (T1) + Lo (T2);
568 S1 := Hi (T1) + Hi (T2);
570 -- Adjust quotient digit if it was too high
573 exit when Le3 (S1, S2, S3, D (J + 1), D (J + 2), D (J + 3));
574 Qd (J + 1) := Qd (J + 1) - 1;
575 Sub3 (S1, S2, S3, 0, Zhi, Zlo);
578 -- Now subtract S1&S2&S3 from D1&D2&D3 ready for next step
580 Sub3 (D (J + 1), D (J + 2), D (J + 3), S1, S2, S3);
583 -- The two quotient digits are now set, and the remainder of the
584 -- scaled division is in D3&D4. To get the remainder for the
585 -- original unscaled division, we rescale this dividend.
587 -- We rescale the divisor as well, to make the proper comparison
588 -- for rounding below.
590 Qu := Qd (1) & Qd (2);
591 Ru := Shift_Right (D (3) & D (4), Scale);
592 Zu := Shift_Right (Zu, Scale);
595 -- Deal with rounding case
597 if Round and then Ru > (Zu - Uns64'(1)) / Uns64'(2) then
598 Qu := Qu + Uns64 (1);
601 -- Set final signs (RM 4.5.5(27-30))
603 -- Case of dividend (X * Y) sign positive
605 if (X >= 0 and then Y >= 0)
606 or else (X < 0 and then Y < 0)
608 R := To_Pos_Int (Ru);
611 Q := To_Pos_Int (Qu);
613 Q := To_Neg_Int (Qu);
616 -- Case of dividend (X * Y) sign negative
619 R := To_Neg_Int (Ru);
622 Q := To_Neg_Int (Qu);
624 Q := To_Pos_Int (Qu);
633 procedure Sub3 (X1, X2, X3 : in out Uns32; Y1, Y2, Y3 : in Uns32) is
653 -------------------------------
654 -- Subtract_With_Ovflo_Check --
655 -------------------------------
657 function Subtract_With_Ovflo_Check (X, Y : Int64) return Int64 is
658 R : constant Int64 := To_Int (To_Uns (X) - To_Uns (Y));
662 if Y > 0 or else R >= 0 then
667 if Y <= 0 or else R < 0 then
673 end Subtract_With_Ovflo_Check;
679 function To_Neg_Int (A : Uns64) return Int64 is
680 R : constant Int64 := -To_Int (A);
694 function To_Pos_Int (A : Uns64) return Int64 is
695 R : constant Int64 := To_Int (A);