+++ /dev/null
-------------------------------------------------------------------------------
--- --
--- GNAT RUN-TIME COMPONENTS --
--- --
--- S Y S T E M . A R I T H _ 6 4 --
--- --
--- B o d y --
--- --
--- $Revision: 1.2.10.1 $
--- --
--- Copyright (C) 1992-2001 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 GNAT.Exceptions; use GNAT.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;
- function "*" (A : Uns64; B : Uns32) return Uns64;
- pragma Inline ("*");
- -- Length doubling multiplications
-
- 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 : Uns64; B : Uns32) return Uns64 is
- begin
- return 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;
projects / msp430-gcc.git / blobdiff