aboutsummaryrefslogtreecommitdiff
path: root/src/base/math/pow.c
diff options
context:
space:
mode:
authorNicholas Noll <nbnoll@eml.cc>2021-12-05 16:16:21 -0800
committerNicholas Noll <nbnoll@eml.cc>2021-12-05 16:16:21 -0800
commit07e77936d535e58b0aeb4f2a11400c1050556739 (patch)
treefb50fad6436ecbc159505bee73a92cd289b99a07 /src/base/math/pow.c
parentb48327d357e0818d1a6ae2a064cfa7d1567e1242 (diff)
Feat: added math library
Used Plan9's libc as starting point. This cleans up dangling references due to loss of libc.
Diffstat (limited to 'src/base/math/pow.c')
-rw-r--r--src/base/math/pow.c69
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);
+}