| ------------------------------------------------------------------------------ |
| -- -- |
| -- GNAT RUN-TIME COMPONENTS -- |
| -- -- |
| -- S Y S T E M . A R I T H _ 6 4 -- |
| -- -- |
| -- B o d y -- |
| -- -- |
| -- Copyright (C) 1992-2002 Free Software Foundation, Inc. -- |
| -- -- |
| -- GNAT is free software; you can redistribute it and/or modify it under -- |
| -- terms of the GNU General Public License as published by the Free Soft- -- |
| -- ware Foundation; either version 2, or (at your option) any later ver- -- |
| -- sion. GNAT is distributed in the hope that it will be useful, but WITH- -- |
| -- OUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY -- |
| -- or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License -- |
| -- for more details. You should have received a copy of the GNU General -- |
| -- Public License distributed with GNAT; see file COPYING. If not, write -- |
| -- to the Free Software Foundation, 59 Temple Place - Suite 330, Boston, -- |
| -- MA 02111-1307, USA. -- |
| -- -- |
| -- As a special exception, if other files instantiate generics from this -- |
| -- unit, or you link this unit with other files to produce an executable, -- |
| -- this unit does not by itself cause the resulting executable to be -- |
| -- covered by the GNU General Public License. This exception does not -- |
| -- however invalidate any other reasons why the executable file might be -- |
| -- covered by the GNU Public License. -- |
| -- -- |
| -- GNAT was originally developed by the GNAT team at New York University. -- |
| -- Extensive contributions were provided by Ada Core Technologies Inc. -- |
| -- -- |
| ------------------------------------------------------------------------------ |
| |
| with System.Pure_Exceptions; use System.Pure_Exceptions; |
| |
| with Interfaces; use Interfaces; |
| with Unchecked_Conversion; |
| |
| package body System.Arith_64 is |
| |
| pragma Suppress (Overflow_Check); |
| pragma Suppress (Range_Check); |
| |
| subtype Uns64 is Unsigned_64; |
| function To_Uns is new Unchecked_Conversion (Int64, Uns64); |
| function To_Int is new Unchecked_Conversion (Uns64, Int64); |
| |
| subtype Uns32 is Unsigned_32; |
| |
| ----------------------- |
| -- Local Subprograms -- |
| ----------------------- |
| |
| function "+" (A, B : Uns32) return Uns64; |
| function "+" (A : Uns64; B : Uns32) return Uns64; |
| pragma Inline ("+"); |
| -- Length doubling additions |
| |
| function "-" (A : Uns64; B : Uns32) return Uns64; |
| pragma Inline ("-"); |
| -- Length doubling subtraction |
| |
| function "*" (A, B : Uns32) return Uns64; |
| pragma Inline ("*"); |
| -- Length doubling multiplication |
| |
| function "/" (A : Uns64; B : Uns32) return Uns64; |
| pragma Inline ("/"); |
| -- Length doubling division |
| |
| function "rem" (A : Uns64; B : Uns32) return Uns64; |
| pragma Inline ("rem"); |
| -- Length doubling remainder |
| |
| function "&" (Hi, Lo : Uns32) return Uns64; |
| pragma Inline ("&"); |
| -- Concatenate hi, lo values to form 64-bit result |
| |
| function Lo (A : Uns64) return Uns32; |
| pragma Inline (Lo); |
| -- Low order half of 64-bit value |
| |
| function Hi (A : Uns64) return Uns32; |
| pragma Inline (Hi); |
| -- High order half of 64-bit value |
| |
| function To_Neg_Int (A : Uns64) return Int64; |
| -- Convert to negative integer equivalent. If the input is in the range |
| -- 0 .. 2 ** 63, then the corresponding negative signed integer (obtained |
| -- by negating the given value) is returned, otherwise constraint error |
| -- is raised. |
| |
| function To_Pos_Int (A : Uns64) return Int64; |
| -- Convert to positive integer equivalent. If the input is in the range |
| -- 0 .. 2 ** 63-1, then the corresponding non-negative signed integer is |
| -- returned, otherwise constraint error is raised. |
| |
| procedure Raise_Error; |
| pragma No_Return (Raise_Error); |
| -- Raise constraint error with appropriate message |
| |
| --------- |
| -- "&" -- |
| --------- |
| |
| function "&" (Hi, Lo : Uns32) return Uns64 is |
| begin |
| return Shift_Left (Uns64 (Hi), 32) or Uns64 (Lo); |
| end "&"; |
| |
| --------- |
| -- "*" -- |
| --------- |
| |
| function "*" (A, B : Uns32) return Uns64 is |
| begin |
| return Uns64 (A) * Uns64 (B); |
| end "*"; |
| |
| --------- |
| -- "+" -- |
| --------- |
| |
| function "+" (A, B : Uns32) return Uns64 is |
| begin |
| return Uns64 (A) + Uns64 (B); |
| end "+"; |
| |
| function "+" (A : Uns64; B : Uns32) return Uns64 is |
| begin |
| return A + Uns64 (B); |
| end "+"; |
| |
| --------- |
| -- "-" -- |
| --------- |
| |
| function "-" (A : Uns64; B : Uns32) return Uns64 is |
| begin |
| return A - Uns64 (B); |
| end "-"; |
| |
| --------- |
| -- "/" -- |
| --------- |
| |
| function "/" (A : Uns64; B : Uns32) return Uns64 is |
| begin |
| return A / Uns64 (B); |
| end "/"; |
| |
| ----------- |
| -- "rem" -- |
| ----------- |
| |
| function "rem" (A : Uns64; B : Uns32) return Uns64 is |
| begin |
| return A rem Uns64 (B); |
| end "rem"; |
| |
| -------------------------- |
| -- Add_With_Ovflo_Check -- |
| -------------------------- |
| |
| function Add_With_Ovflo_Check (X, Y : Int64) return Int64 is |
| R : constant Int64 := To_Int (To_Uns (X) + To_Uns (Y)); |
| |
| begin |
| if X >= 0 then |
| if Y < 0 or else R >= 0 then |
| return R; |
| end if; |
| |
| else -- X < 0 |
| if Y > 0 or else R < 0 then |
| return R; |
| end if; |
| end if; |
| |
| Raise_Error; |
| end Add_With_Ovflo_Check; |
| |
| ------------------- |
| -- Double_Divide -- |
| ------------------- |
| |
| procedure Double_Divide |
| (X, Y, Z : Int64; |
| Q, R : out Int64; |
| Round : Boolean) |
| is |
| Xu : constant Uns64 := To_Uns (abs X); |
| Yu : constant Uns64 := To_Uns (abs Y); |
| |
| Yhi : constant Uns32 := Hi (Yu); |
| Ylo : constant Uns32 := Lo (Yu); |
| |
| Zu : constant Uns64 := To_Uns (abs Z); |
| Zhi : constant Uns32 := Hi (Zu); |
| Zlo : constant Uns32 := Lo (Zu); |
| |
| T1, T2 : Uns64; |
| Du, Qu, Ru : Uns64; |
| Den_Pos : Boolean; |
| |
| begin |
| if Yu = 0 or else Zu = 0 then |
| Raise_Error; |
| end if; |
| |
| -- Compute Y * Z. Note that if the result overflows 64 bits unsigned, |
| -- then the rounded result is clearly zero (since the dividend is at |
| -- most 2**63 - 1, the extra bit of precision is nice here!) |
| |
| if Yhi /= 0 then |
| if Zhi /= 0 then |
| Q := 0; |
| R := X; |
| return; |
| else |
| T2 := Yhi * Zlo; |
| end if; |
| |
| else |
| if Zhi /= 0 then |
| T2 := Ylo * Zhi; |
| else |
| T2 := 0; |
| end if; |
| end if; |
| |
| T1 := Ylo * Zlo; |
| T2 := T2 + Hi (T1); |
| |
| if Hi (T2) /= 0 then |
| Q := 0; |
| R := X; |
| return; |
| end if; |
| |
| Du := Lo (T2) & Lo (T1); |
| Qu := Xu / Du; |
| Ru := Xu rem Du; |
| |
| -- Deal with rounding case |
| |
| if Round and then Ru > (Du - Uns64'(1)) / Uns64'(2) then |
| Qu := Qu + Uns64'(1); |
| end if; |
| |
| -- Set final signs (RM 4.5.5(27-30)) |
| |
| Den_Pos := (Y < 0) = (Z < 0); |
| |
| -- Case of dividend (X) sign positive |
| |
| if X >= 0 then |
| R := To_Int (Ru); |
| |
| if Den_Pos then |
| Q := To_Int (Qu); |
| else |
| Q := -To_Int (Qu); |
| end if; |
| |
| -- Case of dividend (X) sign negative |
| |
| else |
| R := -To_Int (Ru); |
| |
| if Den_Pos then |
| Q := -To_Int (Qu); |
| else |
| Q := To_Int (Qu); |
| end if; |
| end if; |
| end Double_Divide; |
| |
| -------- |
| -- Hi -- |
| -------- |
| |
| function Hi (A : Uns64) return Uns32 is |
| begin |
| return Uns32 (Shift_Right (A, 32)); |
| end Hi; |
| |
| -------- |
| -- Lo -- |
| -------- |
| |
| function Lo (A : Uns64) return Uns32 is |
| begin |
| return Uns32 (A and 16#FFFF_FFFF#); |
| end Lo; |
| |
| ------------------------------- |
| -- Multiply_With_Ovflo_Check -- |
| ------------------------------- |
| |
| function Multiply_With_Ovflo_Check (X, Y : Int64) return Int64 is |
| Xu : constant Uns64 := To_Uns (abs X); |
| Xhi : constant Uns32 := Hi (Xu); |
| Xlo : constant Uns32 := Lo (Xu); |
| |
| Yu : constant Uns64 := To_Uns (abs Y); |
| Yhi : constant Uns32 := Hi (Yu); |
| Ylo : constant Uns32 := Lo (Yu); |
| |
| T1, T2 : Uns64; |
| |
| begin |
| if Xhi /= 0 then |
| if Yhi /= 0 then |
| Raise_Error; |
| else |
| T2 := Xhi * Ylo; |
| end if; |
| |
| elsif Yhi /= 0 then |
| T2 := Xlo * Yhi; |
| |
| else -- Yhi = Xhi = 0 |
| T2 := 0; |
| end if; |
| |
| -- Here we have T2 set to the contribution to the upper half |
| -- of the result from the upper halves of the input values. |
| |
| T1 := Xlo * Ylo; |
| T2 := T2 + Hi (T1); |
| |
| if Hi (T2) /= 0 then |
| Raise_Error; |
| end if; |
| |
| T2 := Lo (T2) & Lo (T1); |
| |
| if X >= 0 then |
| if Y >= 0 then |
| return To_Pos_Int (T2); |
| else |
| return To_Neg_Int (T2); |
| end if; |
| else -- X < 0 |
| if Y < 0 then |
| return To_Pos_Int (T2); |
| else |
| return To_Neg_Int (T2); |
| end if; |
| end if; |
| |
| end Multiply_With_Ovflo_Check; |
| |
| ----------------- |
| -- Raise_Error -- |
| ----------------- |
| |
| procedure Raise_Error is |
| begin |
| Raise_Exception (CE, "64-bit arithmetic overflow"); |
| end Raise_Error; |
| |
| ------------------- |
| -- Scaled_Divide -- |
| ------------------- |
| |
| procedure Scaled_Divide |
| (X, Y, Z : Int64; |
| Q, R : out Int64; |
| Round : Boolean) |
| is |
| Xu : constant Uns64 := To_Uns (abs X); |
| Xhi : constant Uns32 := Hi (Xu); |
| Xlo : constant Uns32 := Lo (Xu); |
| |
| Yu : constant Uns64 := To_Uns (abs Y); |
| Yhi : constant Uns32 := Hi (Yu); |
| Ylo : constant Uns32 := Lo (Yu); |
| |
| Zu : Uns64 := To_Uns (abs Z); |
| Zhi : Uns32 := Hi (Zu); |
| Zlo : Uns32 := Lo (Zu); |
| |
| D1, D2, D3, D4 : Uns32; |
| -- The dividend, four digits (D1 is high order) |
| |
| Q1, Q2 : Uns32; |
| -- The quotient, two digits (Q1 is high order) |
| |
| S1, S2, S3 : Uns32; |
| -- Value to subtract, three digits (S1 is high order) |
| |
| Qu : Uns64; |
| Ru : Uns64; |
| -- Unsigned quotient and remainder |
| |
| Scale : Natural; |
| -- Scaling factor used for multiple-precision divide. Dividend and |
| -- Divisor are multiplied by 2 ** Scale, and the final remainder |
| -- is divided by the scaling factor. The reason for this scaling |
| -- is to allow more accurate estimation of quotient digits. |
| |
| T1, T2, T3 : Uns64; |
| -- Temporary values |
| |
| begin |
| -- First do the multiplication, giving the four digit dividend |
| |
| T1 := Xlo * Ylo; |
| D4 := Lo (T1); |
| D3 := Hi (T1); |
| |
| if Yhi /= 0 then |
| T1 := Xlo * Yhi; |
| T2 := D3 + Lo (T1); |
| D3 := Lo (T2); |
| D2 := Hi (T1) + Hi (T2); |
| |
| if Xhi /= 0 then |
| T1 := Xhi * Ylo; |
| T2 := D3 + Lo (T1); |
| D3 := Lo (T2); |
| T3 := D2 + Hi (T1); |
| T3 := T3 + Hi (T2); |
| D2 := Lo (T3); |
| D1 := Hi (T3); |
| |
| T1 := (D1 & D2) + Uns64'(Xhi * Yhi); |
| D1 := Hi (T1); |
| D2 := Lo (T1); |
| |
| else |
| D1 := 0; |
| end if; |
| |
| else |
| if Xhi /= 0 then |
| T1 := Xhi * Ylo; |
| T2 := D3 + Lo (T1); |
| D3 := Lo (T2); |
| D2 := Hi (T1) + Hi (T2); |
| |
| else |
| D2 := 0; |
| end if; |
| |
| D1 := 0; |
| end if; |
| |
| -- Now it is time for the dreaded multiple precision division. First |
| -- an easy case, check for the simple case of a one digit divisor. |
| |
| if Zhi = 0 then |
| if D1 /= 0 or else D2 >= Zlo then |
| Raise_Error; |
| |
| -- Here we are dividing at most three digits by one digit |
| |
| else |
| T1 := D2 & D3; |
| T2 := Lo (T1 rem Zlo) & D4; |
| |
| Qu := Lo (T1 / Zlo) & Lo (T2 / Zlo); |
| Ru := T2 rem Zlo; |
| end if; |
| |
| -- If divisor is double digit and too large, raise error |
| |
| elsif (D1 & D2) >= Zu then |
| Raise_Error; |
| |
| -- This is the complex case where we definitely have a double digit |
| -- divisor and a dividend of at least three digits. We use the classical |
| -- multiple division algorithm (see section (4.3.1) of Knuth's "The Art |
| -- of Computer Programming", Vol. 2 for a description (algorithm D). |
| |
| else |
| -- First normalize the divisor so that it has the leading bit on. |
| -- We do this by finding the appropriate left shift amount. |
| |
| Scale := 0; |
| |
| if (Zhi and 16#FFFF0000#) = 0 then |
| Scale := 16; |
| Zu := Shift_Left (Zu, 16); |
| end if; |
| |
| if (Hi (Zu) and 16#FF00_0000#) = 0 then |
| Scale := Scale + 8; |
| Zu := Shift_Left (Zu, 8); |
| end if; |
| |
| if (Hi (Zu) and 16#F000_0000#) = 0 then |
| Scale := Scale + 4; |
| Zu := Shift_Left (Zu, 4); |
| end if; |
| |
| if (Hi (Zu) and 16#C000_0000#) = 0 then |
| Scale := Scale + 2; |
| Zu := Shift_Left (Zu, 2); |
| end if; |
| |
| if (Hi (Zu) and 16#8000_0000#) = 0 then |
| Scale := Scale + 1; |
| Zu := Shift_Left (Zu, 1); |
| end if; |
| |
| Zhi := Hi (Zu); |
| Zlo := Lo (Zu); |
| |
| -- Note that when we scale up the dividend, it still fits in four |
| -- digits, since we already tested for overflow, and scaling does |
| -- not change the invariant that (D1 & D2) >= Zu. |
| |
| T1 := Shift_Left (D1 & D2, Scale); |
| D1 := Hi (T1); |
| T2 := Shift_Left (0 & D3, Scale); |
| D2 := Lo (T1) or Hi (T2); |
| T3 := Shift_Left (0 & D4, Scale); |
| D3 := Lo (T2) or Hi (T3); |
| D4 := Lo (T3); |
| |
| -- Compute first quotient digit. We have to divide three digits by |
| -- two digits, and we estimate the quotient by dividing the leading |
| -- two digits by the leading digit. Given the scaling we did above |
| -- which ensured the first bit of the divisor is set, this gives an |
| -- estimate of the quotient that is at most two too high. |
| |
| if D1 = Zhi then |
| Q1 := 2 ** 32 - 1; |
| else |
| Q1 := Lo ((D1 & D2) / Zhi); |
| end if; |
| |
| -- Compute amount to subtract |
| |
| T1 := Q1 * Zlo; |
| T2 := Q1 * Zhi; |
| S3 := Lo (T1); |
| T1 := Hi (T1) + Lo (T2); |
| S2 := Lo (T1); |
| S1 := Hi (T1) + Hi (T2); |
| |
| -- Adjust quotient digit if it was too high |
| |
| loop |
| exit when S1 < D1; |
| |
| if S1 = D1 then |
| exit when S2 < D2; |
| |
| if S2 = D2 then |
| exit when S3 <= D3; |
| end if; |
| end if; |
| |
| Q1 := Q1 - 1; |
| |
| T1 := (S2 & S3) - Zlo; |
| S3 := Lo (T1); |
| T1 := (S1 & S2) - Zhi; |
| S2 := Lo (T1); |
| S1 := Hi (T1); |
| end loop; |
| |
| -- Subtract from dividend (note: do not bother to set D1 to |
| -- zero, since it is no longer needed in the calculation). |
| |
| T1 := (D2 & D3) - S3; |
| D3 := Lo (T1); |
| T1 := (D1 & Hi (T1)) - S2; |
| D2 := Lo (T1); |
| |
| -- Compute second quotient digit in same manner |
| |
| if D2 = Zhi then |
| Q2 := 2 ** 32 - 1; |
| else |
| Q2 := Lo ((D2 & D3) / Zhi); |
| end if; |
| |
| T1 := Q2 * Zlo; |
| T2 := Q2 * Zhi; |
| S3 := Lo (T1); |
| T1 := Hi (T1) + Lo (T2); |
| S2 := Lo (T1); |
| S1 := Hi (T1) + Hi (T2); |
| |
| loop |
| exit when S1 < D2; |
| |
| if S1 = D2 then |
| exit when S2 < D3; |
| |
| if S2 = D3 then |
| exit when S3 <= D4; |
| end if; |
| end if; |
| |
| Q2 := Q2 - 1; |
| |
| T1 := (S2 & S3) - Zlo; |
| S3 := Lo (T1); |
| T1 := (S1 & S2) - Zhi; |
| S2 := Lo (T1); |
| S1 := Hi (T1); |
| end loop; |
| |
| T1 := (D3 & D4) - S3; |
| D4 := Lo (T1); |
| T1 := (D2 & Hi (T1)) - S2; |
| D3 := Lo (T1); |
| |
| -- The two quotient digits are now set, and the remainder of the |
| -- scaled division is in (D3 & D4). To get the remainder for the |
| -- original unscaled division, we rescale this dividend. |
| -- We rescale the divisor as well, to make the proper comparison |
| -- for rounding below. |
| |
| Qu := Q1 & Q2; |
| Ru := Shift_Right (D3 & D4, Scale); |
| Zu := Shift_Right (Zu, Scale); |
| end if; |
| |
| -- Deal with rounding case |
| |
| if Round and then Ru > (Zu - Uns64'(1)) / Uns64'(2) then |
| Qu := Qu + Uns64 (1); |
| end if; |
| |
| -- Set final signs (RM 4.5.5(27-30)) |
| |
| -- Case of dividend (X * Y) sign positive |
| |
| if (X >= 0 and then Y >= 0) |
| or else (X < 0 and then Y < 0) |
| then |
| R := To_Pos_Int (Ru); |
| |
| if Z > 0 then |
| Q := To_Pos_Int (Qu); |
| else |
| Q := To_Neg_Int (Qu); |
| end if; |
| |
| -- Case of dividend (X * Y) sign negative |
| |
| else |
| R := To_Neg_Int (Ru); |
| |
| if Z > 0 then |
| Q := To_Neg_Int (Qu); |
| else |
| Q := To_Pos_Int (Qu); |
| end if; |
| end if; |
| |
| end Scaled_Divide; |
| |
| ------------------------------- |
| -- Subtract_With_Ovflo_Check -- |
| ------------------------------- |
| |
| function Subtract_With_Ovflo_Check (X, Y : Int64) return Int64 is |
| R : constant Int64 := To_Int (To_Uns (X) - To_Uns (Y)); |
| |
| begin |
| if X >= 0 then |
| if Y > 0 or else R >= 0 then |
| return R; |
| end if; |
| |
| else -- X < 0 |
| if Y <= 0 or else R < 0 then |
| return R; |
| end if; |
| end if; |
| |
| Raise_Error; |
| end Subtract_With_Ovflo_Check; |
| |
| ---------------- |
| -- To_Neg_Int -- |
| ---------------- |
| |
| function To_Neg_Int (A : Uns64) return Int64 is |
| R : constant Int64 := -To_Int (A); |
| |
| begin |
| if R <= 0 then |
| return R; |
| else |
| Raise_Error; |
| end if; |
| end To_Neg_Int; |
| |
| ---------------- |
| -- To_Pos_Int -- |
| ---------------- |
| |
| function To_Pos_Int (A : Uns64) return Int64 is |
| R : constant Int64 := To_Int (A); |
| |
| begin |
| if R >= 0 then |
| return R; |
| else |
| Raise_Error; |
| end if; |
| end To_Pos_Int; |
| |
| end System.Arith_64; |
