[Ada] Implement part of System.Fat_Gen more efficiently

gcc/ada/

	* libgnat/s-fatgen.ads (Compose): Add pragma Inline.
	(Copy_Sign): Likewise.
	(Exponent): Likewise.
	(Fraction): Likewise.
	* libgnat/s-fatgen.adb: Remove with clause for System, add
	with and use clauses for System.Unsigned_Types.
	Add pragma Warnings (Off) for non-static constants.
	Remove precomputed tables of powers of radix and add a few
	constants describing the floating-point format.
	(Gradual_Scaling): Delete.
	(Copy_Sign): Reimplement directly.
	(Decompose): Likewise.
	(Scaling): Likewise.
	(Pred): Speed up.
	(Succ): Likewise.
	(Truncation): Tidy up.
	(Valid): Move constants to library level.

2 files changed, 267 insertions, 329 deletions

diff --git a/gcc/ada/libgnat/s-fatgen.adb b/gcc/ada/libgnat/s-fatgen.adb
index e590198..bebe737 100644
--- a/gcc/ada/libgnat/s-fatgen.adb
+++ b/gcc/ada/libgnat/s-fatgen.adb
@@ -29,51 +29,124 @@
 --                                                                          --
 ------------------------------------------------------------------------------
 
---  The implementation here is portable to any IEEE implementation. It does
---  not handle nonbinary radix, and also assumes that model numbers and
---  machine numbers are basically identical, which is not true of all possible
---  floating-point implementations. On a non-IEEE machine, this body must be
---  specialized appropriately, or better still, its generic instantiations
---  should be replaced by efficient machine-specific code.
+--  This implementation is portable to any IEEE implementation. It does not
+--  handle nonbinary radix, and also assumes that model numbers and machine
+--  numbers are basically identical, which is not true of all possible
+--  floating-point implementations.
 
 with Ada.Unchecked_Conversion;
-with System;
+with System.Unsigned_Types; use System.Unsigned_Types;
 
