[Ada] Implement tiered support for floating-point exponentiation

Message ID 20210506075812.GA125513@adacore.com
State New
Headers show
Series
  • [Ada] Implement tiered support for floating-point exponentiation
Related show

Commit Message

Pierre-Marie de Rodat May 6, 2021, 7:58 a.m.
This changes the implementation of floating-point exponentiation from
using Long_Long_Float for all floating-point types to using a base type
tailored to the type being operated on.

The computation is done in double precision internally, which gives more
accurate results for Long_Float and Long_Long_Float.

Tested on x86_64-pc-linux-gnu, committed on trunk

gcc/ada/

	* Makefile.rtl (GNATRTL_NONTASKING_OBJS): Add s-exponr, s-exnflt
	and s-exnlfl.
	* exp_ch4.adb (Expand_N_Op_Expon): Use RE_Exn_Float for Short_Float.
	* rtsfind.ads (RTU_Id): Add System_Exn_Flt and System_Exn_LFlt.
	(RE_Id): Adjust entries for RE_Exn_Float and RE_Exn_Long_Float.
	(RE_Unit_Table): Likewise.
	* libgnat/s-exnflt.ads: New file.
	* libgnat/s-exnlfl.ads: Likewise.
	* libgnat/s-exnllf.ads: Change to mere instantiation.
	* libgnat/s-exnllf.adb: Move implementation to...
	* libgnat/s-exponr.ads: New generic unit.
	* libgnat/s-exponr.adb: ...here and also make it generic.
	(Expon): Do the computation in double precision internally.

Patch

diff --git a/gcc/ada/Makefile.rtl b/gcc/ada/Makefile.rtl
--- a/gcc/ada/Makefile.rtl
+++ b/gcc/ada/Makefile.rtl
@@ -583,6 +583,8 @@  GNATRTL_NONTASKING_OBJS= \
   s-exctab$(objext) \
   s-exctra$(objext) \
   s-exnint$(objext) \
+  s-exnflt$(objext) \
+  s-exnlfl$(objext) \
   s-exnllf$(objext) \
   s-exnlli$(objext) \
   s-expint$(objext) \
@@ -590,6 +592,7 @@  GNATRTL_NONTASKING_OBJS= \
   s-expllu$(objext) \
   s-expmod$(objext) \
   s-exponn$(objext) \
+  s-exponr$(objext) \
   s-expont$(objext) \
   s-exponu$(objext) \
   s-expuns$(objext) \


diff --git a/gcc/ada/exp_ch4.adb b/gcc/ada/exp_ch4.adb
--- a/gcc/ada/exp_ch4.adb
+++ b/gcc/ada/exp_ch4.adb
@@ -9081,15 +9081,12 @@  package body Exp_Ch4 is
       --  overflow), and if there is an infinity generated and a range check
       --  is required, the check will fail anyway.
 
-      --  Historical note: we used to convert everything to Long_Long_Float
-      --  and call a single common routine, but this had the undesirable effect
-      --  of giving different results for small static exponent values and the
-      --  same dynamic values.
-
       else
          pragma Assert (Is_Floating_Point_Type (Rtyp));
 
-         if Rtyp = Standard_Float then
+         --  Short_Float and Float are the same type for GNAT
+
+         if Rtyp = Standard_Short_Float or else Rtyp = Standard_Float then
             Etyp := Standard_Float;
             Rent := RE_Exn_Float;
 


diff --git /dev/null b/gcc/ada/libgnat/s-exnflt.ads
new file mode 100644
--- /dev/null
+++ b/gcc/ada/libgnat/s-exnflt.ads
@@ -0,0 +1,41 @@ 
+------------------------------------------------------------------------------
+--                                                                          --
+--                         GNAT RUN-TIME COMPONENTS                         --
+--                                                                          --
+--                       S Y S T E M . E X N _ F L T                        --
+--                                                                          --
+--                                 S p e c                                  --
+--                                                                          --
+--            Copyright (C) 2021, 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.      --
+--                                                                          --
+------------------------------------------------------------------------------
+
+--  Float exponentiation (checks off)
+
+with System.Exponr;
+
+package System.Exn_Flt is
+
+   function Exn_Float is new Exponr (Float);
+   pragma Pure_Function (Exn_Float);
+
+end System.Exn_Flt;


