fortran: Add IEEE_SIGNBIT and IEEE_FMA functions

The functions are added to the IEEE_ARITHMETIC module, but
are entirely expanded in the front-end, using GCC built-ins.

2022-08-31  Francois-Xavier Coudert <fxcoudert@gcc.gnu.org>

	PR fortran/95644

gcc/fortran/
	* f95-lang.cc (gfc_init_builtin_functions): Declare FMA
	built-ins.
	* mathbuiltins.def: Declare FMA built-ins.
	* trans-intrinsic.cc (conv_intrinsic_ieee_fma): New function.
	(conv_intrinsic_ieee_signbit): New function.
	(gfc_build_intrinsic_lib_fndecls): Add cases for FMA and
	SIGNBIT.

gcc/testsuite/
	* gfortran.dg/ieee/fma_1.f90: New test.
	* gfortran.dg/ieee/signbit_1.f90: New test.

libgfortran/
	* ieee/ieee_arithmetic.F90: Add IEEE_SIGNBIT and IEEE_FMA.

5 files changed, 332 insertions, 2 deletions

diff --git a/gcc/fortran/f95-lang.cc b/gcc/fortran/f95-lang.cc
index 10ac8a9..ff4bf80 100644
--- a/gcc/fortran/f95-lang.cc
+++ b/gcc/fortran/f95-lang.cc
@@ -1281,6 +1281,22 @@ gfc_init_builtin_functions (void)
 		      "__builtin_assume_aligned",
 		      ATTR_CONST_NOTHROW_LEAF_LIST);
 
+  ftype = build_function_type_list (long_double_type_node, long_double_type_node,
+				    long_double_type_node, long_double_type_node,
+				    NULL_TREE);
+  gfc_define_builtin ("__builtin_fmal", ftype, BUILT_IN_FMAL,
+		      "fmal", ATTR_CONST_NOTHROW_LEAF_LIST);
+  ftype = build_function_type_list (double_type_node, double_type_node,
+				    double_type_node, double_type_node,
+				    NULL_TREE);
+  gfc_define_builtin ("__builtin_fma", ftype, BUILT_IN_FMA,
+		      "fma", ATTR_CONST_NOTHROW_LEAF_LIST);
+  ftype = build_function_type_list (float_type_node, float_type_node,
+				    float_type_node, float_type_node,
+				    NULL_TREE);
+  gfc_define_builtin ("__builtin_fmaf", ftype, BUILT_IN_FMAF,
+		      "fmaf", ATTR_CONST_NOTHROW_LEAF_LIST);
+
   gfc_define_builtin ("__emutls_get_address",
 		      builtin_types[BT_FN_PTR_PTR],
 		      BUILT_IN_EMUTLS_GET_ADDRESS,
diff --git a/gcc/fortran/mathbuiltins.def b/gcc/fortran/mathbuiltins.def
index 615214e..9d55c34 100644
--- a/gcc/fortran/mathbuiltins.def
+++ b/gcc/fortran/mathbuiltins.def
@@ -60,6 +60,7 @@ OTHER_BUILTIN (CABS,      "cabs",      cabs,    true)
 OTHER_BUILTIN (COPYSIGN,  "copysign",  2,       true)
 OTHER_BUILTIN (CPOW,      "cpow",      cpow,    true)
 OTHER_BUILTIN (FABS,      "fabs",      1,       true)
+OTHER_BUILTIN (FMA,       "fma",       3,       true)
 OTHER_BUILTIN (FMOD,      "fmod",      2,       true)
 OTHER_BUILTIN (FREXP,     "frexp",     frexp,   false)
 OTHER_BUILTIN (LOGB,      "logb",      1,       true)
diff --git a/gcc/fortran/trans-intrinsic.cc b/gcc/fortran/trans-intrinsic.cc
index ec116ff..bb93802 100644
--- a/gcc/fortran/trans-intrinsic.cc
+++ b/gcc/fortran/trans-intrinsic.cc
@@ -695,7 +695,7 @@ gfc_build_intrinsic_lib_fndecls (void)
        C99-like library functions.  For now, we only handle _Float128
        q-suffixed or IEC 60559 f128-suffixed functions.  */
 
-    tree type, complex_type, func_1, func_2, func_cabs, func_frexp;
+    tree type, complex_type, func_1, func_2, func_3, func_cabs, func_frexp;
     tree func_iround, func_lround, func_llround, func_scalbn, func_cpow;
 
     memset (quad_decls, 0, sizeof(tree) * (END_BUILTINS + 1));
@@ -715,6 +715,8 @@ gfc_build_intrinsic_lib_fndecls (void)
 					     type, NULL_TREE);
     /* type (*) (type, type) */
     func_2 = build_function_type_list (type, type, type, NULL_TREE);
