[PATCH 4/6] MIPS/math: Implement optimized fmaximum/fminmum(,_mag)(,f)

[YunQiang Su <syq@gcc.gnu.org>]

MIPSr6 defines max/maxa/min/mina instructions, which have slight
different with fmaximum/fminimum: return the number instead of NAN if
NUM vs NAN.  Detecting NAN is required before these instructions.

Another problem of the generic implemention for MIPS is that, it
use copysign for +0/-0 problem.  In fact we can use GET_HIGH_WORD or
GET_FLOAT_WORD, and then determine by the sign bit:

  int32_t xi;
  GET_HIGH_WORD (xi, x);
  return (xi < 0 ? y : x);

GET_HIGH_WORD/GET_FLOAT_WORD are much more friendly to MIPS FPU.
We can use `mfhc1` to GET_HIGH_WORD, and `mfc1` to GET_FLOAT_WORD.

Since the abs.fmt instructions will signal if one operand is qNAN or
sNaN, M_FABS, aka __builtin_fabs will issue at least 4 instructions
(mfc1/ext/ins/mtc1).  So detecting NaN is also required.

	* sysdeps/mips/ieee754/s_fmaximum.c
	* sysdeps/mips/ieee754/s_fmaximum_mag.c
	* sysdeps/mips/ieee754/s_fmaximum_magf.c
	* sysdeps/mips/ieee754/s_fmaximumf.c
	* sysdeps/mips/ieee754/s_fminimum.c
	* sysdeps/mips/ieee754/s_fminimum_mag.c
	* sysdeps/mips/ieee754/s_fminimum_magf.c
	* sysdeps/mips/ieee754/s_fminimumf.c

Signed-off-by: YunQiang Su <syq@gcc.gnu.org>
---
 sysdeps/mips/ieee754/s_fmaximum.c      | 48 ++++++++++++++++++++++
 sysdeps/mips/ieee754/s_fmaximum_mag.c  | 57 ++++++++++++++++++++++++++
 sysdeps/mips/ieee754/s_fmaximum_magf.c | 55 +++++++++++++++++++++++++
 sysdeps/mips/ieee754/s_fmaximumf.c     | 46 +++++++++++++++++++++
 sysdeps/mips/ieee754/s_fminimum.c      | 48 ++++++++++++++++++++++
 sysdeps/mips/ieee754/s_fminimum_mag.c  | 57 ++++++++++++++++++++++++++
 sysdeps/mips/ieee754/s_fminimum_magf.c | 55 +++++++++++++++++++++++++
 sysdeps/mips/ieee754/s_fminimumf.c     | 46 +++++++++++++++++++++
 8 files changed, 412 insertions(+)
 create mode 100644 sysdeps/mips/ieee754/s_fmaximum.c
 create mode 100644 sysdeps/mips/ieee754/s_fmaximum_mag.c
 create mode 100644 sysdeps/mips/ieee754/s_fmaximum_magf.c
 create mode 100644 sysdeps/mips/ieee754/s_fmaximumf.c
 create mode 100644 sysdeps/mips/ieee754/s_fminimum.c
 create mode 100644 sysdeps/mips/ieee754/s_fminimum_mag.c
 create mode 100644 sysdeps/mips/ieee754/s_fminimum_magf.c
 create mode 100644 sysdeps/mips/ieee754/s_fminimumf.c

