1  /* Copyright (C) 1997-2023 Free Software Foundation, Inc.
       2     This file is part of the GNU C Library.
       3  
       4     The GNU C Library is free software; you can redistribute it and/or
       5     modify it under the terms of the GNU Lesser General Public
       6     License as published by the Free Software Foundation; either
       7     version 2.1 of the License, or (at your option) any later version.
       8  
       9     The GNU C Library is distributed in the hope that it will be useful,
      10     but WITHOUT ANY WARRANTY; without even the implied warranty of
      11     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
      12     Lesser General Public License for more details.
      13  
      14     You should have received a copy of the GNU Lesser General Public
      15     License along with the GNU C Library.  If not, see
      16     <https://www.gnu.org/licenses/>.  */
      17  
      18  #include <math.h>
      19  #include <math_private.h>
      20  #include "mathimpl.h"
      21  #include <libm-alias-finite.h>
      22  
      23  #ifndef SUFF
      24  #define SUFF
      25  #endif
      26  #ifndef float_type
      27  #define float_type double
      28  #endif
      29  
      30  #define CONCATX(a,b) __CONCAT(a,b)
      31  #define s(name) CONCATX(name,SUFF)
      32  #define m81(func) __m81_u(s(func))
      33  
      34  float_type
      35  s(__ieee754_pow) (float_type x, float_type y)
      36  {
      37    float_type z;
      38    float_type ax;
      39    unsigned long x_cond, y_cond;
      40  
      41    y_cond = __m81_test (y);
      42    if (y_cond & __M81_COND_ZERO)
      43      return 1.0;
      44    if (y_cond & __M81_COND_NAN)
      45      return x == 1.0 ? x : x + y;
      46  
      47    x_cond = __m81_test (x);
      48    if (x_cond & __M81_COND_NAN)
      49      return x + y;
      50  
      51    if (y_cond & __M81_COND_INF)
      52      {
      53        ax = s(fabs) (x);
      54        if (ax == 1.0)
      55  	return ax;
      56        if (ax > 1.0)
      57  	return y_cond & __M81_COND_NEG ? 0 : y;
      58        else
      59  	return y_cond & __M81_COND_NEG ? -y : 0;
      60      }
      61  
      62    if (s(fabs) (y) == 1.0)
      63      return y_cond & __M81_COND_NEG ? 1 / x : x;
      64  
      65    if (y == 2)
      66      return x * x;
      67    if (y == 0.5 && !(x_cond & __M81_COND_NEG))
      68      return m81(sqrt) (x);
      69  
      70    if (x == 10.0)
      71      {
      72        __asm ("ftentox%.x %1, %0" : "=f" (z) : "f" (y));
      73        return z;
      74      }
      75    if (x == 2.0)
      76      {
      77        __asm ("ftwotox%.x %1, %0" : "=f" (z) : "f" (y));
      78        return z;
      79      }
      80  
      81    ax = s(fabs) (x);
      82    if (x_cond & (__M81_COND_INF | __M81_COND_ZERO) || ax == 1.0)
      83      {
      84        z = ax;
      85        if (y_cond & __M81_COND_NEG)
      86  	z = 1 / z;
      87        if (x_cond & __M81_COND_NEG)
      88  	{
      89  	  if (y != m81(__rint) (y))
      90  	    {
      91  	      if (x == -1)
      92  		z = (z - z) / (z - z);
      93  	    }
      94  	  else
      95  	    goto maybe_negate;
      96  	}
      97        return z;
      98      }
      99  
     100    if (x_cond & __M81_COND_NEG)
     101      {
     102        if (y == m81(__rint) (y))
     103  	{
     104  	  z = m81(__ieee754_exp) (y * m81(__ieee754_log) (-x));
     105  	maybe_negate:
     106  	  /* We always use the long double format, since y is already in
     107  	     this format and rounding won't change the result.  */
     108  	  {
     109  	    int32_t exponent;
     110  	    uint32_t i0, i1;
     111  	    GET_LDOUBLE_WORDS (exponent, i0, i1, y);
     112  	    exponent = (exponent & 0x7fff) - 0x3fff;
     113  	    if (exponent <= 31
     114  		? i0 & (1 << (31 - exponent))
     115  		: (exponent <= 63
     116  		   && i1 & (1 << (63 - exponent))))
     117  	      z = -z;
     118  	  }
     119  	}
     120        else
     121  	z = (y - y) / (y - y);
     122      }
     123    else
     124      z = m81(__ieee754_exp) (y * m81(__ieee754_log) (x));
     125    return z;
     126  }
     127  libm_alias_finite (s(__ieee754_pow), s (__pow))