1 files changed, 37 insertions, 0 deletions
diff --git a/heron-expanded.c b/heron-expanded.c
new file mode 100644
index 0000000..bd5d77f
--- /dev/null
+++ b/heron-expanded.c
@@ -0,0 +1,37 @@
+#include "c23-fallback.h"
+#include <stdio.h>
+
+int main(int argc, char* argv[argc+1]) {
+  constexpr double epsP1 = 1.0 + 1E-9;
+  constexpr double epsM1 = 1.0 - 1E-9;
+  for (int i = 1; i < argc; ++i) {
+    /* Read one command line argument as a double. */
+    double const a = strtod(argv[i], nullptr);
+    /* Compute some local constants. */
+    double fact;
+    double alow;
+    double ahig;
+    if (a < 1.0) {
+      fact =  2.0;
+      alow =  a;
+      ahig =  1.0;
+    } else {
+      fact =  0.5;
+      alow =  -a;
+      ahig =  -1.0;
+    }
+    /* A first low quality estimate for the inverse. */
+    double x = 1.0;
+    /* Adapt x until it has the right magnitude. */
+    while (alow*x < ahig) x *= fact;
+    /* We are close, correct with the Heron factor. */
+    for (double prod = a*x;
+         ((prod < epsM1) || (epsP1 < prod));
+         prod = a*x) {
+      x *= (2.0 - prod);
+    }
+    printf("heron: a=%.5e,\tx=%.5e,\ta*x=%.12f\n",
+           a, x, a*x);
+  }
+  return EXIT_SUCCESS;
+}