ref: 1018d7a7091fd5e0433f501a3474cdc3fec322cc
dir: /libkern/pow.c/
#include <lib9.h> double pow(double x, double y) /* return x ^ y (exponentiation) */ { double xy, y1, ye; long i; int ex, ey, flip; if(y == 0.0) return 1.0; flip = 0; if(y < 0.){ y = -y; flip = 1; } y1 = modf(y, &ye); if(y1 != 0.0){ if(x <= 0.) goto zreturn; if(y1 > 0.5) { y1 -= 1.; ye += 1.; } xy = exp(y1 * log(x)); }else xy = 1.0; if(ye > 0x7FFFFFFF){ /* should be ~0UL but compiler can't convert double to ulong */ if(x <= 0.){ zreturn: if(x==0. && !flip) return 0.; return NaN(); } if(flip){ if(y == .5) return 1./sqrt(x); y = -y; }else if(y == .5) return sqrt(x); return exp(y * log(x)); } x = frexp(x, &ex); ey = 0; i = ye; if(i) for(;;){ if(i & 1){ xy *= x; ey += ex; } i >>= 1; if(i == 0) break; x *= x; ex <<= 1; if(x < .5){ x += x; ex -= 1; } } if(flip){ xy = 1. / xy; ey = -ey; } return ldexp(xy, ey); }