Subversion Repositories openrisc

------------------------------------------------------------------------------

--                                                                          --

--                         GNAT RUN-TIME COMPONENTS                         --

--                                                                          --

--                      S Y S T E M . A R I T H _ 6 4                       --

--                                                                          --

--                                 B o d y                                  --

--                                                                          --

--          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 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.                                     --

--                                                                          --

-- 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.               --

--                                                                          --

-- 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. --

-- Extensive contributions were provided by Ada Core Technologies Inc.      --

--                                                                          --

------------------------------------------------------------------------------

with Interfaces; use Interfaces;

with Ada.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 Ada.Unchecked_Conversion (Int64, Uns64);

   function To_Int is new Ada.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, 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 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

   function Hi (A : Uns64) return Uns32;

   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

   --  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 "/";

   -----------

   -- "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

         T2 := (if Zhi /= 0 then Ylo * Zhi else 0);

      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);

      --  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;

      --  Deal with rounding case

      if Round and then Ru > (Du - Uns64'(1)) / Uns64'(2) then

         Qu := Qu + Uns64'(1);

      end if;

      --  Case of dividend (X) sign positive

      if X >= 0 then

         R := To_Int (Ru);

         Q := (if Den_Pos then To_Int (Qu) else -To_Int (Qu));

      --  Case of dividend (X) sign negative

      else

         R := -To_Int (Ru);

         Q := (if Den_Pos then -To_Int (Qu) else To_Int (Qu));

      end if;

   end Double_Divide;

   --------

   -- Hi --

   --------

   function Hi (A : Uns64) return Uns32 is

   begin

      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 --

   --------

   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 Constraint_Error with "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);

      D : array (1 .. 4) of Uns32;

      --  The dividend, four digits (D(1) 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)

      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;

      D (4) := Lo (T1);

      D (3) := Hi (T1);

      if Yhi /= 0 then

         T1 := Xlo * Yhi;

         T2 := D (3) + Lo (T1);

         D (3) := Lo (T2);

         D (2) := Hi (T1) + Hi (T2);

         if Xhi /= 0 then

            T1 := Xhi * Ylo;

            T2 := D (3) + Lo (T1);

            D (3) := Lo (T2);

            T3 := D (2) + Hi (T1);

            T3 := T3 + Hi (T2);

            D (2) := Lo (T3);

            D (1) := Hi (T3);

            T1 := (D (1) & D (2)) + Uns64'(Xhi * Yhi);

            D (1) := Hi (T1);

            D (2) := Lo (T1);

         else

            D (1) := 0;

         end if;

      else

         if Xhi /= 0 then

            T1 := Xhi * Ylo;

            T2 := D (3) + Lo (T1);

            D (3) := Lo (T2);

            D (2) := Hi (T1) + Hi (T2);

         else

            D (2) := 0;

         end if;

         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 D (1) /= 0 or else D (2) >= Zlo then

            Raise_Error;

         --  Here we are dividing at most three digits by one digit

         else

            T1 := D (2) & D (3);

            T2 := Lo (T1 rem Zlo) & D (4);

            Qu := Lo (T1 / Zlo) & Lo (T2 / Zlo);

            Ru := T2 rem Zlo;

         end if;

      --  If divisor is double digit and too large, raise error

      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

      --  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 (D (1) & D (2)) >= Zu.

         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);

         --  Loop to compute quotient digits, runs twice for Qd(1) and Qd(2)

         for J in 0 .. 1 loop

            --  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.

            Qd (J + 1) := (if D (J + 1) = Zhi

                           then 2 ** 32 - 1

                           else Lo ((D (J + 1) & D (J + 2)) / Zhi));

            --  Compute amount to subtract

            T1 := Qd (J + 1) * Zlo;

            T2 := Qd (J + 1) * 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 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;

            --  Now subtract S1&S2&S3 from D1&D2&D3 ready for next step

            Sub3 (D (J + 1), D (J + 2), D (J + 3), S1, S2, S3);

         end loop;

         --  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 := Qd (1) & Qd (2);

         Ru := Shift_Right (D (3) & D (4), 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);

         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;

   ----------

   -- 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;

      X3 := X3 - Y3;

      if Y2 > X2 then

         X1 := X1 - 1;

      end if;

      X2 := X2 - Y2;

      X1 := X1 - Y1;

   end Sub3;

   -------------------------------

   -- 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;