1 #include <math.h> 2 #include <errno.h> 3 #include <limits.h> 4 5 double pow(double x,double y)6pow(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 } 77