Subversion Repositories openrisc

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

--                                                                          --

--                         GNAT COMPILER COMPONENTS                         --

--                                                                          --

--                             E V A L _ F A T                              --

--                                                                          --

--                                 B o d y                                  --

--                                                                          --

--          Copyright (C) 1992-2010, 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.  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 COPYING3.  If not, go to --

-- http://www.gnu.org/licenses for a complete copy of the license.          --

--                                                                          --

-- GNAT was originally developed  by the GNAT team at  New York University. --

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

--                                                                          --

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

with Einfo;    use Einfo;

with Errout;   use Errout;

with Targparm; use Targparm;

package body Eval_Fat is

   Radix : constant Int := 2;

   --  This code is currently only correct for the radix 2 case. We use the

   --  symbolic value Radix where possible to help in the unlikely case of

   --  anyone ever having to adjust this code for another value, and for

   --  documentation purposes.

   --  Another assumption is that the range of the floating-point type is

   --  symmetric around zero.

   type Radix_Power_Table is array (Int range 1 .. 4) of Int;

   Radix_Powers : constant Radix_Power_Table :=

                    (Radix ** 1, Radix ** 2, Radix ** 3, Radix ** 4);

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

   -- Local Subprograms --

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

   procedure Decompose

     (RT       : R;

      X        : T;

      Fraction : out T;

      Exponent : out UI;

      Mode     : Rounding_Mode := Round);

   --  Decomposes a non-zero floating-point number into fraction and exponent

   --  parts. The fraction is in the interval 1.0 / Radix .. T'Pred (1.0) and

   --  uses Rbase = Radix. The result is rounded to a nearest machine number.

   procedure Decompose_Int

     (RT       : R;

      X        : T;

      Fraction : out UI;

      Exponent : out UI;

      Mode     : Rounding_Mode);

   --  This is similar to Decompose, except that the Fraction value returned

   --  is an integer representing the value Fraction * Scale, where Scale is

   --  the value (Machine_Radix_Value (RT) ** Machine_Mantissa_Value (RT)). The

   --  value is obtained by using biased rounding (halfway cases round away

   --  from zero), round to even, a floor operation or a ceiling operation

   --  depending on the setting of Mode (see corresponding descriptions in

   --  Urealp).

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

   -- Adjacent --

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

   function Adjacent (RT : R; X, Towards : T) return T is

   begin

      if Towards = X then

         return X;

      elsif Towards > X then

         return Succ (RT, X);

      else

         return Pred (RT, X);

      end if;

   end Adjacent;

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

   -- Ceiling --

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

   function Ceiling (RT : R; X : T) return T is

      XT : constant T := Truncation (RT, X);

   begin

      if UR_Is_Negative (X) then

         return XT;

      elsif X = XT then

         return X;

      else

         return XT + Ureal_1;

      end if;

   end Ceiling;

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

   -- Compose --

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

   function Compose (RT : R; Fraction : T; Exponent : UI) return T is

      Arg_Frac : T;

      Arg_Exp  : UI;

      pragma Warnings (Off, Arg_Exp);

   begin

      Decompose (RT, Fraction, Arg_Frac, Arg_Exp);

      return Scaling (RT, Arg_Frac, Exponent);

   end Compose;

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

   -- Copy_Sign --

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

   function Copy_Sign (RT : R; Value, Sign : T) return T is

      pragma Warnings (Off, RT);

      Result : T;

   begin

      Result := abs Value;

      if UR_Is_Negative (Sign) then

         return -Result;

      else

         return Result;

      end if;

   end Copy_Sign;

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

   -- Decompose --

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

   procedure Decompose

     (RT       : R;

      X        : T;

      Fraction : out T;

      Exponent : out UI;

      Mode     : Rounding_Mode := Round)

   is

      Int_F : UI;

   begin

      Decompose_Int (RT, abs X, Int_F, Exponent, Mode);

      Fraction := UR_From_Components

       (Num      => Int_F,

        Den      => Machine_Mantissa_Value (RT),

        Rbase    => Radix,

        Negative => False);

      if UR_Is_Negative (X) then

         Fraction := -Fraction;

      end if;

      return;

   end Decompose;

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

   -- Decompose_Int --

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

   --  This procedure should be modified with care, as there are many non-

   --  obvious details that may cause problems that are hard to detect. For

   --  zero arguments, Fraction and Exponent are set to zero. Note that sign

   --  of zero cannot be preserved.

   procedure Decompose_Int

     (RT       : R;

      X        : T;

      Fraction : out UI;

      Exponent : out UI;

      Mode     : Rounding_Mode)

   is

      Base : Int := Rbase (X);

      N    : UI  := abs Numerator (X);

      D    : UI  := Denominator (X);

      N_Times_Radix : UI;

      Even : Boolean;

      --  True iff Fraction is even

      Most_Significant_Digit : constant UI :=

                                 Radix ** (Machine_Mantissa_Value (RT) - 1);

      Uintp_Mark : Uintp.Save_Mark;

      --  The code is divided into blocks that systematically release

      --  intermediate values (this routine generates lots of junk!)

   begin

      if N = Uint_0 then

         Fraction := Uint_0;

         Exponent := Uint_0;

         return;

      end if;

      Calculate_D_And_Exponent_1 : begin

         Uintp_Mark := Mark;

         Exponent := Uint_0;

         --  In cases where Base > 1, the actual denominator is Base**D. For

         --  cases where Base is a power of Radix, use the value 1 for the

         --  Denominator and adjust the exponent.

         --  Note: Exponent has different sign from D, because D is a divisor

         for Power in 1 .. Radix_Powers'Last loop

            if Base = Radix_Powers (Power) then

               Exponent := -D * Power;

               Base := 0;

               D := Uint_1;

               exit;

            end if;

         end loop;

         Release_And_Save (Uintp_Mark, D, Exponent);

      end Calculate_D_And_Exponent_1;

      if Base > 0 then

         Calculate_Exponent : begin

            Uintp_Mark := Mark;

            --  For bases that are a multiple of the Radix, divide the base by

            --  Radix and adjust the Exponent. This will help because D will be

            --  much smaller and faster to process.

            --  This occurs for decimal bases on machines with binary floating-

            --  point for example. When calculating 1E40, with Radix = 2, N

            --  will be 93 bits instead of 133.

            --        N            E

            --      ------  * Radix

            --           D

            --       Base

            --                  N                        E

            --    =  --------------------------  *  Radix

            --                     D        D

            --         (Base/Radix)  * Radix

            --             N                  E-D

            --    =  ---------------  *  Radix

            --                    D

            --        (Base/Radix)

            --  This code is commented out, because it causes numerous

            --  failures in the regression suite. To be studied ???

            while False and then Base > 0 and then Base mod Radix = 0 loop

               Base := Base / Radix;

               Exponent := Exponent + D;

            end loop;

            Release_And_Save (Uintp_Mark, Exponent);

         end Calculate_Exponent;

         --  For remaining bases we must actually compute the exponentiation

         --  Because the exponentiation can be negative, and D must be integer,

         --  the numerator is corrected instead.

         Calculate_N_And_D : begin

            Uintp_Mark := Mark;

            if D < 0 then

               N := N * Base ** (-D);

               D := Uint_1;

            else

               D := Base ** D;

            end if;

            Release_And_Save (Uintp_Mark, N, D);

         end Calculate_N_And_D;

         Base := 0;

      end if;

      --  Now scale N and D so that N / D is a value in the interval [1.0 /

      --  Radix, 1.0) and adjust Exponent accordingly, so the value N / D *

      --  Radix ** Exponent remains unchanged.

      --  Step 1 - Adjust N so N / D >= 1 / Radix, or N = 0

      --  N and D are positive, so N / D >= 1 / Radix implies N * Radix >= D.

      --  As this scaling is not possible for N is Uint_0, zero is handled

      --  explicitly at the start of this subprogram.

      Calculate_N_And_Exponent : begin

         Uintp_Mark := Mark;

         N_Times_Radix := N * Radix;

         while not (N_Times_Radix >= D) loop

            N := N_Times_Radix;

            Exponent := Exponent - 1;

            N_Times_Radix := N * Radix;

         end loop;

         Release_And_Save (Uintp_Mark, N, Exponent);

      end Calculate_N_And_Exponent;

      --  Step 2 - Adjust D so N / D < 1

      --  Scale up D so N / D < 1, so N < D

      Calculate_D_And_Exponent_2 : begin

         Uintp_Mark := Mark;

         while not (N < D) loop

            --  As N / D >= 1, N / (D * Radix) will be at least 1 / Radix, so

            --  the result of Step 1 stays valid

            D := D * Radix;

            Exponent := Exponent + 1;

         end loop;

         Release_And_Save (Uintp_Mark, D, Exponent);

      end Calculate_D_And_Exponent_2;

      --  Here the value N / D is in the range [1.0 / Radix .. 1.0)

      --  Now find the fraction by doing a very simple-minded division until

      --  enough digits have been computed.

      --  This division works for all radices, but is only efficient for a

      --  binary radix. It is just like a manual division algorithm, but

      --  instead of moving the denominator one digit right, we move the

      --  numerator one digit left so the numerator and denominator remain

      --  integral.

      Fraction := Uint_0;

      Even := True;

      Calculate_Fraction_And_N : begin

         Uintp_Mark := Mark;

         loop

            while N >= D loop

               N := N - D;

               Fraction := Fraction + 1;

               Even := not Even;

            end loop;

            --  Stop when the result is in [1.0 / Radix, 1.0)

            exit when Fraction >= Most_Significant_Digit;

            N := N * Radix;

            Fraction := Fraction * Radix;

            Even := True;

         end loop;

         Release_And_Save (Uintp_Mark, Fraction, N);

      end Calculate_Fraction_And_N;

      Calculate_Fraction_And_Exponent : begin

         Uintp_Mark := Mark;

         --  Determine correct rounding based on the remainder which is in

         --  N and the divisor D. The rounding is performed on the absolute

         --  value of X, so Ceiling and Floor need to check for the sign of

         --  X explicitly.

         case Mode is

            when Round_Even =>

               --  This rounding mode should not be used for static

               --  expressions, but only for compile-time evaluation of

               --  non-static expressions.

               if (Even and then N * 2 > D)

                     or else

                  (not Even and then N * 2 >= D)

               then

                  Fraction := Fraction + 1;

               end if;

            when Round   =>

               --  Do not round to even as is done with IEEE arithmetic, but

               --  instead round away from zero when the result is exactly

               --  between two machine numbers. See RM 4.9(38).

               if N * 2 >= D then

                  Fraction := Fraction + 1;

               end if;

            when Ceiling =>

               if N > Uint_0 and then not UR_Is_Negative (X) then

                  Fraction := Fraction + 1;

               end if;

            when Floor   =>

               if N > Uint_0 and then UR_Is_Negative (X) then

                  Fraction := Fraction + 1;

               end if;

         end case;

         --  The result must be normalized to [1.0/Radix, 1.0), so adjust if

         --  the result is 1.0 because of rounding.

         if Fraction = Most_Significant_Digit * Radix then

            Fraction := Most_Significant_Digit;

            Exponent := Exponent + 1;

         end if;

         --  Put back sign after applying the rounding

         if UR_Is_Negative (X) then

            Fraction := -Fraction;

         end if;

         Release_And_Save (Uintp_Mark, Fraction, Exponent);

      end Calculate_Fraction_And_Exponent;

   end Decompose_Int;

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

   -- Exponent --

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

   function Exponent (RT : R; X : T) return UI is

      X_Frac : UI;

      X_Exp  : UI;

      pragma Warnings (Off, X_Frac);

   begin

      Decompose_Int (RT, X, X_Frac, X_Exp, Round_Even);

      return X_Exp;

   end Exponent;

   -----------

   -- Floor --

   -----------

   function Floor (RT : R; X : T) return T is

      XT : constant T := Truncation (RT, X);

   begin

      if UR_Is_Positive (X) then

         return XT;

      elsif XT = X then

         return X;

      else

         return XT - Ureal_1;

      end if;

   end Floor;

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

   -- Fraction --

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

   function Fraction (RT : R; X : T) return T is

      X_Frac : T;

      X_Exp  : UI;

      pragma Warnings (Off, X_Exp);

   begin

      Decompose (RT, X, X_Frac, X_Exp);

      return X_Frac;

   end Fraction;

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

   -- Leading_Part --

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

   function Leading_Part (RT : R; X : T; Radix_Digits : UI) return T is

      RD : constant UI := UI_Min (Radix_Digits, Machine_Mantissa_Value (RT));

      L  : UI;

      Y  : T;

   begin

      L := Exponent (RT, X) - RD;

      Y := UR_From_Uint (UR_Trunc (Scaling (RT, X, -L)));

      return Scaling (RT, Y, L);

   end Leading_Part;

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

   -- Machine --

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

   function Machine

     (RT    : R;

      X     : T;

      Mode  : Rounding_Mode;

      Enode : Node_Id) return T

   is

      X_Frac : T;

      X_Exp  : UI;

      Emin   : constant UI := Machine_Emin_Value (RT);

   begin

      Decompose (RT, X, X_Frac, X_Exp, Mode);

      --  Case of denormalized number or (gradual) underflow

      --  A denormalized number is one with the minimum exponent Emin, but that

      --  breaks the assumption that the first digit of the mantissa is a one.

      --  This allows the first non-zero digit to be in any of the remaining

      --  Mant - 1 spots. The gap between subsequent denormalized numbers is

      --  the same as for the smallest normalized numbers. However, the number

      --  of significant digits left decreases as a result of the mantissa now

      --  having leading seros.

      if X_Exp < Emin then

         declare

            Emin_Den : constant UI := Machine_Emin_Value (RT)

                                        - Machine_Mantissa_Value (RT) + Uint_1;

         begin

            if X_Exp < Emin_Den or not Denorm_On_Target then

               if UR_Is_Negative (X) then

                  Error_Msg_N

                    ("floating-point value underflows to -0.0?", Enode);

                  return Ureal_M_0;

               else

                  Error_Msg_N

                    ("floating-point value underflows to 0.0?", Enode);

                  return Ureal_0;

               end if;

            elsif Denorm_On_Target then

               --  Emin - Mant <= X_Exp < Emin, so result is denormal. Handle

               --  gradual underflow by first computing the number of

               --  significant bits still available for the mantissa and

               --  then truncating the fraction to this number of bits.

               --  If this value is different from the original fraction,

               --  precision is lost due to gradual underflow.

               --  We probably should round here and prevent double rounding as

               --  a result of first rounding to a model number and then to a

               --  machine number. However, this is an extremely rare case that

               --  is not worth the extra complexity. In any case, a warning is

               --  issued in cases where gradual underflow occurs.

               declare

                  Denorm_Sig_Bits : constant UI := X_Exp - Emin_Den + 1;

                  X_Frac_Denorm   : constant T := UR_From_Components

                    (UR_Trunc (Scaling (RT, abs X_Frac, Denorm_Sig_Bits)),

                     Denorm_Sig_Bits,

                     Radix,

                     UR_Is_Negative (X));

               begin

                  if X_Frac_Denorm /= X_Frac then

                     Error_Msg_N

                       ("gradual underflow causes loss of precision?",

                        Enode);

                     X_Frac := X_Frac_Denorm;

                  end if;

               end;

            end if;

         end;

      end if;

      return Scaling (RT, X_Frac, X_Exp);

   end Machine;

   -----------

   -- Model --

   -----------

   function Model (RT : R; X : T) return T is

      X_Frac : T;

      X_Exp  : UI;

   begin

      Decompose (RT, X, X_Frac, X_Exp);

      return Compose (RT, X_Frac, X_Exp);

   end Model;

   ----------

   -- Pred --

   ----------

   function Pred (RT : R; X : T) return T is

   begin

      return -Succ (RT, -X);

   end Pred;

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

   -- Remainder --

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

   function Remainder (RT : R; X, Y : T) return T is

      A        : T;

      B        : T;

      Arg      : T;

      P        : T;

      Arg_Frac : T;

      P_Frac   : T;

      Sign_X   : T;

      IEEE_Rem : T;

      Arg_Exp  : UI;

      P_Exp    : UI;

      K        : UI;

      P_Even   : Boolean;

      pragma Warnings (Off, Arg_Frac);

   begin

      if UR_Is_Positive (X) then

         Sign_X :=  Ureal_1;

      else

         Sign_X := -Ureal_1;

      end if;

      Arg := abs X;

      P   := abs Y;

      if Arg < P then

         P_Even := True;

         IEEE_Rem := Arg;

         P_Exp := Exponent (RT, P);

      else

         --  ??? what about zero cases?

         Decompose (RT, Arg, Arg_Frac, Arg_Exp);

         Decompose (RT, P,   P_Frac,   P_Exp);

         P := Compose (RT, P_Frac, Arg_Exp);

         K := Arg_Exp - P_Exp;

         P_Even := True;

         IEEE_Rem := Arg;

         for Cnt in reverse 0 .. UI_To_Int (K) loop

            if IEEE_Rem >= P then

               P_Even := False;

               IEEE_Rem := IEEE_Rem - P;

            else

               P_Even := True;

            end if;

            P := P * Ureal_Half;

         end loop;

      end if;

      --  That completes the calculation of modulus remainder. The final step

      --  is get the IEEE remainder. Here we compare Rem with (abs Y) / 2.

      if P_Exp >= 0 then

         A := IEEE_Rem;

         B := abs Y * Ureal_Half;

      else

         A := IEEE_Rem * Ureal_2;

         B := abs Y;

      end if;

      if A > B or else (A = B and then not P_Even) then

         IEEE_Rem := IEEE_Rem - abs Y;

      end if;

      return Sign_X * IEEE_Rem;

   end Remainder;

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

   -- Rounding --

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

   function Rounding (RT : R; X : T) return T is

      Result : T;

      Tail   : T;

   begin

      Result := Truncation (RT, abs X);

      Tail   := abs X - Result;

      if Tail >= Ureal_Half  then

         Result := Result + Ureal_1;

      end if;

      if UR_Is_Negative (X) then

         return -Result;

      else

         return Result;

      end if;

   end Rounding;

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

   -- Scaling --

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

   function Scaling (RT : R; X : T; Adjustment : UI) return T is

      pragma Warnings (Off, RT);

   begin

      if Rbase (X) = Radix then

         return UR_From_Components

           (Num      => Numerator (X),

            Den      => Denominator (X) - Adjustment,

            Rbase    => Radix,

            Negative => UR_Is_Negative (X));

      elsif Adjustment >= 0 then

         return X * Radix ** Adjustment;

      else

         return X / Radix ** (-Adjustment);

      end if;

   end Scaling;

   ----------

   -- Succ --

   ----------

   function Succ (RT : R; X : T) return T is

      Emin     : constant UI := Machine_Emin_Value (RT);

      Mantissa : constant UI := Machine_Mantissa_Value (RT);

      Exp      : UI := UI_Max (Emin, Exponent (RT, X));

      Frac     : T;

      New_Frac : T;

   begin

      if UR_Is_Zero (X) then

         Exp := Emin;

      end if;

      --  Set exponent such that the radix point will be directly following the

      --  mantissa after scaling.

      if Denorm_On_Target or Exp /= Emin then

         Exp := Exp - Mantissa;

      else

         Exp := Exp - 1;

      end if;

      Frac := Scaling (RT, X, -Exp);

      New_Frac := Ceiling (RT, Frac);

      if New_Frac = Frac then

         if New_Frac = Scaling (RT, -Ureal_1, Mantissa - 1) then

            New_Frac := New_Frac + Scaling (RT, Ureal_1, Uint_Minus_1);

         else

            New_Frac := New_Frac + Ureal_1;

         end if;

      end if;

      return Scaling (RT, New_Frac, Exp);

   end Succ;

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

   -- Truncation --

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

   function Truncation (RT : R; X : T) return T is

      pragma Warnings (Off, RT);

   begin

      return UR_From_Uint (UR_Trunc (X));

   end Truncation;

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

   -- Unbiased_Rounding --

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

   function Unbiased_Rounding (RT : R; X : T) return T is

      Abs_X  : constant T := abs X;

      Result : T;

      Tail   : T;

   begin

      Result := Truncation (RT, Abs_X);

      Tail   := Abs_X - Result;

      if Tail > Ureal_Half  then

         Result := Result + Ureal_1;

      elsif Tail = Ureal_Half then

         Result := Ureal_2 *

                     Truncation (RT, (Result / Ureal_2) + Ureal_Half);

      end if;

      if UR_Is_Negative (X) then

         return -Result;

      elsif UR_Is_Positive (X) then

         return Result;

      --  For zero case, make sure sign of zero is preserved

      else

         return X;

      end if;

   end Unbiased_Rounding;

end Eval_Fat;