Geek Gadgets 1/ADE-1.bin/ade-dist/gnat-2.06-src.tgz/tar.out/fsf/gnat/ada/s-fatgen.adb

Text File | 1996-09-28 | 17KB | 659 lines

------------------------------------------------------------------------------
--                                                                          --
--                         GNAT COMPILER COMPONENTS                         --
--                                                                          --
--                       S Y S T E M . F A T _ G E N                        --
--                                                                          --
--                                 B o d y                                  --
--                                                                          --
--                            $Revision: 1.5 $                              --
--                                                                          --
-- The GNAT library is free software; you can redistribute it and/or modify --
-- it under terms of the GNU Library General Public License as published by --
-- the Free Software  Foundation; either version 2, or (at your option) any --
-- later version.  The GNAT library is distributed in the hope that it will --
-- be useful, but WITHOUT ANY WARRANTY;  without even  the implied warranty --
-- of MERCHANTABILITY  or  FITNESS FOR  A PARTICULAR PURPOSE.  See the  GNU --
-- Library  General  Public  License for  more  details.  You  should  have --
-- received  a copy of the GNU  Library  General Public License  along with --
-- the GNAT library;  see the file  COPYING.LIB.  If not, write to the Free --
-- Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.        --
--                                                                          --
------------------------------------------------------------------------------

--  The implementation here is portable to any IEEE implementation. It does
--  not handle non-binary 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.

package body System.Fat_Gen is

   Float_Radix        : constant T := T (T'Machine_Radix);
   Float_Radix_Inv    : constant T := 1.0 / Float_Radix;
   Radix_To_M_Minus_1 : constant T := Float_Radix ** (T'Machine_Mantissa - 1);

   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);
   Invrad : constant T := 1.0 / Rad;

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

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

   procedure Decompose (XX : T; Frac : out T; Expo : out UI);
   --  Decomposes a floating-point number into fraction and exponent parts

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

   --  Adjacent( X, 0.0) will be incorrect for abs(X) any power of two
   --  because of a bug in Succ and Pred.  An extra test must be made
   --  for this case ???.

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

   --  Configuration requirement: If this unit is used to implement the
   --  floating-point attributes, then if Signed_Zeros is true, we must
   --  have Machine_Overflows False, and we must generate infinities for
   --  overflow (otherwise we can't implement signed zeroes properly).

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

   begin
      Result := abs Value;

      --  If we have signed zeroes true, then we know that machine overflows
      --  is false, and that infinities are generated for overflow, so we can
      --  test for minus zero by looking at the sign of the corresponding
      --  infinity value (neat trick eh?)

      if Sign = 0.0 and then T'Signed_Zeros then
         if 1.0 / Sign > 0.0 then
            return Result;
         else
            return -Result;
         end if;
      end if;

      --  Handle non-zero cases, and also the zero case if signed zeroes
      --  is false. In the latter case we always treat zero as positive.

      if Sign >= 0.0 then
         return Result;
      else
         return -Result;
      end if;
   end Copy_Sign;

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

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

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

         --  More useful would be defining Expo to be T'Machine_Emin - 1 or
         --  T'Machine_Emin - T'Machine_Mantissa, which would preserve
         --  monotonicity of the exponent fuction..

      --  Check for infinities, transfinites, whatnot.

      elsif X > T'Safe_Last then
         Frac := Invrad;
         Expo := T'Machine_Emax + 1;

      elsif X < T'Safe_First then
         Frac := -Invrad;
         Expo := T'Machine_Emax + 2;    -- how many extra negative values?

      else
         --  Case of nonzero finite x. Essentially, we just multiply
         --  by Rad ** (+-2**N) to reduce the range.

         declare
            Ax : T  := abs X;
            Ex : UI := 0;

         --  Ax * Rad ** Ex is invariant.

         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;

               --  Ax < Rad ** 64

               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;

                  --  Ax < R_Power (N)
               end loop;

               --  1 <= Ax < Rad

               Ax := Ax * Invrad;
               Ex := Ex + 1;

            else
               --  0 < ax < 1

               while Ax < R_Neg_Power (Expbits'Last) loop
                  Ax := Ax * R_Power (Expbits'Last);
                  Ex := Ex - Log_Power (Expbits'Last);
               end loop;

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

            if X > 0.0 then
               Frac := Ax;
            else
               Frac := -Ax;
            end if;

            Expo := Ex;
         end;
      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;

   -----------
   -- 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 >= T'Machine_Mantissa then
         return X;

      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!

   --  This is achieved by using the following array. In fact only the first
   --  element is used, and Machine_Index is always 1, but we make sure the
   --  compiler can't figure this out.

   Machine_Array : array (1 .. 2) of T;
   Machine_Index : Integer;

   function Machine (X : T) return T is
   begin
      Machine_Array (1) := X;
      return Machine_Array (Machine_Index);
   end Machine;

   -----------
   -- 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 Machine (X);
   end Model;

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

   --  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, i.e. decreasing the
   --  exponent by that amount.

   function Pred (X : T) return T is
      X_Frac : T;
      X_Exp  : UI;

   begin
      Decompose (X, X_Frac, X_Exp);
      return X - Scaling (1.0, X_Exp - T'Machine_Mantissa);
   end Pred;

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

   function Remainder (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;

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

      Arg := abs X;
      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 := UI (Arg_Exp) - UI (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 --
   -------------

   --  Return x * rad ** adjustment quickly,
   --  or quietly underflow to zero, or overflow naturally.

   function Scaling (X : T; Adjustment : UI) return T is
   begin
      if X = 0.0 or else Adjustment = 0 then
         return X;
      end if;

      --  Nonzero x. essentially, just multiply repeatedly by Rad ** (+-2**n).

      declare
         Y  : T  := X;
         Ex : UI := Adjustment;

      --  Y * Rad ** Ex is invariant

      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;

            --  -64 < Ex <= 0

            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;

               --  -Log_Power (N) < Ex <= 0
            end loop;

            --  Ex = 0

         else
            --  Ex >= 0

            while Ex >= Log_Power (Expbits'Last) loop
               Y := Y * R_Power (Expbits'Last);
               Ex := Ex - Log_Power (Expbits'Last);
            end loop;

            --  0 <= Ex < 64

            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;

               --  0 <= Ex < Log_Power (N)
            end loop;

            --  Ex = 0
         end if;
         return Y;
      end;
   end Scaling;

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

   --  Similar computation to that of Pred: find value of least significant
   --  bit of given number, and add.

   function Succ (X : T) return T is
      X_Frac : T;
      X_Exp  : UI;

   begin
      Decompose (X, X_Frac, X_Exp);
      return X + Scaling (1.0, X_Exp - T'Machine_Mantissa);
   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 epxonent of N will be normalized and the mantissa shifted
   --  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 >= Radix_To_M_Minus_1 then
         return Machine (X);

      else
         Result := Machine (Radix_To_M_Minus_1 + Result) - Radix_To_M_Minus_1;

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

   ----------------------------
   -- Package Initialization --
   ----------------------------

begin
   --  Set Machine_Index to 1, but not in an obvious manner (see function
   --  Machine to understand why we are behaving in this secretive manner)

   Machine_Index := 2 ** 3;

   for J in 1 .. 3 loop
      Machine_Index := Machine_Index / 2;
   end loop;

end System.Fat_Gen;