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