newlib/libm/mathfp/s_pow.c

   1
   2 /* @(#)z_pow.c 1.0 98/08/13 */
   3
   4 /*
   5 FUNCTION
   6         <<pow>>, <<powf>>---x to the power y
   7 INDEX
   8         pow
   9 INDEX
  10         powf
  11
  12
  13 SYNOPSIS
  14         #include <math.h>
  15         double pow(double <[x]>, double <[y]>);
  16         float pow(float <[x]>, float <[y]>);
  17
  18 DESCRIPTION
  19         <<pow>> and <<powf>> calculate <[x]> raised to the exponent <[y]>.
  20         @tex
  21         (That is, $x^y$.)
  22         @end tex
  23
  24 RETURNS
  25         On success, <<pow>> and <<powf>> return the value calculated.
  26
  27         When the argument values would produce overflow, <<pow>>
  28         returns <<HUGE_VAL>> and set <<errno>> to <<ERANGE>>.  If the
  29         argument <[x]> passed to <<pow>> or <<powf>> is a negative
  30         noninteger, and <[y]> is also not an integer, then <<errno>>
  31         is set to <<EDOM>>.  If <[x]> and <[y]> are both 0, then
  32         <<pow>> and <<powf>> return <<1>>.
  33
  34 PORTABILITY
  35         <<pow>> is ANSI C. <<powf>> is an extension.  */
  36
  37 #include <float.h>
  38 #include "fdlibm.h"
  39 #include "zmath.h"
  40
  41 #ifndef _DOUBLE_IS_32BITS
  42
  43 double pow (double x, double y)
  44 {
  45   double d, k, t, r = 1.0;
  46   int n, sign, exponent_is_even_int = 0;
  47   __uint32_t px;
  48
  49   GET_HIGH_WORD (px, x);
  50
  51   k = modf (y, &d);
  52
  53   if (k == 0.0)
  54     {
  55       /* Exponent y is an integer. */
  56       if (modf (ldexp (y, -1), &t))
  57         {
  58           /* y is odd. */
  59           exponent_is_even_int = 0;
  60         }
  61       else
  62         {
  63           /* y is even. */
  64           exponent_is_even_int = 1;
  65         }
  66     }
  67
  68   if (x == 0.0)
  69     {
  70       if (y <= 0.0)
  71         errno = EDOM;
  72     }
  73   else if ((t = y * log (fabs (x))) >= BIGX)
  74     {
  75       errno = ERANGE;
  76       if (px & 0x80000000)
  77         {
  78           /* x is negative. */
  79           if (k)
  80             {
  81               /* y is not an integer. */
  82               errno = EDOM;
  83               x = 0.0;
  84             }
  85           else if (exponent_is_even_int)
  86             x = z_infinity.d;
  87           else
  88             x = -z_infinity.d;
  89         }
  90       else
  91         {
  92           x = z_infinity.d;
  93         }
  94     }
  95   else if (t < SMALLX)
  96     {
  97       errno = ERANGE;
  98       x = 0.0;
  99     }
 100   else
 101     {
 102       if ( !k && fabs(d) <= 32767 )
 103         {
 104           n = (int) d;
 105
 106           if ((sign = (n < 0)))
 107             n = -n;
 108
 109           while ( n > 0 )
 110             {
 111               if ((unsigned int) n % 2)
 112                 r *= x;
 113               x *= x;
 114               n = (unsigned int) n / 2;
 115             }
 116
 117           if (sign)
 118             r = 1.0 / r;
 119
 120           return r;
 121         }
 122       else
 123         {
 124           if ( px & 0x80000000 )
 125             {
 126               /* x is negative. */
 127               if ( k )
 128                 {
 129                   /* y is not an integer. */
 130                   errno = EDOM;
 131                   return 0.0;
 132                 }
 133             }
 134
 135           x = exp (t);
 136
 137           if (!exponent_is_even_int)
 138             {
 139               if (px & 0x80000000)
 140                 {
 141                   /* y is an odd integer, and x is negative,
 142                      so the result is negative. */
 143                   GET_HIGH_WORD (px, x);
 144                   px |= 0x80000000;
 145                   SET_HIGH_WORD (x, px);
 146                 }
 147             }
 148         }
 149     }
 150
 151   return x;
 152 }
 153
 154 #endif _DOUBLE_IS_32BITS