Implement C23 pown

C23 adds various <math.h> function families originally defined in TS
18661-4.  Add the pown functions, which are like pow but with an
integer exponent.  That exponent has type long long int in C23; it was
intmax_t in TS 18661-4, and as with other interfaces changed after
their initial appearance in the TS, I don't think we need to support
the original version of the interface.  The test inputs are based on
the subset of test inputs for pow that use integer exponents that fit
in long long.

As the first such template implementation that saves and restores the
rounding mode internally (to avoid possible issues with directed
rounding and intermediate overflows or underflows in the wrong
rounding mode), support also needed to be added for using
SET_RESTORE_ROUND* in such template function implementations.  This
required math-type-macros-float128.h to include <fenv_private.h>, so
it can tell whether SET_RESTORE_ROUNDF128 is defined.  In turn, the
include order with <fenv_private.h> included before <math_private.h>
broke loongarch builds, showing up that
sysdeps/loongarch/math_private.h is really a fenv_private.h file
(maybe implemented internally before the consistent split of those
headers in 2018?) and needed to be renamed to fenv_private.h to avoid
errors with duplicate macro definitions if <math_private.h> is
included after <fenv_private.h>.

The underlying implementation uses __ieee754_pow functions (called
more than once in some cases, where the exponent does not fit in the
floating type).  I expect a custom implementation for a given format,
that only handles integer exponents but handles larger exponents
directly, could be faster and more accurate in some cases.

I encourage searching for worst cases for ulps error for these
implementations (necessarily non-exhaustively, given the size of the
input space).

Tested for x86_64 and x86, and with build-many-glibcs.py.

55 files changed, 14133 insertions, 7 deletions

diff --git a/NEWS b/NEWS
index c2a978b460..52cbf1a989 100644
--- a/NEWS
+++ b/NEWS
@@ -14,7 +14,7 @@ Major new features:
   functions for float, double, long double, _FloatN and _FloatNx, and a
   type-generic macro in <tgmath.h>.
 
-  - Power and absolute-value functions: powr, rsqrt.
+  - Power and absolute-value functions: pown, powr, rsqrt.
 
 * On Linux, the pthread_gettid_np function has been added.
 
diff --git a/manual/math.texi b/manual/math.texi
index 7490693c0d..adbe271f36 100644
--- a/manual/math.texi
+++ b/manual/math.texi
@@ -789,6 +789,20 @@ but this is valid for @code{pow} if @var{power} is an integer.
 The @code{powr} functions are from TS 18661-4:2015.
 @end deftypefun
 
+@deftypefun double pown (double @var{base}, long long int @var{power})
+@deftypefunx float pownf (float @var{base}, long long int @var{power})
+@deftypefunx {long double} pownl (long double @var{base}, long long int @var{power})
+@deftypefunx _FloatN pownfN (_Float@var{N} @var{base}, long long int @var{power})
+@deftypefunx _FloatNx pownfNx (_Float@var{N}x @var{base}, long long int @var{power})
+@standards{TS 18661-4:2015, math.h}
+@safety{@prelim{}@mtsafe{}@assafe{}@acsafe{}}
+These return @var{base} raised to @var{power} (an integer).
+
+The @code{pown} functions are from TS 18661-4:2015 (which used
+@code{intmax_t} as the type of @var{power}; the type changed to
+@code{long long int} in C23).
+@end deftypefun
+
 @cindex square root function
 @deftypefun double sqrt (double @var{x})
 @deftypefunx float sqrtf (float @var{x})
diff --git a/math/Makefile b/math/Makefile
index 64d9d4d1e5..6b2331d146 100644
--- a/math/Makefile
+++ b/math/Makefile
@@ -139,6 +139,7 @@ gen-libm-calls = \
   s_log2p1F \
   s_nanF \
   s_nextdownF \
+  s_pownF \
   s_powrF \
   s_rsqrtF \
   s_significandF \
@@ -686,6 +687,7 @@ libm-test-funcs-auto = \
   log1p \
   log2p1 \
   pow \
+  pown \
   powr \
   rsqrt \
   sin \
@@ -1018,6 +1020,7 @@ tgmath3-macros = \
   nexttoward \
   nextup \
   pow \
+  pown \
   powr \
   remainder \
   remquo \
@@ -1442,6 +1445,7 @@ CFLAGS-s_nexttoward.c += -fno-builtin-nexttoward -fno-builtin-nexttowardl
 CFLAGS-s_nexttowardf.c += -fno-builtin-nexttowardf
 CFLAGS-s_nextup.c += -fno-builtin-nextupl
 CFLAGS-e_pow.c += -fno-builtin-powl
+CFLAGS-s_pown.c += -fno-builtin-pownl
 CFLAGS-s_powr.c += -fno-builtin-powrl
 CFLAGS-w_remainder.c += -fno-builtin-remainderl -fno-builtin-dreml
 CFLAGS-s_remquo.c += -fno-builtin-remquol
