Logo ROOT   6.08/07
Reference Guide
tan.h
Go to the documentation of this file.
1 /*
2  * tan.h
3  * The basic idea is to exploit Pade polynomials.
4  * A lot of ideas were inspired by the cephes math library (by Stephen L. Moshier
5  * moshier@na-net.ornl.gov) as well as actual code.
6  * The Cephes library can be found here: http://www.netlib.org/cephes/
7  *
8  * Created on: Jun 23, 2012
9  * Author: Danilo Piparo, Thomas Hauth, Vincenzo Innocente
10  */
11 
12 /*
13  * VDT is free software: you can redistribute it and/or modify
14  * it under the terms of the GNU Lesser Public License as published by
15  * the Free Software Foundation, either version 3 of the License, or
16  * (at your option) any later version.
17  *
18  * This program is distributed in the hope that it will be useful,
19  * but WITHOUT ANY WARRANTY; without even the implied warranty of
20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21  * GNU Lesser Public License for more details.
22  *
23  * You should have received a copy of the GNU Lesser Public License
24  * along with this program. If not, see <http://www.gnu.org/licenses/>.
25  */
26 
27 #ifndef TAN_H_
28 #define TAN_H_
29 
30 #include "vdtcore_common.h"
31 #include "sincos.h"
32 
33 namespace vdt{
34 
35 
36 namespace details{
37 
38 const double PX1tan=-1.30936939181383777646E4;
39 const double PX2tan=1.15351664838587416140E6;
40 const double PX3tan=-1.79565251976484877988E7;
41 
42 const double QX1tan = 1.36812963470692954678E4;
43 const double QX2tan = -1.32089234440210967447E6;
44 const double QX3tan = 2.50083801823357915839E7;
45 const double QX4tan = -5.38695755929454629881E7;
46 
47 const double DP1tan = 7.853981554508209228515625E-1;
48 const double DP2tan = 7.94662735614792836714E-9;
49 const double DP3tan = 3.06161699786838294307E-17;
50 
51 const float DP1Ftan = 0.78515625;
52 const float DP2Ftan = 2.4187564849853515625e-4;
53 const float DP3Ftan = 3.77489497744594108e-8;
54 
55 
56 //------------------------------------------------------------------------------
57 /// Reduce to -45 to 45
58 inline double reduce2quadranttan(double x, int32_t& quad) {
59 
60  x = fabs(x);
61  quad = int( ONEOPIO4 * x ); // always positive, so (int) == std::floor
62  quad = (quad+1) & (~1);
63  const double y = quad;
64  // Extended precision modular arithmetic
65  return ((x - y * DP1tan) - y * DP2tan) - y * DP3tan;
66  }
67 
68 //------------------------------------------------------------------------------
69 /// Reduce to -45 to 45
70 inline float reduce2quadranttan(float x, int32_t& quad) {
71 
72  x = fabs(x);
73  quad = int( ONEOPIO4F * x ); // always positive, so (int) == std::floor
74  quad = (quad+1) & (~1);
75  const float y = quad;
76  // Extended precision modular arithmetic
77  return ((x - y * DP1Ftan) - y * DP2Ftan) - y * DP3Ftan;
78  }
79 
80 }
81 
82 //------------------------------------------------------------------------------
83 /// Double precision tangent implementation
84 inline double fast_tan(double x){
85 
86  const uint64_t sign_mask = details::getSignMask(x);
87 
88  int32_t quad =0;
89  const double z=details::reduce2quadranttan(x,quad);
90 
91  const double zz = z * z;
92 
93  double res=z;
94 
95  if( zz > 1.0e-14 ){
96  double px = details::PX1tan;
97  px *= zz;
98  px += details::PX2tan;
99  px *= zz;
100  px += details::PX3tan;
101 
102  double qx=zz;
103  qx += details::QX1tan;
104  qx *=zz;
105  qx += details::QX2tan;
106  qx *=zz;
107  qx += details::QX3tan;
108  qx *=zz;
109  qx += details::QX4tan;
110 
111  res = z + z * zz * px / qx;
112  }
113 
114  // A no branching way to say: if j&2 res = -1/res. You can!!!
115  quad &=2;
116  quad >>=1;
117  const int32_t alt = quad^1;
118  // Avoid fpe generated by 1/0 if res is 0
119  const double zeroIfXNonZero = (x==0.);
120  res += zeroIfXNonZero;
121  res = quad * (-1./res) + alt * res; // one coeff is one and one is 0!
122 
123  // Again, return 0 if res==0, the correct result otherwhise
124  return details::dpXORuint64(res,sign_mask) * (1.-zeroIfXNonZero);
125 
126 }
127 
128 // Single precision ------------------------------------------------------------
129 
130 inline float fast_tanf(float x){
131  const uint32_t sign_mask = details::getSignMask(x);
132 
133  int32_t quad =0;
134  const float z=details::reduce2quadranttan(x,quad);
135 
136  const float zz = z * z;
137 
138  float res=z;
139 
140  if( zz > 1.0e-14f ){
141  res =
142  ((((( 9.38540185543E-3f * zz
143  + 3.11992232697E-3f) * zz
144  + 2.44301354525E-2f) * zz
145  + 5.34112807005E-2f) * zz
146  + 1.33387994085E-1f) * zz
147  + 3.33331568548E-1f) * zz * z
148  + z;
149  }
150 
151  // A no branching way to say: if j&2 res = -1/res. You can!!!
152  quad &=2;
153  quad >>=1;
154  const int32_t alt = quad^1;
155  // Avoid fpe generated by 1/0 if res is 0
156  const float zeroIfXNonZero = (x==0.f);
157  res += zeroIfXNonZero;
158  res = quad * (-1.f/res) + alt * res; // one coeff is one and one is 0!
159 
160  return details::spXORuint32(res,sign_mask) * (1.f-zeroIfXNonZero);
161 
162 }
163 
164 //------------------------------------------------------------------------------
165 
166 // void tanv(const uint32_t size, double const * __restrict__ iarray, double* __restrict__ oarray);
167 // void fast_tanv(const uint32_t size, double const * __restrict__ iarray, double* __restrict__ oarray);
168 // void tanfv(const uint32_t size, float const * __restrict__ iarray, float* __restrict__ oarray);
169 // void fast_tanfv(const uint32_t size, float const * __restrict__ iarray, float* __restrict__ oarray);
170 
171 //------------------------------------------------------------------------------
172 
173 } //vdt namespace
174 
175 
176 #endif /* COS_H_ */
const float DP1Ftan
Definition: tan.h:51
const double QX3tan
Definition: tan.h:44
float fast_tanf(float x)
Definition: tan.h:130
const float DP2Ftan
Definition: tan.h:52
const double DP2tan
Definition: tan.h:48
const float ONEOPIO4F
double reduce2quadranttan(double x, int32_t &quad)
Reduce to -45 to 45.
Definition: tan.h:58
const double DP1tan
Definition: tan.h:47
const double QX1tan
Definition: tan.h:42
Double_t x[n]
Definition: legend1.C:17
double dpXORuint64(const double x, const uint64_t i)
Makes a XOR of a double and a unsigned long long.
uint64_t getSignMask(const double x)
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
const double QX4tan
Definition: tan.h:45
Double_t E()
Definition: TMath.h:54
const double ONEOPIO4
const double DP3tan
Definition: tan.h:49
double f(double x)
const double PX1tan
Definition: tan.h:38
float spXORuint32(const float x, const uint32_t i)
Makes an OR of a float and a unsigned long.
Double_t y[n]
Definition: legend1.C:17
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
Definition: asin.h:32
you should not use this method at all Int_t Int_t z
Definition: TRolke.cxx:630
const double PX3tan
Definition: tan.h:40
const float DP3Ftan
Definition: tan.h:53
double fast_tan(double x)
Double precision tangent implementation.
Definition: tan.h:84
const double QX2tan
Definition: tan.h:43
const double PX2tan
Definition: tan.h:39