sqrt by Heron's method.

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);
+}