-package body System.Fat_Gen is
-
-   Float_Radix        : constant T := T (T'Machine_Radix);
-   Radix_To_M_Minus_1 : constant T := Float_Radix ** (T'Machine_Mantissa - 1);
+pragma Warnings (Off, "non-static constant in preelaborated unit");
+--  Every constant is static given our instantiation model
 
+package body System.Fat_Gen is
    pragma Assert (T'Machine_Radix = 2);
    --  This version does not handle radix 16
 
-   --  Constants for Decompose and Scaling
+   Rad : constant T := T (T'Machine_Radix);
+   --  Renaming for the machine radix
 
-   Rad    : constant T := T (T'Machine_Radix);
-   Invrad : constant T := 1.0 / Rad;
+   Mantissa : constant Integer := T'Machine_Mantissa;
+   --  Renaming for the machine mantissa
 
-   subtype Expbits is Integer range 0 .. 6;
-   --  2 ** (2 ** 7) might overflow.  How big can radix-16 exponents get?
-
-   Log_Power : constant array (Expbits) of Integer := (1, 2, 4, 8, 16, 32, 64);
-
-   R_Power : constant array (Expbits) of T :=
-     (Rad **  1,
-      Rad **  2,
-      Rad **  4,
-      Rad **  8,
-      Rad ** 16,
-      Rad ** 32,
-      Rad ** 64);
-
-   R_Neg_Power : constant array (Expbits) of T :=
-     (Invrad **  1,
-      Invrad **  2,
-      Invrad **  4,
-      Invrad **  8,
-      Invrad ** 16,
-      Invrad ** 32,
-      Invrad ** 64);
+   Invrad : constant T := 1.0 / Rad;
+   --  Smallest positive mantissa in the canonical form (RM A.5.3(4))
+
+   Small : constant T := Rad ** (T'Machine_Emin - 1);
+   pragma Unreferenced (Small);
+   --  Smallest positive normalized number
+
+   Tiny : constant T := Rad ** (T'Machine_Emin - Mantissa);
+   --  Smallest positive denormalized number
+
+   RM1 : constant T := Rad ** (Mantissa - 1);
+   --  Smallest positive member of the large consecutive integers. It is equal
+   --  to the ratio Small / Tiny, which means that multiplying by it normalizes
+   --  any nonzero denormalized number.
+
+   IEEE_Emin : constant Integer := T'Machine_Emin - 1;
+   IEEE_Emax : constant Integer := T'Machine_Emax - 1;
+   --  The mantissa is a fraction with first digit set in Ada whereas it is
+   --  shifted by 1 digit to the left in the IEEE floating-point format.
+
+   subtype IEEE_Erange is Integer range IEEE_Emin - 1 .. IEEE_Emax + 1;
+   --  The IEEE floating-point format extends the machine range by 1 to the
+   --  left for denormalized numbers and 1 to the right for infinities/NaNs.
+
+   IEEE_Ebias : constant Integer := -(IEEE_Emin - 1);
+   --  The exponent is biased such that denormalized numbers have it zero
+
+   --  The implementation uses a representation type Float_Rep that allows
+   --  direct access to exponent and mantissa of the floating point number.
+
+   --  The Float_Rep type is a simple array of Float_Word elements. This
+   --  representation is chosen to make it possible to size the type based
+   --  on a generic parameter. Since the array size is known at compile
+   --  time, efficient code can still be generated. The size of Float_Word
+   --  elements should be large enough to allow accessing the exponent in
+   --  one read, but small enough so that all floating-point object sizes
+   --  are a multiple of Float_Word'Size.
+
+   --  The following conditions must be met for all possible instantiations
+   --  of the attribute package:
+
+   --    - T'Size is an integral multiple of Float_Word'Size
+
+   --    - The exponent and sign are completely contained in a single
+   --      component of Float_Rep, named Most Significant Word (MSW).
+
+   --    - The sign occupies the most significant bit of the MSW and the
+   --      exponent is in the following bits. The exception is 80-bit
+   --      double extended, where they occupy the low 16-bit halfword.
+
+   --  The low-level primitives Copy_Sign, Decompose, Scaling and Valid are
+   --  implemented by accessing the bit pattern of the floating-point number.
+   --  Only the normalization of denormalized numbers, if any, and the gradual
+   --  underflow are left to the hardware, mainly because there is some leeway
+   --  for the hardware implementation in this area: for example, the MSB of
+   --  the mantissa, which is 1 for normalized numbers and 0 for denormalized
+   --  numbers, may or may not be stored by the hardware.
+
+   Siz : constant := (if System.Word_Size > 32 then 32 else System.Word_Size);
+   type Float_Word is mod 2**Siz;
+
+   N : constant Natural := (T'Size + Siz - 1) / Siz;
+   Rep_Last : constant Natural := Natural'Min (N - 1, (Mantissa + 16) / Siz);
+   --  Determine the number of Float_Words needed for representing the
+   --  entire floating-point value. Do not take into account excessive
+   --  padding, as occurs on IA-64 where 80 bits floats get padded to 128
+   --  bits. In general, the exponent field cannot be larger than 15 bits,
+   --  even for 128-bit floating-point types, so the final format size
+   --  won't be larger than Mantissa + 16.
+
+   type Float_Rep is array (Natural range 0 .. N - 1) of Float_Word;
+   pragma Suppress_Initialization (Float_Rep);
+   --  This pragma suppresses the generation of an initialization procedure
+   --  for type Float_Rep when operating in Initialize/Normalize_Scalars mode.
+
+   MSW : constant Natural := Rep_Last * Standard'Default_Bit_Order;
+   --  Finding the location of the Exponent_Word is a bit tricky. In general
+   --  we assume Word_Order = Bit_Order.
+
+   Exp_Factor : constant Float_Word :=
+                  (if Mantissa = 64
+                   then 1
+                   else 2**(Siz - 1) / Float_Word (IEEE_Emax - IEEE_Emin + 3));
+   --  Factor that the extracted exponent needs to be divided by to be in
+   --  range 0 .. IEEE_Emax - IEEE_Emin + 2. The special case is 80-bit
+   --  double extended, where the exponent starts the 3rd float word.
+
+   Exp_Mask : constant Float_Word :=
+                Float_Word (IEEE_Emax - IEEE_Emin + 2) * Exp_Factor;
+   --  Value needed to mask out the exponent field. This assumes that the
+   --  range 0 .. IEEE_Emax - IEEE_Emin + 2 contains 2**N values, for some
+   --  N in Natural.
+
+   Sign_Mask : constant Float_Word :=
+                 (if Mantissa = 64 then 2**15 else 2**(Siz - 1));
+   --  Value needed to mask out the sign field. The special case is 80-bit
+   --  double extended, where the exponent starts the 3rd float word.
 
    -----------------------
    -- Local Subprograms --
@@ -85,11 +158,6 @@ package body System.Fat_Gen is
    --  the sign of the exponent. The absolute value of Frac is in the range
    --  0.0 <= Frac < 1.0. If Frac = 0.0 or -0.0, then Expo is always zero.
 
-   function Gradual_Scaling (Adjustment : UI) return T;
-   --  Like Scaling with a first argument of 1.0, but returns the smallest
-   --  denormal rather than zero when the adjustment is smaller than
-   --  Machine_Emin. Used for Succ and Pred.
-
    --------------
    -- Adjacent --
    --------------
@@ -139,19 +207,22 @@ package body System.Fat_Gen is
    ---------------
 
    function Copy_Sign (Value, Sign : T) return T is
-      Result : T;
+      S : constant T := T'Machine (Sign);
 
-      function Is_Negative (V : T) return Boolean;
-      pragma Import (Intrinsic, Is_Negative);
+      Rep_S : Float_Rep;
+      for Rep_S'Address use S'Address;
+      --  Rep_S is a view of the Sign parameter
 
-   begin
-      Result := abs Value;
+      V : T := T'Machine (Value);
 
-      if Is_Negative (Sign) then
-         return -Result;
-      else
-         return Result;
-      end if;
+      Rep_V : Float_Rep;
+      for Rep_V'Address use V'Address;
+      --  Rep_V is a view of the Value parameter
+
+   begin
+      Rep_V (MSW) :=
+        (Rep_V (MSW) and not Sign_Mask) or (Rep_S (MSW) and Sign_Mask);
+      return V;
    end Copy_Sign;
 
    ---------------
@@ -159,94 +230,53 @@ package body System.Fat_Gen is
    ---------------
 
    procedure Decompose (XX : T; Frac : out T; Expo : out UI) is
-      X : constant T := T'Machine (XX);
-
-   begin
-      if X = 0.0 then
-
-         --  The normalized exponent of zero is zero, see RM A.5.2(15)
+      X : T := T'Machine (XX);
 
-         Frac := X;
-         Expo := 0;
-
-      --  Check for infinities, transfinites, whatnot
+      Rep : Float_Rep;
+      for Rep'Address use X'Address;
+      --  Rep is a view of the input floating-point parameter
 
-      elsif X > T'Safe_Last then
-         Frac := Invrad;
-         pragma Annotate (CodePeer, Intentional, "dead code",
-                          "Check float range.");
-         Expo := T'Machine_Emax + 1;
+      Exp : constant IEEE_Erange :=
+              Integer ((Rep (MSW) and Exp_Mask) / Exp_Factor) - IEEE_Ebias;
+      --  Mask/Shift X to only get bits from the exponent. Then convert biased
+      --  value to final value.
 
-      elsif X < T'Safe_First then
-         Frac := -Invrad;
-         pragma Annotate (CodePeer, Intentional, "dead code",
-                          "Check float range.");
-         Expo := T'Machine_Emax + 2;    -- how many extra negative values?
+      Minus : constant Boolean := (Rep (MSW) and Sign_Mask) /= 0;
+      --  Mask/Shift X to only get bit from the sign
 
-      else
-         --  Case of nonzero finite x. Essentially, we just multiply
-         --  by Rad ** (+-2**N) to reduce the range.
+   begin
+      --  The normalized exponent of zero is zero, see RM A.5.3(15)
 
-         declare
-            Ax : T  := abs X;
-            Ex : UI := 0;
+      if X = 0.0 then
+         Expo := 0;
+         Frac := X;
 
-         --  Ax * Rad ** Ex is invariant
+      --  Check for infinities and NaNs
 
-         begin
-            if Ax >= 1.0 then
-               while Ax >= R_Power (Expbits'Last) loop
-                  Ax := Ax * R_Neg_Power (Expbits'Last);
-                  Ex := Ex + Log_Power (Expbits'Last);
-               end loop;
+      elsif Exp = IEEE_Emax + 1 then
+         Expo := T'Machine_Emax + 1;
+         Frac := (if Minus then -Invrad else Invrad);
 
-               --  Ax < Rad ** 64
+      --  Check for nonzero denormalized numbers
 
-               for N in reverse Expbits'First .. Expbits'Last - 1 loop
-                  if Ax >= R_Power (N) then
-                     Ax := Ax * R_Neg_Power (N);
-                     Ex := Ex + Log_Power (N);
-                  end if;
+      elsif Exp = IEEE_Emin - 1 then
+         --  Normalize by multiplying by Radix ** (Mantissa - 1)
 
-                  --  Ax < R_Power (N)
+         Decompose (X * RM1, Frac, Expo);
+         Expo := Expo - (Mantissa - 1);
 
-               end loop;
+      --  Case of normalized numbers
 
-               --  1 <= Ax < Rad
+      else
+         --  The Ada exponent is the IEEE exponent plus 1, see above
 
-               Ax := Ax * Invrad;
-               Ex := Ex + 1;
+         Expo := Exp + 1;
 
-            else
-               --  0 < ax < 1
-
-               while Ax < R_Neg_Power (Expbits'Last) loop
-                  Ax := Ax * R_Power (Expbits'Last);
-                  pragma Annotate (CodePeer, Intentional, "dead code",
-                                   "Check float range.");
-                  Ex := Ex - Log_Power (Expbits'Last);
-               end loop;
-               pragma Annotate
-                 (CodePeer, Intentional,
-                  "test always false",
-                  "expected for some instantiations");
-
-               --  Rad ** -64 <= Ax < 1
-
-               for N in reverse Expbits'First .. Expbits'Last - 1 loop
-                  if Ax < R_Neg_Power (N) then
-                     Ax := Ax * R_Power (N);
-                     Ex := Ex - Log_Power (N);
-                  end if;
-
-                  --  R_Neg_Power (N) <= Ax < 1
-
-               end loop;
-            end if;
+         --  Set Ada exponent of X to zero, so we end up with the fraction
 
-            Frac := (if X > 0.0 then Ax else -Ax);
-            Expo := Ex;
-         end;
+         Rep (MSW) := (Rep (MSW) and not Exp_Mask) +
+                        Float_Word (IEEE_Ebias - 1) * Exp_Factor;
+         Frac := X;
       end if;
    end Decompose;
 
@@ -292,38 +322,6 @@ package body System.Fat_Gen is
       return X_Frac;
    end Fraction;
 
-   ---------------------
-   -- Gradual_Scaling --
-   ---------------------
-
-   function Gradual_Scaling  (Adjustment : UI) return T is
-      Y  : T;
-      Y1 : T;
-      Ex : UI := Adjustment;
-
-   begin
-      if Adjustment < T'Machine_Emin - 1 then
-         Y  := 2.0 ** T'Machine_Emin;
-         Y1 := Y;
-         Ex := Ex - T'Machine_Emin;
-         while Ex < 0 loop
-            Y := T'Machine (Y / 2.0);
-
-            if Y = 0.0 then
-               return Y1;
-            end if;
-
-            Ex := Ex + 1;
-            Y1 := Y;
-         end loop;
-
-         return Y1;
-
-      else
-         return Scaling (1.0, Adjustment);
-      end if;
-   end Gradual_Scaling;
-
    ------------------
    -- Leading_Part --
    ------------------
@@ -333,7 +331,7 @@ package body System.Fat_Gen is
       Y, Z : T;
 
    begin
-      if Radix_Digits >= T'Machine_Mantissa then
+      if Radix_Digits >= Mantissa then
          return X;
 
       elsif Radix_Digits <= 0 then
@@ -420,12 +418,11 @@ package body System.Fat_Gen is
       --  Zero has to be treated specially, since its exponent is zero
 
       if X = 0.0 then
-         return -Succ (X);
+         return -Tiny;
 
-      --  Special treatment for most negative number
+      --  Special treatment for largest negative number: raise Constraint_Error
 
       elsif X = T'First then
-
          raise Constraint_Error with "Pred of largest negative number";
 
       --  For infinities, return unchanged
@@ -439,28 +436,33 @@ package body System.Fat_Gen is
 
       --  Subtract from the given number a number equivalent to the value
       --  of its least significant bit. Given that the most significant bit
-      --  represents a value of 1.0 * radix ** (exp - 1), the value we want
-      --  is obtained by shifting this by (mantissa-1) bits to the right,
+      --  represents a value of 1.0 * Radix ** (Exp - 1), the value we want
+      --  is obtained by shifting this by (Mantissa-1) bits to the right,
       --  i.e. decreasing the exponent by that amount.
 
       else
          Decompose (X, X_Frac, X_Exp);
 
-         --  A special case, if the number we had was a positive power of
-         --  two, then we want to subtract half of what we would otherwise
-         --  subtract, since the exponent is going to be reduced.
+         --  For a denormalized number or a normalized number with the lowest
+         --  exponent, just subtract the Tiny.
+
+         if X_Exp <= T'Machine_Emin then
+            return X - Tiny;
 
-         --  Note that X_Frac has the same sign as X, so if X_Frac is 0.5,
-         --  then we know that we have a positive number (and hence a
-         --  positive power of 2).
+         --  A special case, if the number we had was a power of two on the
+         --  positive side of zero, then we want to subtract half of what we
+         --  would have subtracted, since the exponent is going to be reduced.
 
-         if X_Frac = 0.5 then
-            return X - Gradual_Scaling (X_Exp - T'Machine_Mantissa - 1);
+         --  Note that X_Frac has the same sign as X so, if X_Frac is Invrad,
+         --  then we know that we had a power of two on the positive side.
 
-         --  Otherwise the exponent is unchanged
+         elsif X_Frac = Invrad then
+            return X - Scaling (1.0, X_Exp - Mantissa - 1);
+
+         --  Otherwise the adjustment is unchanged
 
          else
-            return X - Gradual_Scaling (X_Exp - T'Machine_Mantissa);
+            return X - Scaling (1.0, X_Exp - Mantissa);
          end if;
       end if;
    end Pred;
@@ -580,70 +582,90 @@ package body System.Fat_Gen is
    -- Scaling --
    -------------
 
-   --  Return x * rad ** adjustment quickly, or quietly underflow to zero,
-   --  or overflow naturally.
-
    function Scaling (X : T; Adjustment : UI) return T is
+      pragma Assert (Mantissa <= 64);
+      --  This implementation handles only 80-bit IEEE Extended or smaller
+
+      XX : T := T'Machine (X);
+
+      Rep : Float_Rep;
+      for Rep'Address use XX'Address;
+      --  Rep is a view of the input floating-point parameter
+
+      Exp : constant IEEE_Erange :=
+              Integer ((Rep (MSW) and Exp_Mask) / Exp_Factor) - IEEE_Ebias;
+      --  Mask/Shift X to only get bits from the exponent. Then convert biased
+      --  value to final value.
+
+      Minus : constant Boolean := (Rep (MSW) and Sign_Mask) /= 0;
+      --  Mask/Shift X to only get bit from the sign
+
+      Expi, Expf : IEEE_Erange;
+
    begin
-      if X = 0.0 or else Adjustment = 0 then
+      --  Check for zero, infinities, NaNs as well as no adjustment
+
+      if X = 0.0 or else Exp = IEEE_Emax + 1 or else Adjustment = 0 then
          return X;
-      end if;
 
-      --  Nonzero x essentially, just multiply repeatedly by Rad ** (+-2**n)
+      --  Check for nonzero denormalized numbers
 
-      declare
-         Y  : T  := X;
-         Ex : UI := Adjustment;
+      elsif Exp = IEEE_Emin - 1 then
+         --  Check for zero result to protect the subtraction below
 
-      --  Y * Rad ** Ex is invariant
+         if Adjustment < -(Mantissa - 1) then
+            XX := 0.0;
+            return (if Minus then -XX else XX);
 
-      begin
-         if Ex < 0 then
-            while Ex <= -Log_Power (Expbits'Last) loop
-               Y := Y * R_Neg_Power (Expbits'Last);
-               Ex := Ex + Log_Power (Expbits'Last);
-            end loop;
+         --  Normalize by multiplying by Radix ** (Mantissa - 1)
 
-            --  -64 < Ex <= 0
+         else
+            return Scaling (XX * RM1, Adjustment - (Mantissa - 1));
+         end if;
 
-            for N in reverse Expbits'First .. Expbits'Last - 1 loop
-               if Ex <= -Log_Power (N) then
-                  Y := Y * R_Neg_Power (N);
-                  Ex := Ex + Log_Power (N);
-               end if;
+      --  Case of normalized numbers
 
-               --  -Log_Power (N) < Ex <= 0
+      else
+         --  Check for overflow
 
-            end loop;
+         if Adjustment > IEEE_Emax - Exp then
+            XX := 0.0;
+            return (if Minus then -1.0 / XX else 1.0 / XX);
 
-            --  Ex = 0
+         --  Check for underflow
 
-         else
-            --  Ex >= 0
+         elsif Adjustment < IEEE_Emin - Exp then
+            --  Check for gradual underflow
 
-            while Ex >= Log_Power (Expbits'Last) loop
-               Y := Y * R_Power (Expbits'Last);
-               Ex := Ex - Log_Power (Expbits'Last);
-            end loop;
+            if T'Denorm
+              and then Adjustment >= IEEE_Emin - (Mantissa - 1) - Exp
+            then
+               Expf := IEEE_Emin;
+               Expi := Exp + Adjustment - Expf;
 
-            --  0 <= Ex < 64
+            --  Case of zero result
 
-            for N in reverse Expbits'First .. Expbits'Last - 1 loop
-               if Ex >= Log_Power (N) then
-                  Y := Y * R_Power (N);
-                  Ex := Ex - Log_Power (N);
-               end if;
+            else
+               XX := 0.0;
+               return (if Minus then -XX else XX);
+            end if;
 
-               --  0 <= Ex < Log_Power (N)
+         --  Case of normalized results
 
-            end loop;
+         else
+            Expf := Exp + Adjustment;
+            Expi := 0;
+         end if;
 
-            --  Ex = 0
+         Rep (MSW) := (Rep (MSW) and not Exp_Mask) +
+                        Float_Word (IEEE_Ebias + Expf) * Exp_Factor;
 
+         if Expi < 0 then
+            XX := XX / T (Long_Long_Unsigned (2) ** (-Expi));
          end if;
 
-         return Y;
-      end;
+         return XX;
+      end if;
    end Scaling;
 
    ----------
@@ -653,34 +675,18 @@ package body System.Fat_Gen is
    function Succ (X : T) return T is
       X_Frac : T;
       X_Exp  : UI;
-      X1, X2 : T;
 
    begin
       --  Treat zero specially since it has a zero exponent
 
       if X = 0.0 then
-         X1 := 2.0 ** T'Machine_Emin;
+         return Tiny;
 
-         --  Following loop generates smallest denormal
-
-         loop
-            X2 := T'Machine (X1 / 2.0);
-            exit when X2 = 0.0;
-            X1 := X2;
-         end loop;
-
-         return X1;
-
-      --  Special treatment for largest positive number
+      --  Special treatment for largest positive number: raise Constraint_Error
 
       elsif X = T'Last then
-
-         --  If not generating infinities, we raise a constraint error
-
          raise Constraint_Error with "Succ of largest positive number";
 
-         --  Otherwise generate a positive infinity
-
       --  For infinities, return unchanged
 
       elsif X < T'First or else X > T'Last then
@@ -690,30 +696,35 @@ package body System.Fat_Gen is
          pragma Annotate (CodePeer, Intentional, "dead code",
                           "Check float range.");
 
-      --  Add to the given number a number equivalent to the value
-      --  of its least significant bit. Given that the most significant bit
-      --  represents a value of 1.0 * radix ** (exp - 1), the value we want
-      --  is obtained by shifting this by (mantissa-1) bits to the right,
+      --  Add to the given number a number equivalent to the value of its
+      --  least significant bit. Given that the most significant bit
+      --  represents a value of 1.0 * Radix ** (Exp - 1), the value we want
+      --  is obtained by shifting this by (Mantissa-1) bits to the right,
       --  i.e. decreasing the exponent by that amount.
 
       else
          Decompose (X, X_Frac, X_Exp);
 
-         --  A special case, if the number we had was a negative power of two,
-         --  then we want to add half of what we would otherwise add, since the
-         --  exponent is going to be reduced.
+         --  For a denormalized number or a normalized number with the lowest
+         --  exponent, just add the Tiny.
+
+         if X_Exp <= T'Machine_Emin then
+            return X + Tiny;
 
-         --  Note that X_Frac has the same sign as X, so if X_Frac is -0.5,
-         --  then we know that we have a negative number (and hence a negative
-         --  power of 2).
+         --  A special case, if the number we had was a power of two on the
+         --  negative side of zero, then we want to add half of what we would
+         --  have added, since the exponent is going to be reduced.
 
-         if X_Frac = -0.5 then
-            return X + Gradual_Scaling (X_Exp - T'Machine_Mantissa - 1);
+         --  Note that X_Frac has the same sign as X, so if X_Frac is -Invrad,
+         --  then we know that we had a power of two on the negative side.
 
-         --  Otherwise the exponent is unchanged
+         elsif X_Frac = -Invrad then
+            return X + Scaling (1.0, X_Exp - Mantissa - 1);
+
+         --  Otherwise the adjustment is unchanged
 
          else
-            return X + Gradual_Scaling (X_Exp - T'Machine_Mantissa);
+            return X + Scaling (1.0, X_Exp - Mantissa);
          end if;
       end if;
    end Succ;
@@ -726,7 +737,7 @@ package body System.Fat_Gen is
 
    --    T'Machine (RM1 + N) - RM1
 
-   --  where N >= 0.0 and RM1 = radix ** (mantissa - 1)
+   --  where N >= 0.0 and RM1 = Radix ** (Mantissa - 1)
 
    --  This works provided that the intermediate result (RM1 + N) does not
    --  have extra precision (which is why we call Machine). When we compute
@@ -743,19 +754,18 @@ package body System.Fat_Gen is
    begin
       Result := abs X;
 
-      if Result >= Radix_To_M_Minus_1 then
+      if Result >= RM1 then
          return T'Machine (X);
 
       else
-         Result :=
-           T'Machine (Radix_To_M_Minus_1 + Result) - Radix_To_M_Minus_1;
+         Result := T'Machine (RM1 + Result) - RM1;
 
          if Result > abs X then
             Result := Result - 1.0;
          end if;
 
          if X > 0.0 then
-            return  Result;
+            return Result;
 
          elsif X < 0.0 then
             return -Result;
@@ -806,86 +816,11 @@ package body System.Fat_Gen is
    -----------
 
    function Valid (X : not null access T) return Boolean is
-      IEEE_Emin : constant Integer := T'Machine_Emin - 1;
-      IEEE_Emax : constant Integer := T'Machine_Emax - 1;
-      --  The mantissa is a fraction with first digit set in Ada whereas it is
-      --  shifted by 1 digit to the left in the IEEE floating-point format.
-
-      subtype IEEE_Erange is Integer range IEEE_Emin - 1 .. IEEE_Emax + 1;
-      --  The IEEE floating-point format extends the machine range by 1 to the
-      --  left for denormalized numbers and 1 to the right for infinities/NaNs.
-
-      IEEE_Ebias : constant Integer := -(IEEE_Emin - 1);
-      --  The exponent is biased such that denormalized numbers have it zero
-
-      --  The implementation uses a representation type Float_Rep that allows
-      --  direct access to exponent and mantissa of the floating point number.
-
-      --  The Float_Rep type is a simple array of Float_Word elements. This
-      --  representation is chosen to make it possible to size the type based
-      --  on a generic parameter. Since the array size is known at compile
-      --  time, efficient code can still be generated. The size of Float_Word
-      --  elements should be large enough to allow accessing the exponent in
-      --  one read, but small enough so that all floating-point object sizes
-      --  are a multiple of Float_Word'Size.
-
-      --  The following conditions must be met for all possible instantiations
-      --  of the attribute package:
-
-      --    - T'Size is an integral multiple of Float_Word'Size
-
-      --    - The exponent and sign are completely contained in a single
-      --      component of Float_Rep, named Most Significant Word (MSW).
-
-      --    - The sign occupies the most significant bit of the MSW and the
-      --      exponent is in the following bits. The exception is 80-bit
-      --      double extended, where they occupy the low 16-bit halfword.
-
-      Siz : constant :=
-              (if System.Word_Size > 32 then 32 else System.Word_Size);
-      type Float_Word is mod 2**Siz;
-
-      N : constant Natural := (T'Size + Siz - 1) / Siz;
-      Rep_Last : constant Natural :=
-                   Natural'Min (N - 1, (T'Machine_Mantissa + 16) / Siz);
-      --  Determine the number of Float_Words needed for representing the
-      --  entire floating-point value. Do not take into account excessive
-      --  padding, as occurs on IA-64 where 80 bits floats get padded to 128
-      --  bits. In general, the exponent field cannot be larger than 15 bits,
-      --  even for 128-bit floating-point types, so the final format size
-      --  won't be larger than T'Machine_Mantissa + 16.
-
-      type Float_Rep is array (Natural range 0 .. N - 1) of Float_Word;
-      pragma Suppress_Initialization (Float_Rep);
-      --  This pragma suppresses the generation of an initialization procedure
-      --  for type Float_Rep when operating in Initialize/Normalize_Scalars
-      --  mode, which would be annoying since Valid has got a pragma Inline.
-
-      MSW : constant Natural := Rep_Last * Standard'Default_Bit_Order;
-      --  Finding the location of the Exponent_Word is a bit tricky. In general
-      --  we assume Word_Order = Bit_Order.
-
-      Exp_Factor : constant Float_Word :=
-                     (if T'Machine_Mantissa = 64
-                       then 1
-                       else 2**(Siz - 1) /
-                              Float_Word (IEEE_Emax - IEEE_Emin + 3));
-      --  Factor that the extracted exponent needs to be divided by to be in
-      --  range 0 .. IEEE_Emax - IEEE_Emin + 2. The special case is 80-bit
-      --  double extended, where the exponent starts the 3rd float word.
-
-      Exp_Mask : constant Float_Word :=
-                   Float_Word (IEEE_Emax - IEEE_Emin + 2) * Exp_Factor;
-      --  Value needed to mask out the exponent field. This assumes that the
-      --  range 0 .. IEEE_Emax - IEEE_Emin + 2 contains 2**N values, for some
-      --  N in Natural.
-
       type Access_T is access all T;
       function To_Address is
         new Ada.Unchecked_Conversion (Access_T, System.Address);
 
       Rep : Float_Rep;
-      pragma Import (Ada, Rep);
       for Rep'Address use To_Address (Access_T (X));
       --  Rep is a view of the input floating-point parameter. Note that we
       --  must avoid reading the actual bits of this parameter in float form
diff --git a/gcc/ada/libgnat/s-fatgen.ads b/gcc/ada/libgnat/s-fatgen.ads
index b8bdf2d2..700cfdc 100644
--- a/gcc/ada/libgnat/s-fatgen.ads
+++ b/gcc/ada/libgnat/s-fatgen.ads
@@ -31,9 +31,8 @@
 
 --  This generic package provides a target independent implementation of the
 --  floating-point attributes that denote functions. The implementations here
---  are portable, but very slow. The runtime contains a set of instantiations
---  of this package for all predefined floating-point types, and these should
---  be replaced by efficient assembly language code where possible.
+--  should be portable and reasonably efficient. The runtime contains a set of
+--  instantiations of this package for all predefined floating-point types.
 
 generic
     type T is digits <>;
@@ -107,6 +106,10 @@ package System.Fat_Gen is
    --  floating point register).
 
 private
+   pragma Inline (Compose);
+   pragma Inline (Copy_Sign);
+   pragma Inline (Exponent);
+   pragma Inline (Fraction);
    pragma Inline (Machine);
    pragma Inline (Model);
    pragma Inline (Valid);