Import sources from 2011-03-30 iso image

1 files changed, 76 insertions, 0 deletions

diff --git a/sys/src/ape/lib/ap/math/pow.c b/sys/src/ape/lib/ap/math/pow.c
new file mode 100755
index 000000000..9512a649b
--- /dev/null
+++ b/sys/src/ape/lib/ap/math/pow.c
@@ -0,0 +1,76 @@
+#include <math.h>
+#include <errno.h>
+#include <limits.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 > LONG_MAX){
+		if(x <= 0.){
+ zreturn:
+			if(x || flip)
+				errno = EDOM;
+			return 0.;
+		}
+		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;
+	}
+	errno = 0;
+	x = ldexp(xy, ey);
+	if(errno && ey < 0) {
+		errno = 0;
+		x = 0.;
+	}
+	return x;
+}