From e5888a1ffdae813d7575f5fb02275c6bb07e5199 Mon Sep 17 00:00:00 2001
From: Taru Karttunen <taruti@taruti.net>
Date: Wed, 30 Mar 2011 15:46:40 +0300
Subject: Import sources from 2011-03-30 iso image

---
 sys/src/ape/lib/ap/plan9/sqrt.c | 56 +++++++++++++++++++++++++++++++++++++++++
 1 file changed, 56 insertions(+)
 create mode 100755 sys/src/ape/lib/ap/plan9/sqrt.c

(limited to 'sys/src/ape/lib/ap/plan9/sqrt.c')

diff --git a/sys/src/ape/lib/ap/plan9/sqrt.c b/sys/src/ape/lib/ap/plan9/sqrt.c
new file mode 100755
index 000000000..98974360d
--- /dev/null
+++ b/sys/src/ape/lib/ap/plan9/sqrt.c
@@ -0,0 +1,56 @@
+/*
+	sqrt returns the square root of its floating
+	point argument. Newton's method.
+
+	calls frexp
+*/
+
+#include <math.h>
+#include <errno.h>
+#define _RESEARCH_SOURCE
+#include <float.h>
+
+double
+sqrt(double arg)
+{
+	double x, temp;
+	int exp, i;
+
+	if(isInf(arg, 1))
+		return arg;
+	if(arg <= 0) {
+		if(arg < 0)
+			errno = EDOM;
+		return 0;
+	}
+	x = frexp(arg, &exp);
+	while(x < 0.5) {
+		x *= 2;
+		exp--;
+	}
+	/*
+	 * NOTE
+	 * this wont work on 1's comp
+	 */
+	if(exp & 1) {
+		x *= 2;
+		exp--;
+	}
+	temp = 0.5 * (1.0+x);
+
+	while(exp > 60) {
+		temp *= (1L<<30);
+		exp -= 60;
+	}
+	while(exp < -60) {
+		temp /= (1L<<30);
+		exp += 60;
+	}
+	if(exp >= 0)
+		temp *= 1L << (exp/2);
+	else
+		temp /= 1L << (-exp/2);
+	for(i=0; i<=4; i++)
+		temp = 0.5*(temp + arg/temp);
+	return temp;
+}
-- 
cgit v1.2.3