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
|
// -*-C++-*-
#ifndef MATHFUNCS_LOG_H
#define MATHFUNCS_LOG_H
#include "mathfuncs_base.h"
#include <cmath>
namespace vecmathlib {
template <typename realvec_t>
realvec_t mathfuncs<realvec_t>::vml_log2(realvec_t x) {
// Algorithm inspired by SLEEF 2.80
// Rescale
intvec_t ilogb_x = ilogb(x * RV(M_SQRT2));
x = ldexp(x, -ilogb_x);
VML_ASSERT(all(x >= RV(M_SQRT1_2) && x <= RV(M_SQRT2)));
realvec_t y = (x - RV(1.0)) / (x + RV(1.0));
realvec_t y2 = y * y;
realvec_t r;
switch (sizeof(real_t)) {
case 4:
// float, error=7.09807175879142775648452461821e-8
r = RV(0.59723611417135718739797302426);
r = mad(r, y2, RV(0.961524413175528426101613434));
r = mad(r, y2, RV(2.88539097665498228703236701));
break;
case 8:
#ifdef VML_HAVE_FP_CONTRACT
// double, error=1.48294180185938512675770096324e-16
r = RV(0.243683403415639178527756320773);
r = mad(r, y2, RV(0.26136626803870009948502658));
r = mad(r, y2, RV(0.320619429891299265439389));
r = mad(r, y2, RV(0.4121983452028499242926));
r = mad(r, y2, RV(0.577078017761894161436));
r = mad(r, y2, RV(0.96179669392233355927));
r = mad(r, y2, RV(2.8853900817779295236));
#else
// double, error=2.1410114030383689267772704676e-14
r = RV(0.283751646449323373643963474845);
r = mad(r, y2, RV(0.31983138095551191299118812));
r = mad(r, y2, RV(0.412211603844146279666022));
r = mad(r, y2, RV(0.5770779098948940070516));
r = mad(r, y2, RV(0.961796694295973716912));
r = mad(r, y2, RV(2.885390081777562819196));
#endif
break;
default:
__builtin_unreachable();
}
r *= y;
// 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) {
// TODO: Check SLEEF 2.80 algorithm
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
|
