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