@@ -1578,6 +1582,7 @@ CFLAGS-s_nextafter.c += -fno-builtin-nextafterf32x -fno-builtin-nextafterf64
 CFLAGS-s_nextdown.c += -fno-builtin-nextdownf32x -fno-builtin-nextdownf64
 CFLAGS-s_nextup.c += -fno-builtin-nextupf32x -fno-builtin-nextupf64
 CFLAGS-e_pow.c += -fno-builtin-powf32x -fno-builtin-powf64
+CFLAGS-s_pown.c += -fno-builtin-pownf32x -fno-builtin-pownf64
 CFLAGS-s_powr.c += -fno-builtin-powrf32x -fno-builtin-powrf64
 CFLAGS-w_remainder.c += -fno-builtin-remainderf32x -fno-builtin-remainderf64
 CFLAGS-s_remquo.c += -fno-builtin-remquof32x -fno-builtin-remquof64
@@ -1705,6 +1710,7 @@ CFLAGS-s_nextafterf.c += -fno-builtin-nextafterf32
 CFLAGS-s_nextdownf.c += -fno-builtin-nextdownf32
 CFLAGS-s_nextupf.c += -fno-builtin-nextupf32
 CFLAGS-e_powf.c += -fno-builtin-powf32
+CFLAGS-s_pownf.c += -fno-builtin-pownf32
 CFLAGS-s_powrf.c += -fno-builtin-powrf32
 CFLAGS-w_remainderf.c += -fno-builtin-remainderf32
 CFLAGS-s_remquof.c += -fno-builtin-remquof32
diff --git a/math/Versions b/math/Versions
index f92707545a..36c0f23da6 100644
--- a/math/Versions
+++ b/math/Versions
@@ -673,9 +673,11 @@ libm {
   }
   GLIBC_2.42 {
     # Functions not involving _Float64x or _Float128, for all configurations.
+    pown; pownf; pownl; pownf32; pownf64; pownf32x;
     powr; powrf; powrl; powrf32; powrf64; powrf32x;
     rsqrt; rsqrtf; rsqrtl; rsqrtf32; rsqrtf64; rsqrtf32x;
     # Functions involving _Float64x or _Float128, for some configurations.
+    pownf64x; pownf128;
     powrf64x; powrf128;
     rsqrtf64x; rsqrtf128;
   }
diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in
index cb5d6443b1..7cb8a896ab 100644
--- a/math/auto-libm-test-in
+++ b/math/auto-libm-test-in
@@ -8151,6 +8151,383 @@ pow 0x1.059c76p+0 0x1.ff80bep+11
 pow 0x1.7ac7cp+5 23
 pow -0x1.7ac7cp+5 23
 
+pown 0 0
+pown 0 -0
+pown -0 0
+pown -0 -0
+
+pown 10 0
+pown 10 -0
+pown -10 0
+pown -10 -0
+
+pown 1 1
+pown 1 -1
+pown 1 0x4fffffffffffffff
+
+# pown (x, +-0) == 1.
+pown 32.75 0
+pown 32.75 -0
+pown -32.75 0
+pown -32.75 -0
+pown 0x1p72 0
+pown 0x1p72 -0
+pown 0x1p-72 0
+pown 0x1p-72 -0
+
+pown 0 1
+pown 0 11
+
+pown -0 1
+pown -0 11
+
+pown 0 2
+
+pown -0 2
+
+# pown (+0, y) == +0 for y an odd integer > 0.
+pown 0.0 27
+pown 0.0 0xffffff
+pown 0.0 0x1fffffffffffff
+
+# pown (-0, y) == -0 for y an odd integer > 0.
+pown -0 27
+pown -0 0xffffff
+pown -0 0x1fffffe
+pown -0 0x1fffffffffffff
+pown -0 0x3ffffffffffffe
+pown -0 0x7fffffffffffffff
+
+# pown (+0, y) == +0 for y > 0 and not an odd integer.
+pown 0.0 4
+pown 0.0 0x1000000
+
+# pown (-0, y) == +0 for y > 0 and not an odd integer.
+pown -0 4
+pown -0 0x1000000
+
+pown 2 4
+pown 256 8
+
+pown -1.0 -0xffffff
+pown -1.0 -0x1fffffe
+pown -1.0 -0x1fffffffffffff
+pown -1.0 -0x3ffffffffffffe
+pown -1.0 -0x7fffffffffffffff
+
+pown -1.0 0xffffff
+pown -1.0 0x1fffffe
+pown -1.0 0x1fffffffffffff
+pown -1.0 0x3ffffffffffffe
+pown -1.0 0x7fffffffffffffff
+
+pown -2.0 126
+pown -2.0 127
+pown -2.0 -126
+pown -2.0 -127
+
+pown -2.0 -0xffffff
+pown -2.0 -0x1fffffe
+pown -2.0 -0x1fffffffffffff
+pown -2.0 -0x3ffffffffffffe
+pown -2.0 -0x7fffffffffffffff
+
+pown -2.0 0xffffff
+pown -2.0 0x1fffffe
+pown -2.0 0x1fffffffffffff
+pown -2.0 0x3ffffffffffffe
+pown -2.0 0x7fffffffffffffff
+
+pown -max -2
+pown -max -3
+pown -max 2
+pown -max 3
+
+pown -max -0xffffff
+pown -max