MeteoIODoc  MeteoIODoc-2.6.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MathOptim.h
Go to the documentation of this file.
1 /***********************************************************************************/
2 /* Copyright 2012 WSL Institute for Snow and Avalanche Research SLF-DAVOS */
3 /***********************************************************************************/
4 /* This file is part of MeteoIO.
5  MeteoIO is free software: you can redistribute it and/or modify
6  it under the terms of the GNU Lesser General Public License as published by
7  the Free Software Foundation, either version 3 of the License, or
8  (at your option) any later version.
9 
10  MeteoIO is distributed in the hope that it will be useful,
11  but WITHOUT ANY WARRANTY; without even the implied warranty of
12  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  GNU Lesser General Public License for more details.
14 
15  You should have received a copy of the GNU Lesser General Public License
16  along with MeteoIO. If not, see <http://www.gnu.org/licenses/>.
17 */
18 #ifndef MATHOPTIM_H
19 #define MATHOPTIM_H
20 
21 #include <stdint.h>
22 #include <cmath>
23 #include <string.h>
24 
25 //Quake3 fast 1/x² approximation
26 // For Magic Derivation see: Chris Lomont http://www.lomont.org/Math/Papers/2003/InvSqrt.pdf
27 // Credited to Greg Walsh.
28 // 32 Bit float magic number - for 64 bits doubles: 0x5fe6ec85e7de30da
29 #define SQRT_MAGIC_D 0x5f3759df
30 #define SQRT_MAGIC_F 0x5f375a86
31 
32 namespace mio {
33 
34 namespace Optim {
35 
44  inline long int round(const double& x) {
45  if (x>=0.) return static_cast<long int>( x+.5 );
46  else return static_cast<long int>( x-.5 );
47  }
48 
57  inline long int floor(const double& x) {
58  const long int xi = static_cast<long int>(x);
59  if (x >= 0 || static_cast<double>(xi) == x) return xi ;
60  else return xi - 1 ;
61  }
62 
71  inline long int ceil(const double& x) {
72  const long int xi = static_cast<long int>(x);
73  if (x <= 0 || static_cast<double>(xi) == x) return xi ;
74  else return xi + 1 ;
75  }
76 
77  inline double intPart(const double &x) {
78  double intpart;
79  modf(x, &intpart);
80  return intpart;
81  }
82 
83  inline double fracPart(const double &x) {
84  double intpart;
85  return modf(x, &intpart);
86  }
87 
88  #ifdef _MSC_VER
89  #pragma warning( push ) //for Visual C++
90  #pragma warning(disable:4244) //Visual C++ righhtfully complains... but this behavior is what we want!
91  #endif
92  //maximum relative error is <1.7% while computation time for sqrt is <1/4. At 0, returns a large number
93  //on a large scale interpolation test on TA, max relative error is 1e-6
94  inline float invSqrt(const float x) {
95  const float xhalf = 0.5f*x;
96 
97  union { // get bits for floating value
98  float x;
99  int i;
100  } u;
101  u.x = x;
102  u.i = SQRT_MAGIC_F - (u.i >> 1); // gives initial guess y0
103  return u.x*(1.5f - xhalf*u.x*u.x);// Newton step, repeating increases accuracy
104  }
105 
106  #ifdef __clang__
107  #pragma clang diagnostic push
108  #pragma clang diagnostic ignored "-Wdouble-promotion"
109  #endif
110  inline double invSqrt(const double x) {
111  const double xhalf = 0.5f*x;
112 
113  union { // get bits for floating value
114  float x;
115  int i;
116  } u;
117  u.x = static_cast<float>(x);
118  u.i = SQRT_MAGIC_D - (u.i >> 1); // gives initial guess y0
119  return u.x*(1.5f - xhalf*u.x*u.x);// Newton step, repeating increases accuracy
120  }
121  #ifdef __clang__
122  #pragma clang diagnostic pop
123  #endif
124 
125  #ifdef _MSC_VER
126  #pragma warning( pop ) //for Visual C++, restore previous warnings behavior
127  #endif
128 
129  inline float fastSqrt_Q3(const float x) {
130  return x * invSqrt(x);
131  }
132 
133  inline double fastSqrt_Q3(const double x) {
134  return x * invSqrt(x);
135  }
136 
137  inline double pow2(const double& val) {return (val*val);}
138  inline double pow3(const double& val) {return (val*val*val);}
139  inline double pow4(const double& val) {return (val*val*val*val);}
140 
141  //please do not use this method directly, call fastPow() instead!
142  inline double fastPowInternal(double a, double b) {
143  //see http://martin.ankerl.com/2012/01/25/optimized-approximative-pow-in-c-and-cpp/
144  // calculate approximation with fraction of the exponent
145  int e = (int) b;
146  union {
147  double d;
148  int x[2];
149  } u = { a };
150  u.x[1] = (int)((b - e) * (u.x[1] - 1072632447) + 1072632447);
151  u.x[0] = 0;
152 
153  // exponentiation by squaring with the exponent's integer part
154  // double r = u.d makes everything much slower, not sure why
155  double r = 1.0;
156  while (e) {
157  if (e & 1) {
158  r *= a;
159  }
160  a *= a;
161  e >>= 1;
162  }
163 
164  return r * u.d;
165  }
166 
179  inline double fastPow(double a, double b) {
180  if (b>0.) {
181  return fastPowInternal(a,b);
182  } else {
183  const double tmp = fastPowInternal(a,-b);
184  return 1./tmp;
185  }
186  }
187 
188  #ifdef __clang__
189  #pragma clang diagnostic push
190  #pragma clang diagnostic ignored "-Wundefined-reinterpret-cast"
191  #endif
192  //see http://metamerist.com/cbrt/cbrt.htm
193  template <int n> inline float nth_rootf(float x) {
194  const bool sgn = (x<0.f)? true : false;
195  if (sgn) x = -x;
196  const int ebits = 8;
197  const int fbits = 23;
198 
199  const int bias = (1 << (ebits-1))-1;
200  int& i = reinterpret_cast<int&>(x);
201  i = (i - (bias << fbits)) / n + (bias << fbits);
202 
203  if (sgn) return -x;
204  else return x;
205  }
206 
207  template <int n> inline double nth_rootd(double x) {
208  const bool sgn = (x<0.)? true : false;
209  if (sgn) x = -x;
210  const int ebits = 11;
211  const int fbits = 52;
212 
213  const int64_t bias = (1 << (ebits-1))-1;
214  int64_t& i = reinterpret_cast<int64_t&>(x);
215  i = (i - (bias << fbits)) / n + (bias << fbits);
216 
217  if (sgn) return -x;
218  else return x;
219  }
220  #ifdef __clang__
221  #pragma clang diagnostic pop
222  #endif
223 
236  inline double cbrt(double x) {
237  const double a = nth_rootd<3>(x);
238  const double a3 = a*a*a;
239  const double b = a * ( (a3 + x) + x) / ( a3 + (a3 + x) );
240  return b; //single iteration, otherwise set a=b and do it again
241  }
242 
253  inline double pow10(double x) {
254  const double a1 = 1.1499196;
255  const double a2 = 0.6774323;
256  const double a3 = 0.2080030;
257  const double a4 = 0.1268089;
258 
259  const double x2 = x*x;
260  const double tmp = 1. + a1*x + a2*x*x + a3*x*x2 + a4*x2*x2;
261  return tmp*tmp;
262  }
263 
264  template <typename T> T fastPow(T p, unsigned q) {
265  T r(1);
266 
267  while (q != 0) {
268  if (q % 2 == 1) { // q is odd
269  r *= p;
270  q--;
271  }
272  p *= p;
273  q /= 2;
274  }
275 
276  return r;
277  }
278 
289  inline double ln_1plusX(double x) {
290  const double a1 = 0.9974442;
291  const double a2 = -.4712839;
292  const double a3 = 0.2256685;
293  const double a4 = -.0587527;
294 
295  const double x2 = x*x;
296  return a1*x + a2*x2 + a3*x*x2 + a4*x2*x2;
297  }
298 
299 }
300 
301 } //end namespace
302 
303 #endif
long int floor(const double &x)
Optimized version of c++ floor() This version works with positive and negative numbers but does not c...
Definition: MathOptim.h:57
#define SQRT_MAGIC_D
Definition: MathOptim.h:29
double pow2(const double &val)
Definition: MathOptim.h:137
double pow4(const double &val)
Definition: MathOptim.h:139
double ln_1plusX(double x)
Optimized version of ln(1+x) This works for 0 <= x <= 1 and offers a theoritical precision of 5e-5 So...
Definition: MathOptim.h:289
double fastPow(double a, double b)
Optimized version of c++ pow() This version works with positive and negative exponents and handles ex...
Definition: MathOptim.h:179
double cbrt(double x)
Optimized version of cubic root This version is based on a single iteration Halley's method (see http...
Definition: MathOptim.h:236
#define SQRT_MAGIC_F
Definition: MathOptim.h:30
float invSqrt(const float x)
Definition: MathOptim.h:94
long int round(const double &x)
Optimized version of c++ round() This version works with positive and negative numbers but does not c...
Definition: MathOptim.h:44
double pow3(const double &val)
Definition: MathOptim.h:138
double fracPart(const double &x)
Definition: MathOptim.h:83
double intPart(const double &x)
Definition: MathOptim.h:77
double fastPowInternal(double a, double b)
Definition: MathOptim.h:142
float fastSqrt_Q3(const float x)
Definition: MathOptim.h:129
double nth_rootd(double x)
Definition: MathOptim.h:207
float nth_rootf(float x)
Definition: MathOptim.h:193
const double e
Definition: Meteoconst.h:66
double pow10(double x)
Optimized version of 10^x This works for 0 <= x <= 1 and offers a theoritical precision of 5e-5 Sourc...
Definition: MathOptim.h:253
long int ceil(const double &x)
Optimized version of c++ ceil() This version works with positive and negative numbers but does not co...
Definition: MathOptim.h:71