ref: 21d0ffbbf001fd209208e33503890ad6847105db
parent: 738f91837bffa7e3ff9ac9323dd91cf9453cde77
author: rodri <rgl@antares-labs.eu>
date: Sat Jan 6 11:15:33 EST 2024
sqrt by Heron's method.
--- /dev/null
+++ b/sqrt.c
@@ -1,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);
+}