blob: a8ffff48791d204a12b4a4f179874f5f66596512 (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
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
|
