add a fast sqrt function for values in between 0 and 1

1 files changed, 22 insertions, 0 deletions

diff --git a/src/misc.c b/src/misc.c
index f705c74..003b170 100644
--- a/src/misc.c
+++ b/src/misc.c
@@ -881,3 +881,25 @@ double fast_acos(double x)
 
     return interp1d(x,xs,ys,LEN(xs));
 }
+
+/* Fast version of sqrt() for `x` between 0 and 1 which uses a lookup table
+ * computed on the first call. */
+double fast_sqrt(double x)
+{
+    size_t i;
+    static int initialized = 0;
+    static double xs[N_SQRT];
+    static double ys[N_SQRT];
+
+    if (!initialized) {
+        for (i = 0; i < LEN(xs); i++) {
+            xs[i] = 1.0*i/(LEN(xs)-1);
+            ys[i] = sqrt(xs[i]);
+        }
+        initialized = 1;
+    }
+
+    if (x > 1.0) return sqrt(x);
+
+    return interp1d(x,xs,ys,LEN(xs));
+}