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