diff --git a/sysdeps/mips/ieee754/s_fmaximum.c b/sysdeps/mips/ieee754/s_fmaximum.c
new file mode 100644
index 0000000000..5a1e6a0313
--- /dev/null
+++ b/sysdeps/mips/ieee754/s_fmaximum.c
@@ -0,0 +1,48 @@
+/* fmaximum().  MIPS version.
+   Copyright (C) 2024 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C 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
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, see
+   <https://www.gnu.org/licenses/>.  */
+
+#include <math.h>
+#include <libm-alias-double.h>
+#include <math_private.h>
+
+double
+__fmaximum (double x, double y)
+{
+#if __mips_isa_rev >= 6 && defined(__mips_hard_float)                         \
+    && !defined(__mips_single_float)
+  /* MAX.d returns NUM if NUM vs qNAN.  */
+  if (isunordered (x, y))
+    return x + y;
+  double ret;
+  asm volatile("max.d	%0, %1, %2" : "=f"(ret) : "f"(x), "f"(y));
+  return ret;
+#else
+  if (isgreater (x, y))
+    return x;
+  else if (isless (x, y))
+    return y;
+  if (isunordered (x, y))
+    return x + y;
+
+  int32_t xi;
+  GET_HIGH_WORD (xi, x);
+  return (xi < 0 ? y : x);
+#endif
+}
+
+libm_alias_double (__fmaximum, fmaximum)
diff --git a/sysdeps/mips/ieee754/s_fmaximum_mag.c b/sysdeps/mips/ieee754/s_fmaximum_mag.c
new file mode 100644
index 0000000000..0eac275167
--- /dev/null
+++ b/sysdeps/mips/ieee754/s_fmaximum_mag.c
@@ -0,0 +1,57 @@
+/* fmaximum_mag().  MIPS version.
+   Copyright (C) 2024 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C 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
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, see
+   <https://www.gnu.org/licenses/>.  */
+
+#include <math.h>
+#include <libm-alias-double.h>
+#include <math_private.h>
+#include <math-type-macros-double.h>
+
+double
+__fmaximum_mag (double x, double y)
+{
+  /* MAXA.d return NUM if NUM vs qNAN.  ABS.d signals both sNAN and qNAN on
+     pre-R5.  */
+  if (isunordered (x, y))
+    return x + y;
+#if __mips_isa_rev >= 6 && defined(__mips_hard_float)                         \
+    && !defined(__mips_single_float)
+  double ret;
+  asm volatile("maxa.d	%0, %1, %2" : "=f"(ret) : "f"(x), "f"(y));
+  return ret;
+#else
+  double ax;
+  double ay;
+#  if defined(__mips_hard_float) && !defined(__mips_single_float)
+  asm volatile("abs.d	%0, %1" : "=f"(ax) : "f"(x));
+  asm volatile("abs.d	%0, %1" : "=f"(ay) : "f"(y));
+#  else
+  ax = M_FABS (x);
+  ay = M_FABS (y);
+#  endif
+  if (isgreater (ax, ay))
+    return x;
+  else if (isless (ax, ay))
+    return y;
+
+  int32_t xi;
+  GET_HIGH_WORD (xi, x);
+  return (xi < 0 ? y : x);
+#endif
+}
+
+libm_alias_double (__fmaximum_mag, fmaximum_mag)
diff --git a/sysdeps/mips/ieee754/s_fmaximum_magf.c b/sysdeps/mips/ieee754/s_fmaximum_magf.c
new file mode 100644
index 0000000000..dd871bac07
--- /dev/null
+++ b/sysdeps/mips/ieee754/s_fmaximum_magf.c
@@ -0,0 +1,55 @@
+/* fmaximum_magf().  MIPS version.
+   Copyright (C) 2024 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C 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
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, see
+   <https://www.gnu.org/licenses/>.  */
+
+#include <math.h>
+#include <libm-alias-float.h>
+#include <math-type-macros-float.h>
+
+float
+__fmaximum_magf (float x, float y)
+{
+  /* MAXA.s return NUM if NUM vs qNAN.  ABS.s signals both sNAN and qNAN on
+     pre-R5.  */
+  if (isunordered (x, y))
+    return x + y;
+#if __mips_isa_rev >= 6 && defined(__mips_hard_float)
+  float ret;
+  asm volatile("maxa.s	%0, %1, %2" : "=f"(ret) : "f"(x), "f"(y));
+  return ret;
+#else
+  float ax;
+  float ay;
+#  if defined(__mips_hard_float)
+  asm volatile("abs.s	%0, %1" : "=f"(ax) : "f"(x));
+  asm volatile("abs.s	%0, %1" : "=f"(ay) : "f"(y));
+#  else
+  ax = M_FABS (x);
+  ay = M_FABS (y);
+#  endif
+  if (isgreater (ax, ay))
+    return x;
+  else if (isless (ax, ay))
+    return y;
+
+  int32_t xi;
+  GET_FLOAT_WORD (xi, x);
+  return (xi < 0 ? y : x);
+#endif
+}
+
+libm_alias_float (__fmaximum_mag, fmaximum_mag)
diff --git a/sysdeps/mips/ieee754/s_fmaximumf.c b/sysdeps/mips/ieee754/s_fmaximumf.c
new file mode 100644
index 0000000000..a266ee76b6
--- /dev/null
+++ b/sysdeps/mips/ieee754/s_fmaximumf.c
@@ -0,0 +1,46 @@
+/* fmaximumf().  MIPS version.
+   Copyright (C) 2024 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C 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
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, see
+   <https://www.gnu.org/licenses/>.  */
+
+#include <math.h>
+#include <libm-alias-float.h>
+
+float
+__fmaximumf (float x, float y)
+{
+#if __mips_isa_rev >= 6 && defined(__mips_hard_float)
+  /* MAX.s returns NUM if NUM vs qNAN.  */
+  if (isunordered (x, y))
+    return x + y;
+  float ret;
+  asm volatile("max.s	%0, %1, %2" : "=f"(ret) : "f"(x), "f"(y));
+  return ret;
+#else
+  if (isgreater (x, y))
+    return x;
+  else if (isless (x, y))
+    return y;
+  if (isunordered (x, y))
+    return x + y;
+
+  int32_t xi;
+  GET_FLOAT_WORD (xi, x);
+  return (xi < 0 ? y : x);
+#endif
+}
+
+libm_alias_float (__fmaximum, fmaximum)
diff --git a/sysdeps/mips/ieee754/s_fminimum.c b/sysdeps/mips/ieee754/s_fminimum.c
new file mode 100644
index 0000000000..083da390ae
--- /dev/null
+++ b/sysdeps/mips/ieee754/s_fminimum.c
@@ -0,0 +1,48 @@
+/* fminimum().  MIPS version.
+   Copyright (C) 2024 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C 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
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, see
+   <https://www.gnu.org/licenses/>.  */
+
+#include <math.h>
+#include <libm-alias-double.h>
+#include <math_private.h>
+
+double
+__fminimum (double x, double y)
+{
+#if __mips_isa_rev >= 6 && defined(__mips_hard_float)                         \
+    && !defined(__mips_single_float)
+  /* MIN.d returns NUM if NUM vs qNAN.  */
+  if (isunordered (x, y))
+    return x + y;
+  double ret;
+  asm volatile("min.d	%0, %1, %2" : "=f"(ret) : "f"(x), "f"(y));
+  return ret;
+#else
+  if (isgreater (x, y))
+    return y;
+  else if (isless (x, y))
+    return x;
+  if (isunordered (x, y))
+    return x + y;
+
+  int32_t xi;
+  GET_HIGH_WORD (xi, x);
+  return (xi < 0 ? x : y);
+#endif
+}
+
+libm_alias_double (__fminimum, fminimum)
diff --git a/sysdeps/mips/ieee754/s_fminimum_mag.c b/sysdeps/mips/ieee754/s_fminimum_mag.c
new file mode 100644
index 0000000000..7adaa1c279
--- /dev/null
+++ b/sysdeps/mips/ieee754/s_fminimum_mag.c
@@ -0,0 +1,57 @@
+/* fminimum_mag().  MIPS version.
+   Copyright (C) 2024 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C 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
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, see
+   <https://www.gnu.org/licenses/>.  */
+
+#include <math.h>
+#include <libm-alias-double.h>
+#include <math_private.h>
+#include <math-type-macros-double.h>
+
+double
+__fminimum_mag (double x, double y)
+{
+  /* MINA.d return NUM if NUM vs qNAN.  ABS.d signals both sNAN and qNAN on
+     pre-R5.  */
+  if (isunordered (x, y))
+    return x + y;
+#if __mips_isa_rev >= 6 && defined(__mips_hard_float)                         \
+    && !defined(__mips_single_float)
+  double ret;
+  asm volatile("mina.d	%0, %1, %2" : "=f"(ret) : "f"(x), "f"(y));
+  return ret;
+#else
+  double ax;
+  double ay;
+#  if defined(__mips_hard_float) && !defined(__mips_single_float)
+  asm volatile("abs.d	%0, %1" : "=f"(ax) : "f"(x));
+  asm volatile("abs.d	%0, %1" : "=f"(ay) : "f"(y));
+#  else
+  ax = M_FABS (x);
+  ay = M_FABS (y);
+#  endif
+  if (isgreater (ax, ay))
+    return y;
+  else if (isless (ax, ay))
+    return x;
+
+  int32_t xi;
+  GET_HIGH_WORD (xi, x);
+  return (xi < 0 ? x : y);
+#endif
+}
+
+libm_alias_double (__fminimum_mag, fminimum_mag)
diff --git a/sysdeps/mips/ieee754/s_fminimum_magf.c b/sysdeps/mips/ieee754/s_fminimum_magf.c
new file mode 100644
index 0000000000..6839e2914d
--- /dev/null
+++ b/sysdeps/mips/ieee754/s_fminimum_magf.c
@@ -0,0 +1,55 @@
+/* fminimum_magf().  MIPS version.
+   Copyright (C) 2024 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C 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
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, see
+   <https://www.gnu.org/licenses/>.  */
+
+#include <math.h>
+#include <libm-alias-float.h>
+#include <math-type-macros-float.h>
+
+float
+__fminimum_magf (float x, float y)
+{
+  /* MAXA.s return NUM if NUM vs qNAN.  ABS.s signals both sNAN and qNAN on
+     pre-R5.  */
+  if (isunordered (x, y))
+    return x + y;
+#if __mips_isa_rev >= 6 && defined(__mips_hard_float)
+  float ret;
+  asm volatile("mina.s	%0, %1, %2" : "=f"(ret) : "f"(x), "f"(y));
+  return ret;
+#else
+  float ax;
+  float ay;
+#  if defined(__mips_hard_float)
+  asm volatile("abs.s	%0, %1" : "=f"(ax) : "f"(x));
+  asm volatile("abs.s	%0, %1" : "=f"(ay) : "f"(y));
+#  else
+  ax = M_FABS (x);
+  ay = M_FABS (y);
+#  endif
+  if (isgreater (ax, ay))
+    return y;
+  else if (isless (ax, ay))
+    return x;
+
+  int32_t xi;
+  GET_FLOAT_WORD (xi, x);
+  return (xi < 0 ? x : y);
+#endif
+}
+
+libm_alias_float (__fminimum_mag, fminimum_mag)
diff --git a/sysdeps/mips/ieee754/s_fminimumf.c b/sysdeps/mips/ieee754/s_fminimumf.c
new file mode 100644
index 0000000000..f37ca1c23b
--- /dev/null
+++ b/sysdeps/mips/ieee754/s_fminimumf.c
@@ -0,0 +1,46 @@
+/* fminimumf().  MIPS version.
+   Copyright (C) 2024 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C 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
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, see
+   <https://www.gnu.org/licenses/>.  */
+
+#include <math.h>
+#include <libm-alias-float.h>
+
+float
+__fminimumf (float x, float y)
+{
+#if __mips_isa_rev >= 6 && defined(__mips_hard_float)
+  /* MIN.s returns NUM if NUM vs qNAN.  */
+  if (isunordered (x, y))
+    return x + y;
+  float ret;
+  asm volatile("min.s	%0, %1, %2" : "=f"(ret) : "f"(x), "f"(y));
+  return ret;
+#else
+  if (isgreater (x, y))
+    return y;
+  else if (isless (x, y))
+    return x;
+  if (isunordered (x, y))
+    return x + y;
+
+  int xi;
+  GET_FLOAT_WORD (xi, x);
+  return (xi < 0 ? x : y);
+#endif
+}
+
+libm_alias_float (__fminimum, fminimum)
-- 
2.39.2