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
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
|
/*
* tan.h
* The basic idea is to exploit Pade polynomials.
* A lot of ideas were inspired by the cephes math library (by Stephen L. Moshier
* moshier@na-net.ornl.gov) as well as actual code.
* The Cephes library can be found here: http://www.netlib.org/cephes/
*
* Created on: Jun 23, 2012
* Author: Danilo Piparo, Thomas Hauth, Vincenzo Innocente
*/
/*
* VDT is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Lesser Public License for more details.
*
* You should have received a copy of the GNU Lesser Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef TAN_H_
#define TAN_H_
#include "vdtcore_common.h"
#include "sincos.h"
namespace vdt{
namespace details{
const double PX1tan=-1.30936939181383777646E4;
const double PX2tan=1.15351664838587416140E6;
const double PX3tan=-1.79565251976484877988E7;
const double QX1tan = 1.36812963470692954678E4;
const double QX2tan = -1.32089234440210967447E6;
const double QX3tan = 2.50083801823357915839E7;
const double QX4tan = -5.38695755929454629881E7;
const double DP1tan = 7.853981554508209228515625E-1;
const double DP2tan = 7.94662735614792836714E-9;
const double DP3tan = 3.06161699786838294307E-17;
const float DP1Ftan = 0.78515625;
const float DP2Ftan = 2.4187564849853515625e-4;
const float DP3Ftan = 3.77489497744594108e-8;
//------------------------------------------------------------------------------
/// Reduce to -45 to 45
inline double reduce2quadranttan(double x, int32_t& quad) {
x = fabs(x);
quad = int( ONEOPIO4 * x ); // always positive, so (int) == std::floor
quad = (quad+1) & (~1);
const double y = quad;
// Extended precision modular arithmetic
return ((x - y * DP1tan) - y * DP2tan) - y * DP3tan;
}
//------------------------------------------------------------------------------
/// Reduce to -45 to 45
inline float reduce2quadranttan(float x, int32_t& quad) {
x = fabs(x);
quad = int( ONEOPIO4F * x ); // always positive, so (int) == std::floor
quad = (quad+1) & (~1);
const float y = quad;
// Extended precision modular arithmetic
return ((x - y * DP1Ftan) - y * DP2Ftan) - y * DP3Ftan;
}
}
//------------------------------------------------------------------------------
/// Double precision tangent implementation
inline double fast_tan(double x){
const uint64_t sign_mask = details::getSignMask(x);
int32_t quad =0;
const double z=details::reduce2quadranttan(x,quad);
const double zz = z * z;
double res=z;
if( zz > 1.0e-14 ){
double px = details::PX1tan;
px *= zz;
px += details::PX2tan;
px *= zz;
px += details::PX3tan;
double qx=zz;
qx += details::QX1tan;
qx *=zz;
qx += details::QX2tan;
qx *=zz;
qx += details::QX3tan;
qx *=zz;
qx += details::QX4tan;
res = z + z * zz * px / qx;
}
// A no branching way to say: if j&2 res = -1/res. You can!!!
quad &=2;
quad >>=1;
const int32_t alt = quad^1;
// Avoid fpe generated by 1/0 if res is 0
const double zeroIfXNonZero = (x==0.);
res += zeroIfXNonZero;
res = quad * (-1./res) + alt * res; // one coeff is one and one is 0!
// Again, return 0 if res==0, the correct result otherwhise
return details::dpXORuint64(res,sign_mask) * (1.-zeroIfXNonZero);
}
// Single precision ------------------------------------------------------------
inline float fast_tanf(float x){
const uint32_t sign_mask = details::getSignMask(x);
int32_t quad =0;
const float z=details::reduce2quadranttan(x,quad);
const float zz = z * z;
float res=z;
if( zz > 1.0e-14f ){
res =
((((( 9.38540185543E-3f * zz
+ 3.11992232697E-3f) * zz
+ 2.44301354525E-2f) * zz
+ 5.34112807005E-2f) * zz
+ 1.33387994085E-1f) * zz
+ 3.33331568548E-1f) * zz * z
+ z;
}
// A no branching way to say: if j&2 res = -1/res. You can!!!
quad &=2;
quad >>=1;
const int32_t alt = quad^1;
// Avoid fpe generated by 1/0 if res is 0
const float zeroIfXNonZero = (x==0.f);
res += zeroIfXNonZero;
res = quad * (-1.f/res) + alt * res; // one coeff is one and one is 0!
return details::spXORuint32(res,sign_mask) * (1.f-zeroIfXNonZero);
}
//------------------------------------------------------------------------------
void tanv(const uint32_t size, double const * __restrict__ iarray, double* __restrict__ oarray);
void fast_tanv(const uint32_t size, double const * __restrict__ iarray, double* __restrict__ oarray);
void tanfv(const uint32_t size, float const * __restrict__ iarray, float* __restrict__ oarray);
void fast_tanfv(const uint32_t size, float const * __restrict__ iarray, float* __restrict__ oarray);
//------------------------------------------------------------------------------
} //vdt namespace
#endif /* TAN_H_ */
|
