1 ------------------------------------------------------------------------------
3 -- GNAT RUN-TIME COMPONENTS --
5 -- S Y S T E M . A R I T H _ 6 4 --
9 -- Copyright (C) 1992-2002 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, 59 Temple Place - Suite 330, Boston, --
20 -- MA 02111-1307, 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 : Uns64; B : Uns32) return Uns64;
61 -- Length doubling subtraction
63 function "*" (A, B : Uns32) return Uns64;
65 -- Length doubling multiplication
67 function "/" (A : Uns64; B : Uns32) return Uns64;
69 -- Length doubling division
71 function "rem" (A : Uns64; B : Uns32) return Uns64;
72 pragma Inline ("rem");
73 -- Length doubling remainder
75 function "&" (Hi, Lo : Uns32) return Uns64;
77 -- Concatenate hi, lo values to form 64-bit result
79 function Lo (A : Uns64) return Uns32;
81 -- Low order half of 64-bit value
83 function Hi (A : Uns64) return Uns32;
85 -- High order half of 64-bit value
87 function To_Neg_Int (A : Uns64) return Int64;
88 -- Convert to negative integer equivalent. If the input is in the range
89 -- 0 .. 2 ** 63, then the corresponding negative signed integer (obtained
90 -- by negating the given value) is returned, otherwise constraint error
93 function To_Pos_Int (A : Uns64) return Int64;
94 -- Convert to positive integer equivalent. If the input is in the range
95 -- 0 .. 2 ** 63-1, then the corresponding non-negative signed integer is
96 -- returned, otherwise constraint error is raised.
98 procedure Raise_Error;
99 pragma No_Return (Raise_Error);
100 -- Raise constraint error with appropriate message
106 function "&" (Hi, Lo : Uns32) return Uns64 is
108 return Shift_Left (Uns64 (Hi), 32) or Uns64 (Lo);
115 function "*" (A, B : Uns32) return Uns64 is
117 return Uns64 (A) * Uns64 (B);
124 function "+" (A, B : Uns32) return Uns64 is
126 return Uns64 (A) + Uns64 (B);
129 function "+" (A : Uns64; B : Uns32) return Uns64 is
131 return A + Uns64 (B);
138 function "-" (A : Uns64; B : Uns32) return Uns64 is
140 return A - Uns64 (B);
147 function "/" (A : Uns64; B : Uns32) return Uns64 is
149 return A / Uns64 (B);
156 function "rem" (A : Uns64; B : Uns32) return Uns64 is
158 return A rem Uns64 (B);
161 --------------------------
162 -- Add_With_Ovflo_Check --
163 --------------------------
165 function Add_With_Ovflo_Check (X, Y : Int64) return Int64 is
166 R : constant Int64 := To_Int (To_Uns (X) + To_Uns (Y));
170 if Y < 0 or else R >= 0 then
175 if Y > 0 or else R < 0 then
181 end Add_With_Ovflo_Check;
187 procedure Double_Divide
192 Xu : constant Uns64 := To_Uns (abs X);
193 Yu : constant Uns64 := To_Uns (abs Y);
195 Yhi : constant Uns32 := Hi (Yu);
196 Ylo : constant Uns32 := Lo (Yu);
198 Zu : constant Uns64 := To_Uns (abs Z);
199 Zhi : constant Uns32 := Hi (Zu);
200 Zlo : constant Uns32 := Lo (Zu);
207 if Yu = 0 or else Zu = 0 then
211 -- Compute Y * Z. Note that if the result overflows 64 bits unsigned,
212 -- then the rounded result is clearly zero (since the dividend is at
213 -- most 2**63 - 1, the extra bit of precision is nice here!)
241 Du := Lo (T2) & Lo (T1);
245 -- Deal with rounding case
247 if Round and then Ru > (Du - Uns64'(1)) / Uns64'(2) then
248 Qu := Qu + Uns64'(1);
251 -- Set final signs (RM 4.5.5(27-30))
253 Den_Pos := (Y < 0) = (Z < 0);
255 -- Case of dividend (X) sign positive
266 -- Case of dividend (X) sign negative
283 function Hi (A : Uns64) return Uns32 is
285 return Uns32 (Shift_Right (A, 32));
292 function Lo (A : Uns64) return Uns32 is
294 return Uns32 (A and 16#FFFF_FFFF#);
297 -------------------------------
298 -- Multiply_With_Ovflo_Check --
299 -------------------------------
301 function Multiply_With_Ovflo_Check (X, Y : Int64) return Int64 is
302 Xu : constant Uns64 := To_Uns (abs X);
303 Xhi : constant Uns32 := Hi (Xu);
304 Xlo : constant Uns32 := Lo (Xu);
306 Yu : constant Uns64 := To_Uns (abs Y);
307 Yhi : constant Uns32 := Hi (Yu);
308 Ylo : constant Uns32 := Lo (Yu);
323 else -- Yhi = Xhi = 0
327 -- Here we have T2 set to the contribution to the upper half
328 -- of the result from the upper halves of the input values.
337 T2 := Lo (T2) & Lo (T1);
341 return To_Pos_Int (T2);
343 return To_Neg_Int (T2);
347 return To_Pos_Int (T2);
349 return To_Neg_Int (T2);
353 end Multiply_With_Ovflo_Check;
359 procedure Raise_Error is
361 Raise_Exception (CE, "64-bit arithmetic overflow");
368 procedure Scaled_Divide
373 Xu : constant Uns64 := To_Uns (abs X);
374 Xhi : constant Uns32 := Hi (Xu);
375 Xlo : constant Uns32 := Lo (Xu);
377 Yu : constant Uns64 := To_Uns (abs Y);
378 Yhi : constant Uns32 := Hi (Yu);
379 Ylo : constant Uns32 := Lo (Yu);
381 Zu : Uns64 := To_Uns (abs Z);
382 Zhi : Uns32 := Hi (Zu);
383 Zlo : Uns32 := Lo (Zu);
385 D1, D2, D3, D4 : Uns32;
386 -- The dividend, four digits (D1 is high order)
389 -- The quotient, two digits (Q1 is high order)
392 -- Value to subtract, three digits (S1 is high order)
396 -- Unsigned quotient and remainder
399 -- Scaling factor used for multiple-precision divide. Dividend and
400 -- Divisor are multiplied by 2 ** Scale, and the final remainder
401 -- is divided by the scaling factor. The reason for this scaling
402 -- is to allow more accurate estimation of quotient digits.
408 -- First do the multiplication, giving the four digit dividend
418 D2 := Hi (T1) + Hi (T2);
429 T1 := (D1 & D2) + Uns64'(Xhi * Yhi);
442 D2 := Hi (T1) + Hi (T2);
451 -- Now it is time for the dreaded multiple precision division. First
452 -- an easy case, check for the simple case of a one digit divisor.
455 if D1 /= 0 or else D2 >= Zlo then
458 -- Here we are dividing at most three digits by one digit
462 T2 := Lo (T1 rem Zlo) & D4;
464 Qu := Lo (T1 / Zlo) & Lo (T2 / Zlo);
468 -- If divisor is double digit and too large, raise error
470 elsif (D1 & D2) >= Zu then
473 -- This is the complex case where we definitely have a double digit
474 -- divisor and a dividend of at least three digits. We use the classical
475 -- multiple division algorithm (see section (4.3.1) of Knuth's "The Art
476 -- of Computer Programming", Vol. 2 for a description (algorithm D).
479 -- First normalize the divisor so that it has the leading bit on.
480 -- We do this by finding the appropriate left shift amount.
484 if (Zhi and 16#FFFF0000#) = 0 then
486 Zu := Shift_Left (Zu, 16);
489 if (Hi (Zu) and 16#FF00_0000#) = 0 then
491 Zu := Shift_Left (Zu, 8);
494 if (Hi (Zu) and 16#F000_0000#) = 0 then
496 Zu := Shift_Left (Zu, 4);
499 if (Hi (Zu) and 16#C000_0000#) = 0 then
501 Zu := Shift_Left (Zu, 2);
504 if (Hi (Zu) and 16#8000_0000#) = 0 then
506 Zu := Shift_Left (Zu, 1);
512 -- Note that when we scale up the dividend, it still fits in four
513 -- digits, since we already tested for overflow, and scaling does
514 -- not change the invariant that (D1 & D2) >= Zu.
516 T1 := Shift_Left (D1 & D2, Scale);
518 T2 := Shift_Left (0 & D3, Scale);
519 D2 := Lo (T1) or Hi (T2);
520 T3 := Shift_Left (0 & D4, Scale);
521 D3 := Lo (T2) or Hi (T3);
524 -- Compute first quotient digit. We have to divide three digits by
525 -- two digits, and we estimate the quotient by dividing the leading
526 -- two digits by the leading digit. Given the scaling we did above
527 -- which ensured the first bit of the divisor is set, this gives an
528 -- estimate of the quotient that is at most two too high.
533 Q1 := Lo ((D1 & D2) / Zhi);
536 -- Compute amount to subtract
541 T1 := Hi (T1) + Lo (T2);
543 S1 := Hi (T1) + Hi (T2);
545 -- Adjust quotient digit if it was too high
560 T1 := (S2 & S3) - Zlo;
562 T1 := (S1 & S2) - Zhi;
567 -- Subtract from dividend (note: do not bother to set D1 to
568 -- zero, since it is no longer needed in the calculation).
570 T1 := (D2 & D3) - S3;
572 T1 := (D1 & Hi (T1)) - S2;
575 -- Compute second quotient digit in same manner
580 Q2 := Lo ((D2 & D3) / Zhi);
586 T1 := Hi (T1) + Lo (T2);
588 S1 := Hi (T1) + Hi (T2);
603 T1 := (S2 & S3) - Zlo;
605 T1 := (S1 & S2) - Zhi;
610 T1 := (D3 & D4) - S3;
612 T1 := (D2 & Hi (T1)) - S2;
615 -- The two quotient digits are now set, and the remainder of the
616 -- scaled division is in (D3 & D4). To get the remainder for the
617 -- original unscaled division, we rescale this dividend.
618 -- We rescale the divisor as well, to make the proper comparison
619 -- for rounding below.
622 Ru := Shift_Right (D3 & D4, Scale);
623 Zu := Shift_Right (Zu, Scale);
626 -- Deal with rounding case
628 if Round and then Ru > (Zu - Uns64'(1)) / Uns64'(2) then
629 Qu := Qu + Uns64 (1);
632 -- Set final signs (RM 4.5.5(27-30))
634 -- Case of dividend (X * Y) sign positive
636 if (X >= 0 and then Y >= 0)
637 or else (X < 0 and then Y < 0)
639 R := To_Pos_Int (Ru);
642 Q := To_Pos_Int (Qu);
644 Q := To_Neg_Int (Qu);
647 -- Case of dividend (X * Y) sign negative
650 R := To_Neg_Int (Ru);
653 Q := To_Neg_Int (Qu);
655 Q := To_Pos_Int (Qu);
661 -------------------------------
662 -- Subtract_With_Ovflo_Check --
663 -------------------------------
665 function Subtract_With_Ovflo_Check (X, Y : Int64) return Int64 is
666 R : constant Int64 := To_Int (To_Uns (X) - To_Uns (Y));
670 if Y > 0 or else R >= 0 then
675 if Y <= 0 or else R < 0 then
681 end Subtract_With_Ovflo_Check;
687 function To_Neg_Int (A : Uns64) return Int64 is
688 R : constant Int64 := -To_Int (A);
702 function To_Pos_Int (A : Uns64) return Int64 is
703 R : constant Int64 := To_Int (A);