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