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
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
|
// -*-C++-*-
#ifndef MATHFUNCS_LOG_H
#define MATHFUNCS_LOG_H
#include "mathfuncs_base.h"
#include <cassert>
#include <cmath>
namespace vecmathlib {
template<typename realvec_t>
realvec_t mathfuncs<realvec_t>::vml_log2(realvec_t x)
{
// Rescale
assert(all(x > RV(0.0)));
intvec_t ilogb_x = ilogb(x);
x = scalbn(x, -ilogb_x);
assert(all(x >= RV(1.0) && x < RV(2.0)));
// Approximate
// for |x|>0.01: (*) log(x) = Sum[n=1,nmax,n%2==1] 2/n ((x-1) / (x+1))^n
// else: (**) log(x) = Sum[n=1,nmax] (-1)^(n+1) 1/n (x-1)^n
assert(all(x >= RV(1.0) && x < RV(2.0)));
// nmax max_error of (*)
// 5 5.9e-5
// 7 1.3e-6
// 9 2.9e-8
// 15 4.4e-13
// 17 1.1e-14
// 19 3.0e-16
int const nmax = sizeof(real_t)==4 ? 9 : 17;
x *= RV(M_SQRT1_2); // shift range to increase accuracy
realvec_t xm1 = x - RV(1.0);
boolvec_t near1 = fabs(xm1) < RV(0.0001); // epsilon^(1/niters)
// for (*)
realvec_t xm1_xp1 = xm1 / (x + RV(1.0));
realvec_t xm1_xp1_2 = xm1_xp1 * xm1_xp1;
// for (**)
realvec_t mxm1 = - xm1;
realvec_t y = ifthen(near1, xm1, RV(2.0) * xm1_xp1);
realvec_t yf = ifthen(near1, mxm1, xm1_xp1_2);
y *= RV(M_LOG2E);
realvec_t r = y;
for (int n=3, nn=2; n<nmax; n+=2, ++nn) {
y *= yf;
r += y * ifthen(near1, RV(R(1.0) / R(nn)), RV(R(1.0) / R(n)));
}
r += RV(0.5); // correct result for range shift
// Undo rescaling
r += convert_float(ilogb_x);
return r;
}
template<typename realvec_t>
inline
realvec_t mathfuncs<realvec_t>::vml_log(realvec_t x)
{
return log2(x) * RV(M_LN2);
}
template<typename realvec_t>
inline
realvec_t mathfuncs<realvec_t>::vml_log10(realvec_t x)
{
return log(x) * RV(M_LOG10E);
}
template<typename realvec_t>
inline
realvec_t mathfuncs<realvec_t>::vml_log1p(realvec_t x)
{
return log(RV(1.0) + x);
#if 0
// Goldberg, theorem 4
realvec_t x1 = RV(1.0) + x;
x1.barrier();
return ifthen(x1 == x, x, x * log(x1) / (x1 - RV(1.0)));
#endif
}
}; // namespace vecmathlib
#endif // #ifndef MATHFUNCS_LOG_H
|
