summaryrefslogtreecommitdiff
path: root/sys/src/libc/sparc/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/sparc/sqrt.c
Import sources from 2011-03-30 iso image
Diffstat (limited to 'sys/src/libc/sparc/sqrt.c')
-rwxr-xr-xsys/src/libc/sparc/sqrt.c103
1 files changed, 103 insertions, 0 deletions
diff --git a/sys/src/libc/sparc/sqrt.c b/sys/src/libc/sparc/sqrt.c
new file mode 100755
index 000000000..fa27c35ef
--- /dev/null
+++ b/sys/src/libc/sparc/sqrt.c
@@ -0,0 +1,103 @@
+#include <u.h>
+#include <libc.h>
+
+static long sqtab[64] =
+{
+ 0x6cdb2, 0x726d4, 0x77ea3, 0x7d52f, 0x82a85, 0x87eb1, 0x8d1c0, 0x923bd,
+ 0x974b2, 0x9c4a8, 0xa13a9, 0xa61be, 0xaaeee, 0xafb41, 0xb46bf, 0xb916e,
+ 0xbdb55, 0xc247a, 0xc6ce3, 0xcb495, 0xcfb95, 0xd41ea, 0xd8796, 0xdcca0,
+ 0xe110c, 0xe54dd, 0xe9818, 0xedac0, 0xf1cd9, 0xf5e67, 0xf9f6e, 0xfdfef,
+ 0x01fe0, 0x05ee6, 0x09cfd, 0x0da30, 0x11687, 0x1520c, 0x18cc8, 0x1c6c1,
+ 0x20000, 0x2388a, 0x27068, 0x2a79e, 0x2de32, 0x3142b, 0x3498c, 0x37e5b,
+ 0x3b29d, 0x3e655, 0x41989, 0x44c3b, 0x47e70, 0x4b02b, 0x4e16f, 0x51241,
+ 0x542a2, 0x57296, 0x5a220, 0x5d142, 0x60000, 0x62e5a, 0x65c55, 0x689f2,
+};
+
+double
+sqrt(double arg)
+{
+ int e, ms;
+ double a, t;
+ union
+ {
+ double d;
+ struct
+ {
+ long ms;
+ long ls;
+ };
+ } u;
+
+ u.d = arg;
+ ms = u.ms;
+
+ /*
+ * sign extend the mantissa with
+ * exponent. result should be > 0 for
+ * normal case.
+ */
+ e = ms >> 20;
+ if(e <= 0) {
+ if(e == 0)
+ return 0;
+ return NaN();
+ }
+
+ /*
+ * pick up arg/4 by adjusting exponent
+ */
+ u.ms = ms - (2 << 20);
+ a = u.d;
+
+ /*
+ * use 5 bits of mantissa and 1 bit
+ * of exponent to form table index.
+ * insert exponent/2 - 1.
+ */
+ e = (((e - 1023) >> 1) + 1022) << 20;
+ u.ms = *(long*)((char*)sqtab + ((ms >> 13) & 0xfc)) | e;
+ u.ls = 0;
+
+ /*
+ * three laps of newton
+ */
+ e = 1 << 20;
+ t = u.d;
+ u.d = t + a/t;
+ u.ms -= e; /* u.d /= 2; */
+ t = u.d;
+ u.d = t + a/t;
+ u.ms -= e; /* u.d /= 2; */
+ t = u.d;
+
+ return t + a/t;
+}
+
+/*
+ * this is the program that generated the table.
+ * it calls sqrt by some other means.
+ *
+ * void
+ * main(void)
+ * {
+ * int i;
+ * union U
+ * {
+ * double d;
+ * struct
+ * {
+ * long ms;
+ * long ls;
+ * };
+ * } u;
+ *
+ * for(i=0; i<64; i++) {
+ * u.ms = (i<<15) | 0x3fe04000;
+ * u.ls = 0;
+ * u.d = sqrt(u.d);
+ * print(" 0x%.5lux,", u.ms & 0xfffff);
+ * }
+ * print("\n");
+ * exits(0);
+ * }
+ */