diff --git /dev/null b/gcc/ada/libgnat/s-exnlfl.ads
new file mode 100644
--- /dev/null
+++ b/gcc/ada/libgnat/s-exnlfl.ads
@@ -0,0 +1,41 @@ 
+------------------------------------------------------------------------------
+--                                                                          --
+--                         GNAT RUN-TIME COMPONENTS                         --
+--                                                                          --
+--                      S Y S T E M . E X N _ L F L T                       --
+--                                                                          --
+--                                 S p e c                                  --
+--                                                                          --
+--            Copyright (C) 2021, 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.      --
+--                                                                          --
+------------------------------------------------------------------------------
+
+--  Long_Float exponentiation (checks off)
+
+with System.Exponr;
+
+package System.Exn_LFlt is
+
+   function Exn_Long_Float is new Exponr (Long_Float);
+   pragma Pure_Function (Exn_Long_Float);
+
+end System.Exn_LFlt;


diff --git a/gcc/ada/libgnat/s-exnllf.adb b/gcc/ada/libgnat/s-exnllf.adb
--- a/gcc/ada/libgnat/s-exnllf.adb
+++ b/gcc/ada/libgnat/s-exnllf.adb
@@ -29,154 +29,8 @@ 
 --                                                                          --
 ------------------------------------------------------------------------------
 
---  Note: the reason for treating exponents in the range 0 .. 4 specially is
---  to ensure identical results to the static inline expansion in the case of
---  a compile time known exponent in this range. The use of Float'Machine and
---  Long_Float'Machine is to avoid unwanted extra precision in the results.
+--  This package does not require a body, since it is an instantiation. We
+--  provide a dummy file containing a No_Body pragma so that previous versions
+--  of the body (which did exist) will not interfere.
 
---  Note that for a negative exponent in Left ** Right, we compute the result
---  as:
-
---     1.0 / (Left ** (-Right))
-
---  Note that the case of Left being zero is not special, it will simply result
---  in a division by zero at the end, yielding a correctly signed infinity, or
---  possibly generating an overflow.
-
---  Note on overflow: This coding assumes that the target generates infinities
---  with standard IEEE semantics. If this is not the case, then the code
---  for negative exponent may raise Constraint_Error. This follows the
---  implementation permission given in RM 4.5.6(12).
-
-package body System.Exn_LLF is
-
-   subtype Negative is Integer range Integer'First .. -1;
-
-   function Exp
-     (Left  : Long_Long_Float;
-      Right : Natural) return Long_Long_Float;
-   --  Common routine used if Right is greater or equal to 5
-
-   ---------------
-   -- Exn_Float --
-   ---------------
-
-   function Exn_Float
-     (Left  : Float;
-      Right : Integer) return Float
-   is
-      Temp : Float;
-   begin
-      case Right is
-         when 0 =>
-            return 1.0;
-         when 1 =>
-            return Left;
-         when 2 =>
-            return Float'Machine (Left * Left);
-         when 3 =>
-            return Float'Machine (Left * Left * Left);
-         when 4 =>
-            Temp := Float'Machine (Left * Left);
-            return Float'Machine (Temp * Temp);
-         when Negative =>
-            return Float'Machine (1.0 / Exn_Float (Left, -Right));
-         when others =>
-            return
-              Float'Machine
-                (Float (Exp (Long_Long_Float (Left), Right)));
-      end case;
-   end Exn_Float;
-
-   --------------------
-   -- Exn_Long_Float --
-   --------------------
-
-   function Exn_Long_Float
-     (Left  : Long_Float;
-      Right : Integer) return Long_Float
-   is
-      Temp : Long_Float;
-   begin
-      case Right is
-         when 0 =>
-            return 1.0;
-         when 1 =>
-            return Left;
-         when 2 =>
-            return Long_Float'Machine (Left * Left);
-         when 3 =>
-            return Long_Float'Machine (Left * Left * Left);
-         when 4 =>
-            Temp := Long_Float'Machine (Left * Left);
-            return Long_Float'Machine (Temp * Temp);
-         when Negative =>
-            return Long_Float'Machine (1.0 / Exn_Long_Float (Left, -Right));
-         when others =>
-            return
-              Long_Float'Machine
-                (Long_Float (Exp (Long_Long_Float (Left), Right)));
-      end case;
-   end Exn_Long_Float;
-
-   -------------------------
-   -- Exn_Long_Long_Float --
-   -------------------------
-
-   function Exn_Long_Long_Float
-     (Left  : Long_Long_Float;
-      Right : Integer) return Long_Long_Float
-   is
-      Temp : Long_Long_Float;
-   begin
-      case Right is
-         when 0 =>
-            return 1.0;
-         when 1 =>
-            return Left;
-         when 2 =>
-            return Left * Left;
-         when 3 =>
-            return Left * Left * Left;
-         when 4 =>
-            Temp := Left * Left;
-            return Temp * Temp;
-         when Negative =>
-            return 1.0 / Exn_Long_Long_Float (Left, -Right);
-         when others =>
-            return Exp (Left, Right);
-      end case;
-   end Exn_Long_Long_Float;
-
-   ---------
-   -- Exp --
-   ---------
-
-   function Exp
-     (Left  : Long_Long_Float;
-      Right : Natural) return Long_Long_Float
-   is
-      Result : Long_Long_Float := 1.0;
-      Factor : Long_Long_Float := Left;
-      Exp    : Natural := Right;
-
-   begin
-      --  We use the standard logarithmic approach, Exp gets shifted right
-      --  testing successive low order bits and Factor is the value of the
-      --  base raised to the next power of 2. If the low order bit or Exp is
-      --  set, multiply the result by this factor.
-
-      loop
-         if Exp rem 2 /= 0 then
-            Result := Result * Factor;
-         end if;
-
-         Exp := Exp / 2;
-         exit when Exp = 0;
-         Factor := Factor * Factor;
-      end loop;
-
-      return Result;
-   end Exp;
-
-end System.Exn_LLF;
+pragma No_Body;


