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