+    /* type (*) (type, type, type) */
+    func_3 = build_function_type_list (type, type, type, type, NULL_TREE);
     /* type (*) (type, &int) */
     func_frexp
       = build_function_type_list (type,
@@ -9781,7 +9783,7 @@ conv_ieee_function_args (gfc_se *se, gfc_expr *expr, tree *argarray,
 }
 
 
-/* Generate code for intrinsics IEEE_IS_NAN, IEEE_IS_FINITE,
+/* Generate code for intrinsics IEEE_IS_NAN, IEEE_IS_FINITE
    and IEEE_UNORDERED, which translate directly to GCC type-generic
    built-ins.  */
 
@@ -9801,6 +9803,23 @@ conv_intrinsic_ieee_builtin (gfc_se * se, gfc_expr * expr,
 }
 
 
+/* Generate code for intrinsics IEEE_SIGNBIT.  */
+
+static void
+conv_intrinsic_ieee_signbit (gfc_se * se, gfc_expr * expr)
+{
+  tree arg, signbit;
+
+  conv_ieee_function_args (se, expr, &arg, 1);
+  signbit = build_call_expr_loc (input_location,
+				 builtin_decl_explicit (BUILT_IN_SIGNBIT),
+				 1, arg);
+  signbit = fold_build2_loc (input_location, NE_EXPR, logical_type_node,
+			     signbit, integer_zero_node);
+  se->expr = fold_convert (gfc_typenode_for_spec (&expr->ts), signbit);
+}
+
+
 /* Generate code for IEEE_IS_NORMAL intrinsic:
      IEEE_IS_NORMAL(x) --> (__builtin_isnormal(x) || x == 0)  */
 
@@ -10207,6 +10226,30 @@ conv_intrinsic_ieee_value (gfc_se *se, gfc_expr *expr)
 }
 
 
+/* Generate code for IEEE_FMA.  */
+
+static void
+conv_intrinsic_ieee_fma (gfc_se * se, gfc_expr * expr)
+{
+  tree args[3], decl, call;
+  int argprec;
+
+  conv_ieee_function_args (se, expr, args, 3);
+
+  /* All three arguments should have the same type.  */
+  gcc_assert (TYPE_PRECISION (TREE_TYPE (args[0])) == TYPE_PRECISION (TREE_TYPE (args[1])));
+  gcc_assert (TYPE_PRECISION (TREE_TYPE (args[0])) == TYPE_PRECISION (TREE_TYPE (args[2])));
+
+  /* Call the type-generic FMA built-in.  */
+  argprec = TYPE_PRECISION (TREE_TYPE (args[0]));
+  decl = builtin_decl_for_precision (BUILT_IN_FMA, argprec);
+  call = build_call_expr_loc_array (input_location, decl, 3, args);
+
+  /* Convert to the final type.  */
+  se->expr = fold_convert (TREE_TYPE (args[0]), call);
+}
+
+
 /* Generate code for an intrinsic function from the IEEE_ARITHMETIC
    module.  */
 
@@ -10221,6 +10264,8 @@ gfc_conv_ieee_arithmetic_function (gfc_se * se, gfc_expr * expr)
     conv_intrinsic_ieee_builtin (se, expr, BUILT_IN_ISFINITE, 1);
   else if (startswith (name, "_gfortran_ieee_unordered"))
     conv_intrinsic_ieee_builtin (se, expr, BUILT_IN_ISUNORDERED, 2);
+  else if (startswith (name, "_gfortran_ieee_signbit"))
+    conv_intrinsic_ieee_signbit (se, expr);
   else if (startswith (name, "_gfortran_ieee_is_normal"))
     conv_intrinsic_ieee_is_normal (se, expr);
   else if (startswith (name, "_gfortran_ieee_is_negative"))
@@ -10241,6 +10286,8 @@ gfc_conv_ieee_arithmetic_function (gfc_se * se, gfc_expr * expr)
     conv_intrinsic_ieee_class (se, expr);
   else if (startswith (name, "ieee_value_") && ISDIGIT (name[11]))
     conv_intrinsic_ieee_value (se, expr);
+  else if (startswith (name, "_gfortran_ieee_fma"))
+    conv_intrinsic_ieee_fma (se, expr);
   else
     /* It is not among the functions we translate directly.  We return
        false, so a library function call is emitted.  */