diff --git a/gcc/ada/libgnat/s-exnllf.ads b/gcc/ada/libgnat/s-exnllf.ads
--- a/gcc/ada/libgnat/s-exnllf.ads
+++ b/gcc/ada/libgnat/s-exnllf.ads
@@ -29,21 +29,13 @@ 
 --                                                                          --
 ------------------------------------------------------------------------------
 
---  [Long_[Long_]]Float exponentiation (checks off)
+--  Long_Long_Float exponentiation (checks off)
 
-package System.Exn_LLF is
-   pragma Pure;
-
-   function Exn_Float
-     (Left  : Float;
-      Right : Integer) return Float;
+with System.Exponr;
 
-   function Exn_Long_Float
-     (Left  : Long_Float;
-      Right : Integer) return Long_Float;
+package System.Exn_LLF is
 
-   function Exn_Long_Long_Float
-     (Left  : Long_Long_Float;
-      Right : Integer) return Long_Long_Float;
+   function Exn_Long_Long_Float is new Exponr (Long_Long_Float);
+   pragma Pure_Function (Exn_Long_Long_Float);
 
 end System.Exn_LLF;


diff --git /dev/null b/gcc/ada/libgnat/s-exponr.adb
new file mode 100644
--- /dev/null
+++ b/gcc/ada/libgnat/s-exponr.adb
@@ -0,0 +1,122 @@ 
+------------------------------------------------------------------------------
+--                                                                          --
+--                         GNAT RUN-TIME COMPONENTS                         --
+--                                                                          --
+--                        S Y S T E M . E X P O N R                         --
+--                                                                          --
+--                                 B o d y                                  --
+--                                                                          --
+--            Copyright (C) 2021, 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.      --
+--                                                                          --
+------------------------------------------------------------------------------
+
+--  Note that the reason for treating exponents in the range 0 .. 4 specially
+--  is to ensure identical results with the static expansion in the case of a
+--  compile-time known exponent in this range; similarly, the use 'Machine is
+--  to avoid unwanted extra precision in the results.
+
+--  For a negative exponent, we compute the result as per RM 4.5.6(11/3):
+
+--     Left ** Right = 1.0 / (Left ** (-Right))
+
+--  Note that the case of Left being zero is not special, it will simply result
+--  in a division by zero at the end, yielding a correctly signed infinity, or
+--  possibly raising an overflow exception.
+
+--  Note on overflow: this coding assumes that the target generates infinities
+--  with standard IEEE semantics. If this is not the case, then the code for
+--  negative exponents may raise Constraint_Error, which is in keeping with the
+--  implementation permission given in RM 4.5.6(12).
+
+with System.Double_Real;
+
+function System.Exponr (Left : Num; Right : Integer) return Num is
+
+   package Double_Real is new System.Double_Real (Num);
+   use type Double_Real.Double_T;
+
+   subtype Double_T is Double_Real.Double_T;
+   --  The double floating-point type
+
+   subtype Negative is Integer range Integer'First .. -1;
+   --  The range of negative exponents
+
+   function Expon (Left : Num; Right : Natural) return Num;
+   --  Routine used if Right is greater than 4
+
+   -----------
+   -- Expon --
+   -----------
+
+   function Expon (Left : Num; Right : Natural) return Num is
+      Result : Double_T := Double_Real.To_Double (1.0);
+      Factor : Double_T := Double_Real.To_Double (Left);
+      Exp    : Natural  := Right;
+
+   begin
+      --  We use the standard logarithmic approach, Exp gets shifted right
+      --  testing successive low order bits and Factor is the value of the
+      --  base raised to the next power of 2. If the low order bit or Exp
+      --  is set, multiply the result by this factor.
+
+      loop
+         if Exp rem 2 /= 0 then
+            Result := Result * Factor;
+            exit when Exp = 1;
+         end if;
+
+         Exp := Exp / 2;
+         Factor := Double_Real.Sqr (Factor);
+      end loop;
+
+      return Double_Real.To_Single (Result);
+   end Expon;
+
+begin
+   case Right is
+      when 0 =>
+         return 1.0;
+
+      when 1 =>
+         return Left;
+
+      when 2 =>
+         return Num'Machine (Left * Left);
+
+      when 3 =>
+         return Num'Machine (Left * Left * Left);
+
+      when 4 =>
+         declare
+            Sqr : constant Num := Num'Machine (Left * Left);
+
+         begin
+            return Num'Machine (Sqr * Sqr);
+         end;
+
+      when Negative =>
+         return Num'Machine (1.0 / Exponr (Left, -Right));
+
+      when others =>
+         return Num'Machine (Expon (Left, Right));
+   end case;
+end System.Exponr;


