summaryrefslogtreecommitdiff
path: root/sys/src/libc/port/sqrt.c
diff options
context:
space:
mode:
authorTaru Karttunen <taruti@taruti.net>2011-03-30 15:46:40 +0300
committerTaru Karttunen <taruti@taruti.net>2011-03-30 15:46:40 +0300
commite5888a1ffdae813d7575f5fb02275c6bb07e5199 (patch)
treed8d51eac403f07814b9e936eed0c9a79195e2450 /sys/src/libc/port/sqrt.c
Import sources from 2011-03-30 iso image
Diffstat (limited to 'sys/src/libc/port/sqrt.c')
-rwxr-xr-xsys/src/libc/port/sqrt.c54
1 files changed, 54 insertions, 0 deletions
diff --git a/sys/src/libc/port/sqrt.c b/sys/src/libc/port/sqrt.c
new file mode 100755
index 000000000..a92efac88
--- /dev/null
+++ b/sys/src/libc/port/sqrt.c
@@ -0,0 +1,54 @@
+/*
+ sqrt returns the square root of its floating
+ point argument. Newton's method.
+
+ calls frexp
+*/
+
+#include <u.h>
+#include <libc.h>
+
+double
+sqrt(double arg)
+{
+ double x, temp;
+ int exp, i;
+
+ if(arg <= 0) {
+ if(arg < 0)
+ return NaN();
+ return 0;
+ }
+ if(isInf(arg, 1))
+ return arg;
+ 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;
+}