path: root/gcc/ada/libgnat/s-fatgen.adb

------------------------------------------------------------------------------
--                                                                          --
--                         GNAT COMPILER COMPONENTS                         --
--                                                                          --
--                       S Y S T E M . F A T _ G E N                        --
--                                                                          --
--                                 B o d y                                  --
--                                                                          --
--          Copyright (C) 1992-2022, 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.      --
--                                                                          --
------------------------------------------------------------------------------

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

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

   Rad : constant T := T (T'Machine_Radix);
   --  Renaming for the machine radix

   Mantissa : constant Integer := T'Machine_Mantissa;
   --  Renaming for the machine mantissa

   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);
   --  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 low-level primitives Copy_Sign, Decompose, Finite_Succ, 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, that is 1 for normalized numbers and 0 for denormalized
   --  numbers, may or may not be stored by the hardware.

   Siz : constant := 16;
   type Float_Word is mod 2**Siz;
   --  We use the GCD of the size of all the supported floating-point formats

   N  : constant Natural := (T'Size + Siz - 1) / Siz;
   NR : constant Natural := (Mantissa + 16 + Siz - 1) / Siz;
   Rep_Last : constant Natural := Natural'Min (N, NR) - 1;
   --  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 :=
                  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

   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 := 2**(Siz - 1);
   --  Value needed to mask out the sign field

   -----------------------
   -- Local Subprograms --
   -----------------------

   procedure Decompose (XX : T; Frac : out T; Expo : out UI);
   --  Decomposes a floating-point number into fraction and exponent parts.
   --  Both results are signed, with Frac having the sign of XX, and UI has
   --  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 Finite_Succ (X : T) return T;
   --  Return the successor of X, a finite number not equal to T'Last

   --------------
   -- Adjacent --
   --------------

   function Adjacent (X, Towards : T) return T is
   begin
      if Towards = X then
         return X;
      elsif Towards > X then
         return Succ (X);
      else
         return Pred (X);
      end if;
   end Adjacent;

   -------------
   -- Ceiling --
   -------------

   function Ceiling (X : T) return T is
      XT : constant T := Truncation (X);
   begin
      if X <= 0.0 then
         return XT;
      elsif X = XT then
         return X;
      else
         return XT + 1.0;
      end if;
   end Ceiling;

   -------------
   -- Compose --
   -------------

   function Compose (Fraction : T; Exponent : UI) return T is
      Arg_Frac : T;
      Arg_Exp  : UI;
   begin
      Decompose (Fraction, Arg_Frac, Arg_Exp);
      return Scaling (Arg_Frac, Exponent);
   end Compose;

   ---------------
   -- Copy_Sign --
   ---------------

   function Copy_Sign (Value, Sign : T) return T is
      S : constant T := T'Machine (Sign);

      Rep_S : Float_Rep;
      for Rep_S'Address use S'Address;
      --  Rep_S is a view of the Sign parameter

      V : T := T'Machine (Value);

      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;

   ---------------
   -- Decompose --
   ---------------

   procedure Decompose (XX : T; Frac : out T; Expo : out UI) is
      X : T := T'Machine (XX);

      Rep : Float_Rep;
      for Rep'Address use X'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

   begin
      --  The normalized exponent of zero is zero, see RM A.5.3(15)

      if X = 0.0 then
         Expo := 0;
         Frac := X;

      --  Check for infinities and NaNs

      elsif Exp = IEEE_Emax + 1 then
         Expo := T'Machine_Emax + 1;
         Frac := (if Minus then -Invrad else Invrad);

      --  Check for nonzero denormalized numbers

      elsif Exp = IEEE_Emin - 1 then
         --  Normalize by multiplying by Radix ** (Mantissa - 1)

         Decompose (X * RM1, Frac, Expo);
         Expo := Expo - (Mantissa - 1);

      --  Case of normalized numbers

      else
         --  The Ada exponent is the IEEE exponent plus 1, see above

         Expo := Exp + 1;

         --  Set Ada exponent of X to zero, so we end up with the fraction

         Rep (MSW) := (Rep (MSW) and not Exp_Mask) +
                        Float_Word (IEEE_Ebias - 1) * Exp_Factor;
         Frac := X;
      end if;
   end Decompose;

   --------------
   -- Exponent --
   --------------

   function Exponent (X : T) return UI is
      X_Frac : T;
      X_Exp  : UI;
   begin
      Decompose (X, X_Frac, X_Exp);
      return X_Exp;
   end Exponent;

   -----------------
   -- Finite_Succ --
   -----------------

   function Finite_Succ (X : T) return T is
      XX : T := T'Machine (X);

      Rep : Float_Rep;
      for Rep'Address use XX'Address;
      --  Rep is a view of the input floating-point parameter

   begin
      --  If the floating-point type does not support denormalized numbers,
      --  there is a couple of problematic values, namely -Small and Zero,
      --  because the increment is equal to Small in these cases.

      if not T'Denorm then
         declare
            Small : constant T := Rad ** (T'Machine_Emin - 1);
            --  Smallest positive normalized number declared here and not at
            --  library level for the sake of the CCG compiler, which cannot
            --  currently compile the constant because the target is C90.

         begin
            if X = -Small then
               XX := 0.0;
               return -XX;
            elsif X = 0.0 then
               return Small;
            end if;
         end;
      end if;

      --  In all the other cases, the increment is equal to 1 in the binary
      --  integer representation of the number if X is nonnegative and equal
      --  to -1 if X is negative.

      if XX >= 0.0 then
         --  First clear the sign of negative Zero

         Rep (MSW) := Rep (MSW) and not Sign_Mask;

         --  Deal with big endian

         if MSW = 0 then
            for J in reverse 0 .. Rep_Last loop
               Rep (J) := Rep (J) + 1;

               --  For 80-bit IEEE Extended, the MSB of the mantissa is stored
               --  so, when it has been flipped, its status must be reanalyzed.

               if Mantissa = 64 and then J = 1 then

                  --  If the MSB changed from denormalized to normalized, then
                  --  keep it normalized since the exponent will be bumped.

                  if Rep (J) = 2**(Siz - 1) then
                     null;

                  --  If the MSB changed from normalized, restore it since we
                  --  cannot denormalize in this context.

                  elsif Rep (J) = 0 then
                     Rep (J) := 2**(Siz - 1);

                  else
                     exit;
                  end if;

               --  In other cases, stop if there is no carry

               else
                  exit when Rep (J) > 0;
               end if;
            end loop;

         --  Deal with little endian

         else
            for J in 0 .. Rep_Last loop
               Rep (J) := Rep (J) + 1;

               --  For 80-bit IEEE Extended, the MSB of the mantissa is stored
               --  so, when it has been flipped, its status must be reanalyzed.

               if Mantissa = 64 and then J = Rep_Last - 1 then

                  --  If the MSB changed from denormalized to normalized, then
                  --  keep it normalized since the exponent will be bumped.

                  if Rep (J) = 2**(Siz - 1) then
                     null;

                  --  If the MSB changed from normalized, restore it since we
                  --  cannot denormalize in this context.

                  elsif Rep (J) = 0 then
                     Rep (J) := 2**(Siz - 1);

                  else
                     exit;
                  end if;

               --  In other cases, stop if there is no carry

               else
                  exit when Rep (J) > 0;
               end if;
            end loop;
         end if;

      else
         if MSW = 0 then
            for J in reverse 0 .. Rep_Last loop
               Rep (J) := Rep (J) - 1;

               --  For 80-bit IEEE Extended, the MSB of the mantissa is stored
               --  so, when it has been flipped, its status must be reanalyzed.

               if Mantissa = 64 and then J = 1 then

                  --  If the MSB changed from normalized to denormalized, then
                  --  keep it normalized if the exponent is not 1.

                  if Rep (J) = 2**(Siz - 1) - 1 then
                     if Rep (0) /= 2**(Siz - 1) + 1 then
                        Rep (J) := 2**Siz - 1;
                     end if;

                  else
                     exit;
                  end if;

               --  In other cases, stop if there is no borrow

               else
                  exit when Rep (J) < 2**Siz - 1;
               end if;
            end loop;

         else
            for J in 0 .. Rep_Last loop
               Rep (J) := Rep (J) - 1;

               --  For 80-bit IEEE Extended, the MSB of the mantissa is stored
               --  so, when it has been flipped, its status must be reanalyzed.

               if Mantissa = 64 and then J = Rep_Last - 1 then

                  --  If the MSB changed from normalized to denormalized, then
                  --  keep it normalized if the exponent is not 1.

                  if Rep (J) = 2**(Siz - 1) - 1 then
                     if Rep (Rep_Last) /= 2**(Siz - 1) + 1 then
                        Rep (J) := 2**Siz - 1;
                     end if;

                  else
                     exit;
                  end if;

               --  In other cases, stop if there is no borrow

               else
                  exit when Rep (J) < 2**Siz - 1;
               end if;
            end loop;
         end if;
      end if;

      return XX;
   end Finite_Succ;

   -----------
   -- Floor --
   -----------

   function Floor (X : T) return T is
      XT : constant T := Truncation (X);
   begin
      if X >= 0.0 then
         return XT;
      elsif XT = X then
         return X;
      else
         return XT - 1.0;
      end if;
   end Floor;

   --------------
   -- Fraction --
   --------------

   function Fraction (X : T) return T is
      X_Frac : T;
      X_Exp  : UI;
   begin
      Decompose (X, X_Frac, X_Exp);
      return X_Frac;
   end Fraction;

   ------------------
   -- Leading_Part --
   ------------------

   function Leading_Part (X : T; Radix_Digits : UI) return T is
      L    : UI;
      Y, Z : T;

   begin
      if Radix_Digits >= Mantissa then
         return X;

      elsif Radix_Digits <= 0 then
         raise Constraint_Error;

      else
         L := Exponent (X) - Radix_Digits;
         Y := Truncation (Scaling (X, -L));
         Z := Scaling (Y, L);
         return Z;
      end if;
   end Leading_Part;

   -------------
   -- Machine --
   -------------

   --  The trick with Machine is to force the compiler to store the result
   --  in memory so that we do not have extra precision used. The compiler
   --  is clever, so we have to outwit its possible optimizations. We do
   --  this by using an intermediate pragma Volatile location.

   function Machine (X : T) return T is
      Temp : T;
      pragma Volatile (Temp);
   begin
      Temp := X;
      return Temp;
   end Machine;

   ----------------------
   -- Machine_Rounding --
   ----------------------

   --  For now, the implementation is identical to that of Rounding, which is
   --  a permissible behavior, but is not the most efficient possible approach.

   function Machine_Rounding (X : T) return T is
      Result : T;
      Tail   : T;

   begin
      Result := Truncation (abs X);
      Tail   := abs X - Result;

      if Tail >= 0.5 then
         Result := Result + 1.0;
      end if;

      if X > 0.0 then
         return Result;

      elsif X < 0.0 then
         return -Result;

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

      else
         return X;
      end if;
   end Machine_Rounding;

   -----------
   -- Model --
   -----------

   --  We treat Model as identical to Machine. This is true of IEEE and other
   --  nice floating-point systems, but not necessarily true of all systems.

   function Model (X : T) return T is
   begin
      return T'Machine (X);
   end Model;

   ----------
   -- Pred --
   ----------

   function Pred (X : T) return T is
   begin
      --  Special treatment for largest negative number: raise Constraint_Error

      if X = T'First then
         raise Constraint_Error with "Pred of largest negative number";

      --  For finite numbers, use the symmetry around zero of floating point

      elsif X > T'First and then X <= T'Last then
         pragma Annotate (CodePeer, Intentional, "test always true",
                          "Check for invalid float");
         pragma Annotate (CodePeer, Intentional, "condition predetermined",
                          "Check for invalid float");
         return -Finite_Succ (-X);

      --  For infinities and NaNs, return unchanged

      else
         return X;
         pragma Annotate (CodePeer, Intentional, "dead code",
                          "Check float range.");
      end if;
   end Pred;

   ---------------
   -- Remainder --
   ---------------

   function Remainder (X, Y : T) return T is
      A        : T;
      B        : T;
      Arg      : T;
      P        : T;
      P_Frac   : T;
      Sign_X   : T;
      IEEE_Rem : T;
      Arg_Exp  : UI;
      P_Exp    : UI;
      K        : UI;
      P_Even   : Boolean;

      Arg_Frac : T;

   begin
      if Y = 0.0 then
         raise Constraint_Error;
      end if;

      if X > 0.0 then
         Sign_X :=  1.0;
         Arg := X;
      else
         Sign_X := -1.0;
         Arg := -X;
      end if;

      P := abs Y;

      if Arg < P then
         P_Even := True;
         IEEE_Rem := Arg;
         P_Exp := Exponent (P);

      else
         Decompose (Arg, Arg_Frac, Arg_Exp);
         Decompose (P,   P_Frac,   P_Exp);

         P := Compose (P_Frac, Arg_Exp);
         K := Arg_Exp - P_Exp;
         P_Even := True;
         IEEE_Rem := Arg;

         for Cnt in reverse 0 .. K loop
            if IEEE_Rem >= P then
               P_Even := False;
               IEEE_Rem := IEEE_Rem - P;
            else
               P_Even := True;
            end if;

            P := P * 0.5;
         end loop;
      end if;

      --  That completes the calculation of modulus remainder. The final
      --  step is get the IEEE remainder. Here we need to compare Rem with
      --  (abs Y) / 2. We must be careful of unrepresentable Y/2 value
      --  caused by subnormal numbers

      if P_Exp >= 0 then
         A := IEEE_Rem;
         B := abs Y * 0.5;

      else
         A := IEEE_Rem * 2.0;
         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 (X : T) return T is
      Result : T;
      Tail   : T;

   begin
      Result := Truncation (abs X);
      Tail   := abs X - Result;

      if Tail >= 0.5 then
         Result := Result + 1.0;
      end if;

      if X > 0.0 then
         return Result;

      elsif X < 0.0 then
         return -Result;

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

      else
         return X;
      end if;
   end Rounding;

   -------------
   -- Scaling --
   -------------

   function Scaling (X : T; Adjustment : UI) return T is
      pragma Assert (Mantissa <= 64);
      --  This implementation handles only 80-bit IEEE Extended or smaller

      package UST renames System.Unsigned_Types;
      use type UST.Long_Long_Unsigned;

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

      --  Check for nonzero denormalized numbers

      elsif Exp = IEEE_Emin - 1 then
         --  Check for zero result to protect the subtraction below

         if Adjustment < -(Mantissa - 1) then
            XX := 0.0;
            return (if Minus then -XX else XX);

         --  Normalize by multiplying by Radix ** (Mantissa - 1)

         else
            return Scaling (XX * RM1, Adjustment - (Mantissa - 1));
         end if;

      --  Case of normalized numbers

      else
         --  Check for overflow

         if Adjustment > IEEE_Emax - Exp then
            --  Optionally raise Constraint_Error as per RM A.5.3(29)

            if T'Machine_Overflows then
               raise Constraint_Error with "Too large exponent";

            else
               XX := 0.0;
               return (if Minus then -1.0 / XX else 1.0 / XX);
               pragma Annotate (CodePeer, Intentional, "overflow check",
                                "Infinity produced");
               pragma Annotate (CodePeer, Intentional, "divide by zero",
                                "Infinity produced");
            end if;

         --  Check for underflow

         elsif Adjustment < IEEE_Emin - Exp then
            --  Check for possibly gradual underflow (up to the hardware)

            if Adjustment >= IEEE_Emin - Mantissa - Exp then
               Expf := IEEE_Emin;
               Expi := Exp + Adjustment - Expf;

            --  Case of zero result

            else
               XX := 0.0;
               return (if Minus then -XX else XX);
            end if;

         --  Case of normalized results

         else
            Expf := Exp + Adjustment;
            Expi := 0;
         end if;

         Rep (MSW) := (Rep (MSW) and not Exp_Mask) +
                        Float_Word (IEEE_Ebias + Expf) * Exp_Factor;

         if Expi < 0 then
            --  Given that Expi >= -Mantissa, only -64 is problematic

            if Expi = -64 then
               pragma Annotate
                 (CodePeer, Intentional, "test always false",
                  "test always false in some instantiations");
               XX := XX / 2.0;
               Expi := -63;
            end if;

            XX := XX / T (UST.Long_Long_Unsigned (2) ** (-Expi));
         end if;

         return XX;
      end if;
   end Scaling;

   ----------
   -- Succ --
   ----------

   function Succ (X : T) return T is
   begin
      --  Special treatment for largest positive number: raise Constraint_Error

      if X = T'Last then
         raise Constraint_Error with "Succ of largest positive number";

      --  For finite numbers, call the specific routine

      elsif X >= T'First and then X < T'Last then
         pragma Annotate (CodePeer, Intentional, "test always true",
                          "Check for invalid float");
         pragma Annotate (CodePeer, Intentional, "condition predetermined",
                          "Check for invalid float");
         return Finite_Succ (X);

      --  For infinities and NaNs, return unchanged

      else
         return X;
         pragma Annotate (CodePeer, Intentional, "dead code",
                          "Check float range.");
      end if;
   end Succ;

   ----------------
   -- Truncation --
   ----------------

   --  The basic approach is to compute

   --    T'Machine (RM1 + N) - RM1

   --  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
   --  RM1 + N, the exponent of N will be normalized and the mantissa shifted
   --  appropriately so the lower order bits, which cannot contribute to the
   --  integer part of N, fall off on the right. When we subtract RM1 again,
   --  the significant bits of N are shifted to the left, and what we have is
   --  an integer, because only the first e bits are different from zero
   --  (assuming binary radix here).

   function Truncation (X : T) return T is
      Result : T;

   begin
      Result := abs X;

      if Result >= RM1 then
         return T'Machine (X);

      else
         Result := T'Machine (RM1 + Result) - RM1;

         if Result > abs X then
            Result := Result - 1.0;
         end if;

         if X > 0.0 then
            return Result;

         elsif X < 0.0 then
            return -Result;

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

         else
            return X;
         end if;
      end if;
   end Truncation;

   -----------------------
   -- Unbiased_Rounding --
   -----------------------

   function Unbiased_Rounding (X : T) return T is
      Abs_X  : constant T := abs X;
      Result : T;
      Tail   : T;

   begin
      Result := Truncation (Abs_X);
      Tail   := Abs_X - Result;

      if Tail > 0.5 then
         Result := Result + 1.0;

      elsif Tail = 0.5 then
         Result := 2.0 * Truncation ((Result / 2.0) + 0.5);
      end if;

      if X > 0.0 then
         return Result;

      elsif X < 0.0 then
         return -Result;

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

      else
         return X;
      end if;
   end Unbiased_Rounding;

   -----------
   -- Valid --
   -----------

   function Valid (X : not null access T) return Boolean is
      type Access_T is access all T;
      function To_Address is
        new Ada.Unchecked_Conversion (Access_T, System.Address);

      Rep : Float_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
      --  since it may be a signalling NaN.

      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.

   begin
      if Exp = IEEE_Emax + 1 then
         --  This is an infinity or a NaN, i.e. always invalid

         return False;

      elsif Exp in IEEE_Emin .. IEEE_Emax then
         --  This is a normalized number, i.e. always valid

         return True;

      else pragma Assert (Exp = IEEE_Emin - 1);
         --  This is a denormalized number, valid if T'Denorm is True or 0.0

         if T'Denorm then
            return True;

         --  Note that we cannot do a direct comparison with 0.0 because the
         --  hardware may evaluate it to True for all denormalized numbers.

         else
            --  First clear the sign bit (the exponent is already zero)

            Rep (MSW) := Rep (MSW) and not Sign_Mask;

            return (for all J in 0 .. Rep_Last => Rep (J) = 0);
         end if;
      end if;
   end Valid;

end System.Fat_Gen;