diff --git a/gcc/testsuite/gfortran.dg/ieee/fma_1.f90 b/gcc/testsuite/gfortran.dg/ieee/fma_1.f90
new file mode 100644
index 0000000..3463642
--- /dev/null
+++ b/gcc/testsuite/gfortran.dg/ieee/fma_1.f90
@@ -0,0 +1,100 @@
+! Test IEEE_FMA
+! { dg-do run }
+
+  use, intrinsic :: ieee_features
+  use, intrinsic :: ieee_exceptions
+  use, intrinsic :: ieee_arithmetic
+  implicit none
+
+  integer :: ex
+
+  real :: sx1, sx2, sx3
+  double precision :: dx1, dx2, dx3
+
+  ! k1 and k2 will be large real kinds, if supported, and single/double
+  ! otherwise
+  integer, parameter :: k1 = &
+    max(ieee_selected_real_kind(precision(0.d0) + 1), kind(0.))
+  integer, parameter :: k2 = &
+    max(ieee_selected_real_kind(precision(0._k1) + 1), kind(0.d0))
+
+  real(kind=k1) :: lx1, lx2, lx3
+  real(kind=k2) :: wx1, wx2, wx3
+
+  ! Float
+
+  sx1 = 3 ; sx2 = 2 ; sx3 = 1
+  if (ieee_fma(sx1, sx2, sx3) /= 7) stop 1
+  sx1 = 0 ; sx2 = 2 ; sx3 = 1
+  if (ieee_fma(sx1, sx2, sx3) /= 1) stop 2
+  sx1 = 3 ; sx2 = 2 ; sx3 = 0
+  if (ieee_fma(sx1, sx2, sx3) /= 6) stop 3
+
+  ex = int(log(rrspacing(real(1, kind(sx1)))) / log(real(2, kind(sx1)))) - 1
+  sx1 = 1 + spacing(real(1, kind(sx1)))
+  sx2 = 2 ; sx2 = sx2 ** ex ; sx2 = sx2 * 3
+  sx3 = -sx2
+
+  print *, sx1 * sx2 + sx3
+  print *, ieee_fma(sx1, sx2, sx3)
+  if (ieee_fma(sx1, sx2, sx3) /= real(3, kind(sx1)) / 2) stop 4
+  !if (ieee_fma(sx1, sx2, sx3) == sx1 * sx2 + sx3) stop 5
+
+  ! Double
+
+  dx1 = 3 ; dx2 = 2 ; dx3 = 1
+  if (ieee_fma(dx1, dx2, dx3) /= 7) stop 1
+  dx1 = 0 ; dx2 = 2 ; dx3 = 1
+  if (ieee_fma(dx1, dx2, dx3) /= 1) stop 2
+  dx1 = 3 ; dx2 = 2 ; dx3 = 0
+  if (ieee_fma(dx1, dx2, dx3) /= 6) stop 3
+
+  ex = int(log(rrspacing(real(1, kind(dx1)))) / log(real(2, kind(dx1)))) - 1
+  dx1 = 1 + spacing(real(1, kind(dx1)))
+  dx2 = 2 ; dx2 = dx2 ** ex ; dx2 = dx2 * 3
+  dx3 = -dx2
+
+  print *, dx1 * dx2 + dx3
+  print *, ieee_fma(dx1, dx2, dx3)
+  if (ieee_fma(dx1, dx2, dx3) /= real(3, kind(dx1)) / 2) stop 4
+  !if (ieee_fma(dx1, dx2, dx3) == dx1 * dx2 + dx3) stop 5
+
+  ! Large kind 1
+
+  lx1 = 3 ; lx2 = 2 ; lx3 = 1
+  if (ieee_fma(lx1, lx2, lx3) /= 7) stop 1
+  lx1 = 0 ; lx2 = 2 ; lx3 = 1
+  if (ieee_fma(lx1, lx2, lx3) /= 1) stop 2
+  lx1 = 3 ; lx2 = 2 ; lx3 = 0
+  if (ieee_fma(lx1, lx2, lx3) /= 6) stop 3
+
+  ex = int(log(rrspacing(real(1, kind(lx1)))) / log(real(2, kind(lx1)))) - 1
+  lx1 = 1 + spacing(real(1, kind(lx1)))
+  lx2 = 2 ; lx2 = lx2 ** ex ; lx2 = lx2 * 3
+  lx3 = -lx2
+
+  print *, lx1 * lx2 + lx3
+  print *, ieee_fma(lx1, lx2, lx3)
+  if (ieee_fma(lx1, lx2, lx3) /= real(3, kind(lx1)) / 2) stop 4
+  if (ieee_fma(lx1, lx2, lx3) == lx1 * lx2 + lx3) stop 5
+
+  ! Large kind 2
+
+  wx1 = 3 ; wx2 = 2 ; wx3 = 1
+  if (ieee_fma(wx1, wx2, wx3) /= 7) stop 1
+  wx1 = 0 ; wx2 = 2 ; wx3 = 1
+  if (ieee_fma(wx1, wx2, wx3) /= 1) stop 2
+  wx1 = 3 ; wx2 = 2 ; wx3 = 0
+  if (ieee_fma(wx1, wx2, wx3) /= 6) stop 3
+
+  ex = int(log(rrspacing(real(1, kind(wx1)))) / log(real(2, kind(wx1)))) - 1
+  wx1 = 1 + spacing(real(1, kind(wx1)))
+  wx2 = 2 ; wx2 = wx2 ** ex ; wx2 = wx2 * 3
+  wx3 = -wx2
+
+  print *, wx1 * wx2 + wx3
+  print *, ieee_fma(wx1, wx2, wx3)
+  if (ieee_fma(wx1, wx2, wx3) /= real(3, kind(wx1)) / 2) stop 4
+  if (ieee_fma(wx1, wx2, wx3) == wx1 * wx2 + wx3) stop 5
+
+end
diff --git a/gcc/testsuite/gfortran.dg/ieee/signbit_1.f90 b/gcc/testsuite/gfortran.dg/ieee/signbit_1.f90
new file mode 100644
index 0000000..5d6e41d
--- /dev/null
+++ b/gcc/testsuite/gfortran.dg/ieee/signbit_1.f90
@@ -0,0 +1,166 @@
+! Test IEEE_SIGNBIT
+! { dg-do run }
+
+  use, intrinsic :: ieee_features
+  use, intrinsic :: ieee_exceptions
+  use, intrinsic :: ieee_arithmetic
+  implicit none
+
+  real :: sx1
+  double precision :: dx1
+
+  ! k1 and k2 will be large real kinds, if supported, and single/double
+  ! otherwise
+  integer, parameter :: k1 = &
+    max(ieee_selected_real_kind(precision(0.d0) + 1), kind(0.))
+  integer, parameter :: k2 = &
+    max(ieee_selected_real_kind(precision(0._k1) + 1), kind(0.d0))
+
+  real(kind=k1) :: xk1
+  real(kind=k2) :: xk2
+
+  ! Float
+
+  sx1 = 1.3
+  if (ieee_signbit(sx1)) stop 1
+  sx1 = huge(sx1)
+  if (ieee_signbit(sx1)) stop 2
+  sx1 = ieee_value(sx1, ieee_positive_inf)
+  if (ieee_signbit(sx1)) stop 3
+  sx1 = tiny(sx1)
+  if (ieee_signbit(sx1)) stop 4
+  sx1 = tiny(sx1)
+  sx1 = sx1 / 101
+  if (ieee_signbit(sx1)) stop 5
+  sx1 = 0
+  if (ieee_signbit(sx1)) stop 6
+  sx1 = ieee_value(sx1, ieee_quiet_nan)
+  if (ieee_signbit(sx1)) stop 7
+
+  sx1 = -1.3
+  if (.not. ieee_signbit(sx1)) stop 8
+  sx1 = -huge(sx1)
+  if (.not. ieee_signbit(sx1)) stop 9
+  sx1 = -ieee_value(sx1, ieee_positive_inf)
+  if (.not. ieee_signbit(sx1)) stop 10
+  sx1 = -tiny(sx1)
+  if (.not. ieee_signbit(sx1)) stop 11
+  sx1 = -tiny(sx1)
+  sx1 = sx1 / 101
+  if (.not. ieee_signbit(sx1)) stop 12
+  sx1 = 0
+  sx1 = -sx1
+  if (.not. ieee_signbit(sx1)) stop 13
+  sx1 = ieee_value(sx1, ieee_quiet_nan)
+  sx1 = -sx1
+  if (.not. ieee_signbit(sx1)) stop 14
+
+  ! Double
+
+  dx1 = 1.3
+  if (ieee_signbit(dx1)) stop 1
+  dx1 = huge(dx1)
+  if (ieee_signbit(dx1)) stop 2
+  dx1 = ieee_value(dx1, ieee_positive_inf)
+  if (ieee_signbit(dx1)) stop 3
+  dx1 = tiny(dx1)
+  if (ieee_signbit(dx1)) stop 4
+  dx1 = tiny(dx1)
+  dx1 = dx1 / 101
+  if (ieee_signbit(dx1)) stop 5
+  dx1 = 0
+  if (ieee_signbit(dx1)) stop 6
+  dx1 = ieee_value(dx1, ieee_quiet_nan)
+  if (ieee_signbit(dx1)) stop 7
+
+  dx1 = -1.3
+  if (.not. ieee_signbit(dx1)) stop 8
+  dx1 = -huge(dx1)
+  if (.not. ieee_signbit(dx1)) stop 9
+  dx1 = -ieee_value(dx1, ieee_positive_inf)
+  if (.not. ieee_signbit(dx1)) stop 10
+  dx1 = -tiny(dx1)
+  if (.not. ieee_signbit(dx1)) stop 11
+  dx1 = -tiny(dx1)
+  dx1 = dx1 / 101
+  if (.not. ieee_signbit(dx1)) stop 12
+  dx1 = 0
+  dx1 = -dx1
+  if (.not. ieee_signbit(dx1)) stop 13
+  dx1 = ieee_value(dx1, ieee_quiet_nan)
+  dx1 = -dx1
+  if (.not. ieee_signbit(dx1)) stop 14
+
+  ! Large kind 1
+
+  xk1 = 1.3
+  if (ieee_signbit(xk1)) stop 1
+  xk1 = huge(xk1)
+  if (ieee_signbit(xk1)) stop 2
+  xk1 = ieee_value(xk1, ieee_positive_inf)
+  if (ieee_signbit(xk1)) stop 3
+  xk1 = tiny(xk1)
+  if (ieee_signbit(xk1)) stop 4
+  xk1 = tiny(xk1)
+  xk1 = xk1 / 101
+  if (ieee_signbit(xk1)) stop 5
+  xk1 = 0
+  if (ieee_signbit(xk1)) stop 6
+  xk1 = ieee_value(xk1, ieee_quiet_nan)
+  if (ieee_signbit(xk1)) stop 7
+
+  xk1 = -1.3
+  if (.not. ieee_signbit(xk1)) stop 8
+  xk1 = -huge(xk1)
+  if (.not. ieee_signbit(xk1)) stop 9
+  xk1 = -ieee_value(xk1, ieee_positive_inf)
+  if (.not. ieee_signbit(xk1)) stop 10
+  xk1 = -tiny(xk1)
+  if (.not. ieee_signbit(xk1)) stop 11
+  xk1 = -tiny(xk1)
+  xk1 = xk1 / 101
+  if (.not. ieee_signbit(xk1)) stop 12
+  xk1 = 0
+  xk1 = -xk1
+  if (.not. ieee_signbit(xk1)) stop 13
+  xk1 = ieee_value(xk1, ieee_quiet_nan)
+  xk1 = -xk1
+  if (.not. ieee_signbit(xk1)) stop 14
+
+  ! Large kind 2
+
+  xk2 = 1.3
+  if (ieee_signbit(xk2)) stop 1
+  xk2 = huge(xk2)
+  if (ieee_signbit(xk2)) stop 2
+  xk2 = ieee_value(xk2, ieee_positive_inf)
+  if (ieee_signbit(xk2)) stop 3
+  xk2 = tiny(xk2)
+  if (ieee_signbit(xk2)) stop 4
+  xk2 = tiny(xk2)
+  xk2 = xk2 / 101
+  if (ieee_signbit(xk2)) stop 5
+  xk2 = 0
+  if (ieee_signbit(xk2)) stop 6
+  xk2 = ieee_value(xk2, ieee_quiet_nan)
+  if (ieee_signbit(xk2)) stop 7
+
+  xk2 = -1.3
+  if (.not. ieee_signbit(xk2)) stop 8
+  xk2 = -huge(xk2)
+  if (.not. ieee_signbit(xk2)) stop 9
+  xk2 = -ieee_value(xk2, ieee_positive_inf)
+  if (.not. ieee_signbit(xk2)) stop 10
+  xk2 = -tiny(xk2)
+  if (.not. ieee_signbit(xk2)) stop 11
+  xk2 = -tiny(xk2)
+  xk2 = xk2 / 101
+  if (.not. ieee_signbit(xk2)) stop 12
+  xk2 = 0
+  xk2 = -xk2
+  if (.not. ieee_signbit(xk2)) stop 13
+  xk2 = ieee_value(xk2, ieee_quiet_nan)
+  xk2 = -xk2
+  if (.not. ieee_signbit(xk2)) stop 14
+
+end