1 files changed, 91 insertions, 0 deletions
diff --git a/math/s_pown_template.c b/math/s_pown_template.c
new file mode 100644
index 0000000..e46bf96
--- /dev/null
+++ b/math/s_pown_template.c
@@ -0,0 +1,91 @@
+/* Return X^Y for integer Y.
+   Copyright (C) 2025 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 <errno.h>
+#include <fenv_private.h>
+#include <limits.h>
+#include <math.h>
+#include <math-barriers.h>
+#include <math-narrow-eval.h>
+#include <math_private.h>
+#include <stdlib.h>
+
+FLOAT
+M_DECL_FUNC (__pown) (FLOAT x, long long int y)
+{
+  FLOAT ret;
+#if M_MANT_DIG >= LLONG_WIDTH - 1
+  ret = M_SUF (__ieee754_pow) (x, y);
+#else
+  if ((y <= 1LL << M_MANT_DIG
+       && y >= -(1LL << M_MANT_DIG))
+      || y == LLONG_MIN)
+    ret = M_SUF (__ieee754_pow) (x, y);
+  else
+    {
+      {
+	/* To avoid problems with overflow and underflow from
+	   intermediate computations occurring in the wrong directed
+	   rounding mode (when X is negative and Y is odd), do these
+	   computations in round-to-nearest mode and correct
+	   overflowing and underflowing results afterwards.  */
+	M_SET_RESTORE_ROUND (FE_TONEAREST);
+	ret = M_LIT (1.0);
+	while (y != 0)
+	  {
+	    long long int absy = llabs (y);
+	    int bits = LLONG_WIDTH - __builtin_clzll (absy);
+	    if (bits <= M_MANT_DIG)
+	      {
+		ret *= M_SUF (__ieee754_pow) (x, y);
+		y = 0;
+	      }
+	    else
+	      {
+		long long int absy_high
+		  = absy & (((1ULL << M_MANT_DIG) - 1) << (bits - M_MANT_DIG));
+		long long int absy_low = absy - absy_high;
+		long long int y_high = y < 0 ? -absy_high : absy_high;
+		ret *= M_SUF (__ieee754_pow) (x, y_high);
+		y = y < 0 ? -absy_low : absy_low;
+	      }
+	  }
+	ret = math_narrow_eval (ret);
+	math_force_eval (ret);
+      }
+      if (isfinite (x) && x != M_LIT (0.0))
+	{
+	  if (isinf (ret))
+	    ret = math_narrow_eval (M_COPYSIGN (M_MAX, ret) * M_MAX);
+	  else if (ret == M_LIT (0.0))
+	    ret = math_narrow_eval (M_COPYSIGN (M_MIN, ret) * M_MIN);
+	}
+    }
+#endif
+  if (!isfinite (ret))
+    {
+      if (isfinite (x))
+	__set_errno (ERANGE);
+    }
+  else if (ret == 0
+	   && isfinite (x)
+	   && x != 0)
+    __set_errno (ERANGE);
+  return ret;
+}
+declare_mgen_alias (__pown, pown);