add a faster version of fast_acos()

This commit updates fast_acos() to use an algorithm I found here:
https://web.archive.org/web/20161223122122/http://http.developer.nvidia.com:80/Cg/acos.html
which is faster than using a lookup table.

1 files changed, 17 insertions, 15 deletions

diff --git a/src/misc.c b/src/misc.c
index f705c74..242fa50 100644
--- a/src/misc.c
+++ b/src/misc.c
@@ -863,21 +863,23 @@ void get_dir(double *dir, double theta, double phi)
     dir[2] = cos_theta;
 }
 
-/* Fast version of acos() which uses a lookup table computed on the first call. */
+/* Fast version of acos() from
+ * https://web.archive.org/web/20161223122122/http://http.developer.nvidia.com:80/Cg/acos.html.
+ *
+ * Note: I'm not sure if this correctly returns nan for `x` outside of [-1,1].
+ * */
 double fast_acos(double x)
 {
-    size_t i;
-    static int initialized = 0;
-    static double xs[N_ACOS];
-    static double ys[N_ACOS];
-
-    if (!initialized) {
-        for (i = 0; i < LEN(xs); i++) {
-            xs[i] = -1.0 + 2.0*i/(LEN(xs)-1);
-            ys[i] = acos(xs[i]);
-        }
-        initialized = 1;
-    }
-
-    return interp1d(x,xs,ys,LEN(xs));
+    double negate = (double) (x < 0);
+    x = fabs(x);
+    double ret = -0.0187293;
+    ret = ret * x;
+    ret = ret + 0.0742610;
+    ret = ret * x;
+    ret = ret - 0.2121144;
+    ret = ret * x;
+    ret = ret + 1.5707288;
+    ret = ret * sqrt(1.0-x);
+    ret = ret - 2 * negate * ret;
+    return negate * 3.14159265358979 + ret;
 }