toolchain/gcc/newlib/libm/mathfp/sf_pow.c

   1
   2 /* @(#)z_powf.c 1.0 98/08/13 */
   3 #include <float.h>
   4 #include "fdlibm.h"
   5 #include "zmath.h"
   6
   7 float powf (float x, float y)
   8 {
   9   float d, k, t, r = 1.0;
  10   int n, sign, exponent_is_even_int = 0;
  11   __int32_t px;
  12
  13   GET_FLOAT_WORD (px, x);
  14
  15   k = modff (y, &d);
  16
  17   if (k == 0.0)
  18     {
  19       /* Exponent y is an integer. */
  20       if (modff (ldexpf (y, -1), &t))
  21         {
  22           /* y is odd. */
  23           exponent_is_even_int = 0;
  24         }
  25       else
  26         {
  27           /* y is even. */
  28           exponent_is_even_int = 1;
  29         }
  30     }
  31
  32   if (x == 0.0 && y <= 0.0)
  33     {
  34       errno = EDOM;
  35     }
  36   else if ((t = y * log (fabsf (x))) >= BIGX)
  37     {
  38       errno = ERANGE;
  39       if (px & 0x80000000)
  40         {
  41           /* x is negative. */
  42           if (k)
  43             {
  44               /* y is not an integer. */
  45               errno = EDOM;
  46               x = 0.0;
  47             }
  48           else if (exponent_is_even_int)
  49             x = z_infinity_f.f;
  50           else
  51             x = -z_infinity_f.f;
  52         }
  53     else
  54       {
  55         x = z_infinity_f.f;
  56       }
  57     }
  58   else if (t < SMALLX)
  59     {
  60       errno = ERANGE;
  61       x = 0.0;
  62     }
  63   else
  64     {
  65       if ( !k && fabsf (d) <= 32767 )
  66         {
  67           n = (int) d;
  68
  69           if ((sign = (n < 0)))
  70             n = -n;
  71
  72           while ( n > 0 )
  73             {
  74               if ((unsigned int) n % 2)
  75                 r *= x;
  76               x *= x;
  77               n = (unsigned int) n / 2;
  78             }
  79
  80           if (sign)
  81             r = 1.0 / r;
  82
  83           return r;
  84         }
  85       else
  86         {
  87           if ( px & 0x80000000 )
  88             {
  89               /* x is negative. */
  90               if (k)
  91                 {
  92                   /* y is not an integer. */
  93                   errno = EDOM;
  94                   return 0.0;
  95                 }
  96             }
  97
  98           x = exp (t);
  99
 100           if (!exponent_is_even_int)
 101             {
 102               if (px & 0x80000000)
 103                 {
 104                   /* y is an odd integer, and x is negative,
 105                      so the result is negative. */
 106                   GET_FLOAT_WORD (px, x);
 107                   px |= 0x80000000;
 108                   SET_FLOAT_WORD (x, px);
 109                 }
 110             }
 111         }
 112     }
 113
 114   return x;
 115 }
 116
 117 #ifdef _DOUBLE_IS_32BITS
 118
 119 double pow (double x, double y)
 120 {
 121   return (double) powf ((float) x, (float) y);
 122 }
 123
 124 #endif /* defined(_DOUBLE_IS_32BITS) */