summaryrefslogtreecommitdiff
path: root/sqrt.c
diff options
context:
space:
mode:
Diffstat (limited to 'sqrt.c')
-rw-r--r--sqrt.c71
1 files changed, 71 insertions, 0 deletions
diff --git a/sqrt.c b/sqrt.c
new file mode 100644
index 0000000..5a78693
--- /dev/null
+++ b/sqrt.c
@@ -0,0 +1,71 @@
+#include <u.h>
+#include <libc.h>
+
+int iters;
+
+/*
+ * Heron's method to compute the √
+ *
+ * iteratively do
+ * x1 = ½(x0 + n/x0)
+ * since
+ * lim M→∞ (xM) = √n
+ */
+//double
+//√(double n)
+//{
+// int i;
+// double x;
+//
+// x = 2;
+// for(i = 0; i < iters; i++)
+// x = 0.5*(x + n/x);
+// return x;
+//}
+
+double
+√(double n)
+{
+ double x0, x;
+
+ if(n == 0)
+ return 0;
+
+ x0 = -1;
+ x = n > 1? n/2: 1;
+ while(x0 != x){
+ x0 = x;
+ x = 0.5*(x0 + n/x0);
+ iters++;
+ }
+ return x;
+}
+
+void
+usage(void)
+{
+ fprint(2, "usage: %s number [prec]\n", argv0);
+ exits("usage");
+}
+
+void
+main(int argc, char *argv[])
+{
+ int prec;
+ double n;
+
+ prec = 10;
+ ARGBEGIN{
+ default: usage();
+ }ARGEND
+ if(argc < 1)
+ usage();
+
+ n = strtod(argv[0], nil);
+ if(n < 0)
+ sysfatal("too complex");
+ if(argc > 2)
+ prec = strtoul(argv[1], nil, 10);
+ print("√%g = %.*f (took %d iterations)\n", n, prec, √(n), iters);
+ exits(nil);
+}