-- --
-- B o d y --
-- --
--- $Revision: 1.16 $
--- --
--- Copyright (C) 1992-2001 Free Software Foundation, Inc. --
+-- Copyright (C) 1992-2009, 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- --
+-- ware Foundation; either version 3, 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. --
+-- or FITNESS FOR A PARTICULAR PURPOSE. --
+-- --
+-- As a special exception under Section 7 of GPL version 3, you are granted --
+-- additional permissions described in the GCC Runtime Library Exception, --
+-- version 3.1, as published by the Free Software Foundation. --
-- --
--- 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. --
+-- You should have received a copy of the GNU General Public License and --
+-- a copy of the GCC Runtime Library Exception along with this program; --
+-- see the files COPYING3 and COPYING.RUNTIME respectively. If not, see --
+-- <http://www.gnu.org/licenses/>. --
-- --
-- GNAT was originally developed by the GNAT team at New York University. --
--- It is now maintained by Ada Core Technologies Inc (http://www.gnat.com). --
+-- Extensive contributions were provided by Ada Core Technologies Inc. --
-- --
------------------------------------------------------------------------------
-with GNAT.Exceptions; use GNAT.Exceptions;
-
with Interfaces; use Interfaces;
-with Unchecked_Conversion;
+with Ada.Unchecked_Conversion;
package body System.Arith_64 is
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);
+ function To_Uns is new Ada.Unchecked_Conversion (Int64, Uns64);
+ function To_Int is new Ada.Unchecked_Conversion (Uns64, Int64);
subtype Uns32 is Unsigned_32;
pragma Inline ("+");
-- Length doubling additions
- function "-" (A : Uns64; B : Uns32) return Uns64;
- pragma Inline ("-");
- -- Length doubling subtraction
-
function "*" (A, B : Uns32) return Uns64;
- function "*" (A : Uns64; B : Uns32) return Uns64;
pragma Inline ("*");
- -- Length doubling multiplications
+ -- Length doubling multiplication
function "/" (A : Uns64; B : Uns32) return Uns64;
pragma Inline ("/");
pragma Inline ("&");
-- Concatenate hi, lo values to form 64-bit result
+ function Le3 (X1, X2, X3 : Uns32; Y1, Y2, Y3 : Uns32) return Boolean;
+ -- Determines if 96 bit value X1&X2&X3 <= Y1&Y2&Y3
+
function Lo (A : Uns64) return Uns32;
pragma Inline (Lo);
-- Low order half of 64-bit value
pragma Inline (Hi);
-- High order half of 64-bit value
+ procedure Sub3 (X1, X2, X3 : in out Uns32; Y1, Y2, Y3 : Uns32);
+ -- Computes X1&X2&X3 := X1&X2&X3 - Y1&Y1&Y3 with mod 2**96 wrap
+
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
return Uns64 (A) * Uns64 (B);
end "*";
- function "*" (A : Uns64; B : Uns32) return Uns64 is
- begin
- return A * Uns64 (B);
- end "*";
-
---------
-- "+" --
---------
end "+";
---------
- -- "-" --
- ---------
-
- function "-" (A : Uns64; B : Uns32) return Uns64 is
- begin
- return A - Uns64 (B);
- end "-";
-
- ---------
-- "/" --
---------
end if;
else
- if Zhi /= 0 then
- T2 := Ylo * Zhi;
- else
- T2 := 0;
- end if;
+ T2 := (if Zhi /= 0 then Ylo * Zhi else 0);
end if;
T1 := Ylo * Zlo;
end if;
Du := Lo (T2) & Lo (T1);
+
+ -- Set final signs (RM 4.5.5(27-30))
+
+ Den_Pos := (Y < 0) = (Z < 0);
+
+ -- Check overflow case of largest negative number divided by 1
+
+ if X = Int64'First and then Du = 1 and then not Den_Pos then
+ Raise_Error;
+ end if;
+
+ -- Perform the actual division
+
Qu := Xu / Du;
Ru := Xu rem Du;
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;
+ Q := (if Den_Pos then To_Int (Qu) else -To_Int (Qu));
-- 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;
+ Q := (if Den_Pos then -To_Int (Qu) else To_Int (Qu));
end if;
end Double_Divide;
return Uns32 (Shift_Right (A, 32));
end Hi;
+ ---------
+ -- Le3 --
+ ---------
+
+ function Le3 (X1, X2, X3 : Uns32; Y1, Y2, Y3 : Uns32) return Boolean is
+ begin
+ if X1 < Y1 then
+ return True;
+ elsif X1 > Y1 then
+ return False;
+ elsif X2 < Y2 then
+ return True;
+ elsif X2 > Y2 then
+ return False;
+ else
+ return X3 <= Y3;
+ end if;
+ end Le3;
+
--------
-- Lo --
--------
T2 := Xhi * Ylo;
end if;
- else
- if Yhi /= 0 then
- T2 := Xlo * Yhi;
- else
- return X * Y;
- 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);
procedure Raise_Error is
begin
- Raise_Exception (CE, "64-bit arithmetic overflow");
+ raise Constraint_Error with "64-bit arithmetic overflow";
end Raise_Error;
-------------------
Zhi : Uns32 := Hi (Zu);
Zlo : Uns32 := Lo (Zu);
- D1, D2, D3, D4 : Uns32;
- -- The dividend, four digits (D1 is high order)
+ D : array (1 .. 4) of Uns32;
+ -- The dividend, four digits (D(1) is high order)
- Q1, Q2 : Uns32;
- -- The quotient, two digits (Q1 is high order)
+ Qd : array (1 .. 2) of Uns32;
+ -- The quotient digits, two digits (Qd(1) is high order)
S1, S2, S3 : Uns32;
-- Value to subtract, three digits (S1 is high order)
-- First do the multiplication, giving the four digit dividend
T1 := Xlo * Ylo;
- D4 := Lo (T1);
- D3 := Hi (T1);
+ D (4) := Lo (T1);
+ D (3) := Hi (T1);
if Yhi /= 0 then
T1 := Xlo * Yhi;
- T2 := D3 + Lo (T1);
- D3 := Lo (T2);
- D2 := Hi (T1) + Hi (T2);
+ T2 := D (3) + Lo (T1);
+ D (3) := Lo (T2);
+ D (2) := Hi (T1) + Hi (T2);
if Xhi /= 0 then
T1 := Xhi * Ylo;
- T2 := D3 + Lo (T1);
- D3 := Lo (T2);
- T3 := D2 + Hi (T1);
+ T2 := D (3) + Lo (T1);
+ D (3) := Lo (T2);
+ T3 := D (2) + Hi (T1);
T3 := T3 + Hi (T2);
- D2 := Lo (T3);
- D1 := Hi (T3);
+ D (2) := Lo (T3);
+ D (1) := Hi (T3);
- T1 := (D1 & D2) + Uns64'(Xhi * Yhi);
- D1 := Hi (T1);
- D2 := Lo (T1);
+ T1 := (D (1) & D (2)) + Uns64'(Xhi * Yhi);
+ D (1) := Hi (T1);
+ D (2) := Lo (T1);
else
- D1 := 0;
+ D (1) := 0;
end if;
else
if Xhi /= 0 then
T1 := Xhi * Ylo;
- T2 := D3 + Lo (T1);
- D3 := Lo (T2);
- D2 := Hi (T1) + Hi (T2);
+ T2 := D (3) + Lo (T1);
+ D (3) := Lo (T2);
+ D (2) := Hi (T1) + Hi (T2);
else
- D2 := 0;
+ D (2) := 0;
end if;
- D1 := 0;
+ D (1) := 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
+ if D (1) /= 0 or else D (2) >= 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;
+ T1 := D (2) & D (3);
+ T2 := Lo (T1 rem Zlo) & D (4);
Qu := Lo (T1 / Zlo) & Lo (T2 / Zlo);
Ru := T2 rem Zlo;
-- If divisor is double digit and too large, raise error
- elsif (D1 & D2) >= Zu then
+ elsif (D (1) & D (2)) >= 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
+ -- multiple division algorithm (see section (4.3.1) of Knuth's "The Art
-- of Computer Programming", Vol. 2 for a description (algorithm D).
else
-- 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
+ -- not change the invariant that (D (1) & D (2)) >= Zu.
- T1 := Q1 * Zlo;
- T2 := Q1 * Zhi;
- S3 := Lo (T1);
- T1 := Hi (T1) + Lo (T2);
- S2 := Lo (T1);
- S1 := Hi (T1) + Hi (T2);
+ T1 := Shift_Left (D (1) & D (2), Scale);
+ D (1) := Hi (T1);
+ T2 := Shift_Left (0 & D (3), Scale);
+ D (2) := Lo (T1) or Hi (T2);
+ T3 := Shift_Left (0 & D (4), Scale);
+ D (3) := Lo (T2) or Hi (T3);
+ D (4) := Lo (T3);
- -- Adjust quotient digit if it was too high
+ -- Loop to compute quotient digits, runs twice for Qd(1) and Qd(2)
- loop
- exit when S1 < D1;
+ for J in 0 .. 1 loop
- if S1 = D1 then
- exit when S2 < D2;
+ -- Compute next quotient digit. We have to divide three digits by
+ -- two digits. 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 S2 = D2 then
- exit when S3 <= D3;
- end if;
- end if;
+ Qd (J + 1) := (if D (J + 1) = Zhi
+ then 2 ** 32 - 1
+ else Lo ((D (J + 1) & D (J + 2)) / Zhi));
- Q1 := Q1 - 1;
+ -- Compute amount to subtract
- T1 := (S2 & S3) - Zlo;
+ T1 := Qd (J + 1) * Zlo;
+ T2 := Qd (J + 1) * Zhi;
S3 := Lo (T1);
- T1 := (S1 & S2) - Zhi;
+ T1 := Hi (T1) + Lo (T2);
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;
+ S1 := Hi (T1) + Hi (T2);
- T1 := Q2 * Zlo;
- T2 := Q2 * 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 < D2;
+ loop
+ exit when Le3 (S1, S2, S3, D (J + 1), D (J + 2), D (J + 3));
+ Qd (J + 1) := Qd (J + 1) - 1;
+ Sub3 (S1, S2, S3, 0, Zhi, Zlo);
+ end loop;
- if S1 = D2 then
- exit when S2 < D3;
+ -- Now subtract S1&S2&S3 from D1&D2&D3 ready for next step
- 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);
+ Sub3 (D (J + 1), D (J + 2), D (J + 3), S1, S2, S3);
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
+ -- 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);
+ Qu := Qd (1) & Qd (2);
+ Ru := Shift_Right (D (3) & D (4), Scale);
Zu := Shift_Right (Zu, Scale);
end if;
-- Case of dividend (X * Y) sign positive
- if (X >= 0 and then Y >= 0)
- or else (X < 0 and then Y < 0)
- then
+ 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;
+ Q := (if Z > 0 then To_Pos_Int (Qu) else To_Neg_Int (Qu));
-- Case of dividend (X * Y) sign negative
else
R := To_Neg_Int (Ru);
+ Q := (if Z > 0 then To_Neg_Int (Qu) else To_Pos_Int (Qu));
+ end if;
+ end Scaled_Divide;
- if Z > 0 then
- Q := To_Neg_Int (Qu);
- else
- Q := To_Pos_Int (Qu);
+ ----------
+ -- Sub3 --
+ ----------
+
+ procedure Sub3 (X1, X2, X3 : in out Uns32; Y1, Y2, Y3 : Uns32) is
+ begin
+ if Y3 > X3 then
+ if X2 = 0 then
+ X1 := X1 - 1;
end if;
+
+ X2 := X2 - 1;
end if;
- end Scaled_Divide;
+ X3 := X3 - Y3;
+
+ if Y2 > X2 then
+ X1 := X1 - 1;
+ end if;
+
+ X2 := X2 - Y2;
+ X1 := X1 - Y1;
+ end Sub3;
-------------------------------
-- Subtract_With_Ovflo_Check --
OSDN Git Service
