1 files changed, 83 insertions, 0 deletions
diff --git a/mathfuncs_asin.h b/mathfuncs_asin.h
new file mode 100644
index 0000000..a8ffff4
--- /dev/null
+++ b/mathfuncs_asin.h
@@ -0,0 +1,83 @@
+// -*-C++-*-
+
+#ifndef MATHFUNCS_ASIN_H
+#define MATHFUNCS_ASIN_H
+
+#include "mathfuncs_base.h"
+
+#include <cassert>
+#include <cmath>
+
+
+
+namespace vecmathlib {
+  
+  template<typename realvec_t>
+  realvec_t mathfuncs<realvec_t>::vml_atan(realvec_t x)
+  {
+    // Handle negative values
+    realvec_t x0 = x;
+    x = fabs(x);
+    
+    // Reduce range using 1/x identity
+    assert(all(x >= RV(0.0)));
+    boolvec_t gt_one = x > RV(1.0);
+    x = ifthen(gt_one, rcp(x), x);
+    
+    // Reduce range again using half-angle formula; see
+    // <https://en.wikipedia.org/wiki/Inverse_trigonometric_functions>.
+    // This is necessary for good convergence below.
+    x = x / (RV(1.0) + sqrt(RV(1.0) + x*x));
+    
+    // Taylor expansion; see
+    // <https://en.wikipedia.org/wiki/Inverse_trigonometric_functions>.
+    assert(all(x >= RV(0.0) && x <= RV(0.5)));
+    int const nmax = 30;        // ???
+    realvec_t y = x / (RV(1.0) + x*x);
+    realvec_t x2 = x * y;
+    realvec_t r = y;
+    for (int n=3; n<nmax; n+=2) {
+      y *= RV(R(n-1) / R(n)) * x2;
+      r += y;
+    }
+    
+    // Undo second range reduction
+    r = RV(2.0) * r;
+    
+    // Undo range reduction
+    r = ifthen(gt_one, RV(M_PI_2) - r, r);
+    
+    // Handle negative values
+    r = copysign(r, x0);
+    
+    return r;
+  }
+  
+  
+  
+  template<typename realvec_t>
+  realvec_t mathfuncs<realvec_t>::vml_acos(realvec_t x)
+  {
+    // Handle negative values
+    boolvec_t is_negative = signbit(x);
+    x = fabs(x);
+    
+    realvec_t r = RV(2.0) * atan(sqrt(RV(1.0) - x*x) / (RV(1.0) + x));
+    
+    // Handle negative values
+    r = ifthen(is_negative, RV(M_PI) - r, r);
+    
+    return r;
+  }
+  
+  
+  
+  template<typename realvec_t>
+  realvec_t mathfuncs<realvec_t>::vml_asin(realvec_t x)
+  {
+    return RV(2.0) * atan(x / (RV(1.0) + sqrt(RV(1.0) - x*x)));
+  }
+  
+}; // namespace vecmathlib
+
+#endif  // #ifndef MATHFUNCS_ASIN_H