From 70a440e997b3b34c1b07f4110298cd8078d8f279 Mon Sep 17 00:00:00 2001 From: tlatorre Date: Mon, 29 Jul 2019 12:28:32 -0500 Subject: 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. --- src/misc.c | 32 +++++++++++++++++--------------- 1 file 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; } -- cgit