diff options
Diffstat (limited to 'src/base/math/pow.c')
-rw-r--r-- | src/base/math/pow.c | 69 |
1 files changed, 69 insertions, 0 deletions
diff --git a/src/base/math/pow.c b/src/base/math/pow.c new file mode 100644 index 0000000..912a4e1 --- /dev/null +++ b/src/base/math/pow.c @@ -0,0 +1,69 @@ +#include <u.h> +#include <base.h> + +double +math·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 = math·modf(y, &ye); + if(y1 != 0.0){ + if(x <= 0.) + goto zreturn; + if(y1 > 0.5) { + y1 -= 1.; + ye += 1.; + } + xy = math·exp(y1 * math·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 math·NaN(); + } + if(flip){ + if(y == .5) + return 1./math·sqrt(x); + y = -y; + }else if(y == .5) + return math·sqrt(x); + return math·exp(y * math·log(x)); + } + x = math·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 math·ldexp(xy, ey); +} |