diff --git /dev/null b/gcc/ada/libgnat/s-exponr.ads
new file mode 100644
--- /dev/null
+++ b/gcc/ada/libgnat/s-exponr.ads
@@ -0,0 +1,38 @@ 
+------------------------------------------------------------------------------
+--                                                                          --
+--                         GNAT RUN-TIME COMPONENTS                         --
+--                                                                          --
+--                        S Y S T E M . E X P O N R                         --
+--                                                                          --
+--                                 S p e c                                  --
+--                                                                          --
+--            Copyright (C) 2021, 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.      --
+--                                                                          --
+------------------------------------------------------------------------------
+
+--  Real exponentiation (checks off)
+
+generic
+
+   type Num is digits <>;
+
+function System.Exponr (Left : Num; Right : Integer) return Num;


diff --git a/gcc/ada/rtsfind.ads b/gcc/ada/rtsfind.ads
--- a/gcc/ada/rtsfind.ads
+++ b/gcc/ada/rtsfind.ads
@@ -228,6 +228,8 @@  package Rtsfind is
       System_Exception_Table,
       System_Exceptions_Debug,
       System_Exn_Int,
+      System_Exn_Flt,
+      System_Exn_LFlt,
       System_Exn_LLF,
       System_Exn_LLI,
       System_Exn_LLLI,
@@ -906,8 +908,10 @@  package Rtsfind is
 
      RE_Exn_Integer,                     -- System.Exn_Int
 
-     RE_Exn_Float,                       -- System.Exn_LLF
-     RE_Exn_Long_Float,                  -- System.Exn_LLF
+     RE_Exn_Float,                       -- System.Exn_Flt
+
+     RE_Exn_Long_Float,                  -- System.Exn_LFlt
+
      RE_Exn_Long_Long_Float,             -- System.Exn_LLF
 
      RE_Exn_Long_Long_Integer,           -- System.Exn_LLI
@@ -2592,8 +2596,10 @@  package Rtsfind is
 
      RE_Exn_Integer                      => System_Exn_Int,
 
-     RE_Exn_Float                        => System_Exn_LLF,
-     RE_Exn_Long_Float                   => System_Exn_LLF,
+     RE_Exn_Float                        => System_Exn_Flt,
+
+     RE_Exn_Long_Float                   => System_Exn_LFlt,
+
      RE_Exn_Long_Long_Float              => System_Exn_LLF,
 
      RE_Exn_Long_Long_Integer            => System_Exn_LLI,