sys/src/ape/lib/ap/math/pow.c

   1 #include <math.h>
   2 #include <errno.h>
   3 #include <limits.h>
   4
   5 double
   6 pow(double x, double y) /* return x ^ y (exponentiation) */
   7 {
   8         double xy, y1, ye;
   9         long i;
  10         int ex, ey, flip;
  11
  12         if(y == 0.0)
  13                 return 1.0;
  14
  15         flip = 0;
  16         if(y < 0.){
  17                 y = -y;
  18                 flip = 1;
  19         }
  20         y1 = modf(y, &ye);
  21         if(y1 != 0.0){
  22                 if(x <= 0.)
  23                         goto zreturn;
  24                 if(y1 > 0.5) {
  25                         y1 -= 1.;
  26                         ye += 1.;
  27                 }
  28                 xy = exp(y1 * log(x));
  29         }else
  30                 xy = 1.0;
  31         if(ye > LONG_MAX){
  32                 if(x <= 0.){
  33  zreturn:
  34                         if(x || flip)
  35                                 errno = EDOM;
  36                         return 0.;
  37                 }
  38                 if(flip){
  39                         if(y == .5)
  40                                 return 1./sqrt(x);
  41                         y = -y;
  42                 }else if(y == .5)
  43                         return sqrt(x);
  44                 return exp(y * log(x));
  45         }
  46         x = frexp(x, &ex);
  47         ey = 0;
  48         i = ye;
  49         if(i)
  50                 for(;;){
  51                         if(i & 1){
  52                                 xy *= x;
  53                                 ey += ex;
  54                         }
  55                         i >>= 1;
  56                         if(i == 0)
  57                                 break;
  58                         x *= x;
  59                         ex <<= 1;
  60                         if(x < .5){
  61                                 x += x;
  62                                 ex -= 1;
  63                         }
  64                 }
  65         if(flip){
  66                 xy = 1. / xy;
  67                 ey = -ey;
  68         }
  69         errno = 0;
  70         x = ldexp(xy, ey);
  71         if(errno && ey < 0) {
  72                 errno = 0;
  73                 x = 0.;
  74         }
  75         return x;
  76 }