Logo ROOT   6.08/07
Reference Guide
sincos.h
Go to the documentation of this file.
1 /*
2  * sincos_common.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 #include "vdtcore_common.h"
28 #include <cmath>
29 #include <limits>
30 
31 #ifndef SINCOS_COMMON_H_
32 #define SINCOS_COMMON_H_
33 
34 namespace vdt{
35 
36 namespace details{
37 
38 // double precision constants
39 
40 const double DP1sc = 7.85398125648498535156E-1;
41 const double DP2sc = 3.77489470793079817668E-8;
42 const double DP3sc = 2.69515142907905952645E-15;
43 
44 const double C1sin = 1.58962301576546568060E-10;
45 const double C2sin =-2.50507477628578072866E-8;
46 const double C3sin = 2.75573136213857245213E-6;
47 const double C4sin =-1.98412698295895385996E-4;
48 const double C5sin = 8.33333333332211858878E-3;
49 const double C6sin =-1.66666666666666307295E-1;
50 
51 const double C1cos =-1.13585365213876817300E-11;
52 const double C2cos = 2.08757008419747316778E-9;
53 const double C3cos =-2.75573141792967388112E-7;
54 const double C4cos = 2.48015872888517045348E-5;
55 const double C5cos =-1.38888888888730564116E-3;
56 const double C6cos = 4.16666666666665929218E-2;
57 
58 const double DP1 = 7.853981554508209228515625E-1;
59 const double DP2 = 7.94662735614792836714E-9;
60 const double DP3 = 3.06161699786838294307E-17;
61 
62 // single precision constants
63 
64 const float DP1F = 0.78515625;
65 const float DP2F = 2.4187564849853515625e-4;
66 const float DP3F = 3.77489497744594108e-8;
67 
68 const float T24M1 = 16777215.;
69 
70 //------------------------------------------------------------------------------
71 
72 inline double get_sin_px(const double x){
73  double px=C1sin;
74  px *= x;
75  px += C2sin;
76  px *= x;
77  px += C3sin;
78  px *= x;
79  px += C4sin;
80  px *= x;
81  px += C5sin;
82  px *= x;
83  px += C6sin;
84  return px;
85 }
86 
87 //------------------------------------------------------------------------------
88 
89 inline double get_cos_px(const double x){
90  double px=C1cos;
91  px *= x;
92  px += C2cos;
93  px *= x;
94  px += C3cos;
95  px *= x;
96  px += C4cos;
97  px *= x;
98  px += C5cos;
99  px *= x;
100  px += C6cos;
101  return px;
102 }
103 
104 
105 //------------------------------------------------------------------------------
106 /// Reduce to 0 to 45
107 inline double reduce2quadrant(double x, int32_t& quad) {
108 
109  x = fabs(x);
110  quad = int (ONEOPIO4 * x); // always positive, so (int) == std::floor
111  quad = (quad+1) & (~1);
112  const double y = double (quad);
113  // Extended precision modular arithmetic
114  return ((x - y * DP1) - y * DP2) - y * DP3;
115  }
116 
117 //------------------------------------------------------------------------------
118 /// Sincos only for -45deg < x < 45deg
119 inline void fast_sincos_m45_45( const double z, double & s, double &c ) {
120 
121  double zz = z * z;
122  s = z + z * zz * get_sin_px(zz);
123  c = 1.0 - zz * .5 + zz * zz * get_cos_px(zz);
124  }
125 
126 
127 //------------------------------------------------------------------------------
128 
129 } // End namespace details
130 
131 /// Double precision sincos
132 inline void fast_sincos( const double xx, double & s, double &c ) {
133  // I have to use doubles to make it vectorise...
134 
135  int j;
136  double x = details::reduce2quadrant(xx,j);
137  const double signS = (j&4);
138 
139  j-=2;
140 
141  const double signC = (j&4);
142  const double poly = j&2;
143 
145 
146  //swap
147  if( poly==0 ) {
148  const double tmp = c;
149  c=s;
150  s=tmp;
151  }
152 
153  if(signC == 0.)
154  c = -c;
155  if(signS != 0.)
156  s = -s;
157  if (xx < 0.)
158  s = -s;
159 
160  }
161 
162 
163 // Single precision functions
164 
165 namespace details {
166 //------------------------------------------------------------------------------
167 /// Reduce to 0 to 45
168 inline float reduce2quadrant(float x, int & quad) {
169  /* make argument positive */
170  x = fabs(x);
171 
172  quad = int (ONEOPIO4F * x); /* integer part of x/PIO4 */
173 
174  quad = (quad+1) & (~1);
175  const float y = float(quad);
176  // quad &=4;
177  // Extended precision modular arithmetic
178  return ((x - y * DP1F) - y * DP2F) - y * DP3F;
179  }
180 
181 
182 //------------------------------------------------------------------------------
183 
184 
185 
186 /// Sincos only for -45deg < x < 45deg
187 inline void fast_sincosf_m45_45( const float x, float & s, float &c ) {
188 
189  float z = x * x;
190 
191  s = (((-1.9515295891E-4f * z
192  + 8.3321608736E-3f) * z
193  - 1.6666654611E-1f) * z * x)
194  + x;
195 
196  c = (( 2.443315711809948E-005f * z
197  - 1.388731625493765E-003f) * z
198  + 4.166664568298827E-002f) * z * z
199  - 0.5f * z + 1.0f;
200  }
201 
202 //------------------------------------------------------------------------------
203 
204 } // end details namespace
205 
206 /// Single precision sincos
207 inline void fast_sincosf( const float xx, float & s, float &c ) {
208 
209 
210  int j;
211  const float x = details::reduce2quadrant(xx,j);
212  int signS = (j&4);
213 
214  j-=2;
215 
216  const int signC = (j&4);
217  const int poly = j&2;
218 
219  float ls,lc;
221 
222  //swap
223  if( poly==0 ) {
224  const float tmp = lc;
225  lc=ls; ls=tmp;
226  }
227 
228  if(signC == 0) lc = -lc;
229  if(signS != 0) ls = -ls;
230  if (xx<0) ls = -ls;
231  c=lc;
232  s=ls;
233  }
234 
235 
236 } // end namespace vdt
237 
238 #endif
double reduce2quadrant(double x, int32_t &quad)
Reduce to 0 to 45.
Definition: sincos.h:107
double get_sin_px(const double x)
Definition: sincos.h:72
const double DP1sc
Definition: sincos.h:40
void fast_sincos(const double xx, double &s, double &c)
Double precision sincos.
Definition: sincos.h:132
return c
const float ONEOPIO4F
const double DP2
Definition: sincos.h:59
const double DP3sc
Definition: sincos.h:42
const double C2sin
Definition: sincos.h:45
const double C3cos
Definition: sincos.h:53
const double DP3
Definition: sincos.h:60
double get_cos_px(const double x)
Definition: sincos.h:89
const double C1sin
Definition: sincos.h:44
const double DP2sc
Definition: sincos.h:41
Double_t x[n]
Definition: legend1.C:17
void fast_sincos_m45_45(const double z, double &s, double &c)
Sincos only for -45deg < x < 45deg.
Definition: sincos.h:119
void fast_sincosf_m45_45(const float x, float &s, float &c)
Sincos only for -45deg < x < 45deg.
Definition: sincos.h:187
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
const float DP2F
Definition: sincos.h:65
const double C5sin
Definition: sincos.h:48
const double C1cos
Definition: sincos.h:51
Double_t E()
Definition: TMath.h:54
const double ONEOPIO4
const double C4cos
Definition: sincos.h:54
const double C5cos
Definition: sincos.h:55
const double DP1
Definition: sincos.h:58
const double C4sin
Definition: sincos.h:47
double f(double x)
Double_t y[n]
Definition: legend1.C:17
const double C6sin
Definition: sincos.h:49
void fast_sincosf(const float xx, float &s, float &c)
Single precision sincos.
Definition: sincos.h:207
Definition: asin.h:32
you should not use this method at all Int_t Int_t z
Definition: TRolke.cxx:630
const float DP1F
Definition: sincos.h:64
const float T24M1
Definition: sincos.h:68
const float DP3F
Definition: sincos.h:66
const double C3sin
Definition: sincos.h:46
const double C2cos
Definition: sincos.h:52
const double C6cos
Definition: sincos.h:56