CLAM-Development  1.4.0
CLAM_Math.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2004 MUSIC TECHNOLOGY GROUP (MTG)
3  * UNIVERSITAT POMPEU FABRA
4  *
5  *
6  * This program is free software; you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation; either version 2 of the License, or
9  * (at your option) any later version.
10  *
11  * This program 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 General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with this program; if not, write to the Free Software
18  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19  *
20  */
21 
22 #ifndef __CLAM_MATH__
23 #define __CLAM_MATH__
24 
25 #include <cmath>
26 #include "DataTypes.hxx"
27 #include "FastRounding.hxx"
28 
29 
30 //The following constants are defined also in OSDefines but only for windows and using the #define preprocessor
31 //directive. It is much better to use const float declarations
32 const float PI_ = 3.1415926535897932384626433832795028841972; /* pi */
33 const float ONE_OVER_PI = (0.3183098861837906661338147750939f);
34 const float TWOPI = (6.2831853071795864769252867665590057683943f); /* 2*pi */
35 const float ONE_OVER_TWOPI = (0.15915494309189535682609381638f);
36 const float PI_2 = (1.5707963267948966192313216916397514420986f); /* pi/2 */
37 const float TWO_OVER_PI = (0.636619772367581332267629550188f);
38 const float LN2 = (0.6931471805599453094172321214581765680755f); /* ln(2) */
39 const float ONE_OVER_LN2 = (1.44269504088896333066907387547f);
40 const float LN10 = (2.3025850929940456840179914546843642076011f); /* ln(10) */
41 const float ONE_OVER_LN10 = (0.43429448190325177635683940025f);
44 const long LONG_OFFSET = 4096L;
45 const float FLOAT_OFFSET = 4096.0;
46 const float HUGE_ = 1.0e8;
47 const float ROOT2 = (1.4142135623730950488016887242096980785697f); /* sqrt(2) */
48 
50 inline float CLAM_sin(register float x)
51 {
52 #ifndef CLAM_OPTIMIZE
53  return (float) sin((double)x);
54 #else
55  x *= ONE_OVER_PI;
56  register float accumulator, xPower, xSquared;
57  register long evenIntPart = ((long)(0.5f*x + 1024.5) - 1024)<<1;
58  x -= (float)evenIntPart;
59  xSquared = x*x;
60  accumulator = 3.14159265358979f*x;
61  xPower = xSquared*x;
62  accumulator += -5.16731953364340f*xPower;
63  xPower *= xSquared;
64  accumulator += 2.54620566822659f*xPower;
65  xPower *= xSquared;
66  accumulator += -0.586027023087261f*xPower;
67  xPower *= xSquared;
68  accumulator += 0.06554823491427f*xPower;
69  return accumulator;
70 #endif
71 }
72 
73 inline float CLAM_cos(register float x)
74  {
75 #ifndef CLAM_OPTIMIZE
76  return (float) cos((double)x);
77 #else
78  x *= ONE_OVER_PI;
79  register float accumulator, xPower, xSquared;
80 
81  register long evenIntPart = ((long)(0.5f*x + 1024.5f) - 1024)<<1;
82  x -= (float)evenIntPart;
83 
84  xSquared = x*x;
85  accumulator = 1.57079632679490f*x; /* series for sin(PI/2*x) */
86  xPower = xSquared*x;
87  accumulator += -0.64596406188166f*xPower;
88  xPower *= xSquared;
89  accumulator += 0.07969158490912f*xPower;
90  xPower *= xSquared;
91  accumulator += -0.00467687997706f*xPower;
92  xPower *= xSquared;
93  accumulator += 0.00015303015470f*xPower;
94  return 1.0f - 2.0f*accumulator*accumulator; /* cos(w) = 1 - 2*(sin(w/2))^2 */
95 #endif
96  }
97 
98 inline float CLAM_atan(register float x)
99  {
100 #ifndef CLAM_OPTIMIZE
101  return (float) atan((double)x);
102 #else
103  register float accumulator, xPower, xSquared, offset;
104 
105  offset = 0.0f;
106 
107  if (x < -1.0f)
108  {
109  offset = -PI_2;
110  x = -1.0f/x;
111  }
112  else if (x > 1.0f)
113  {
114  offset = PI_2;
115  x = -1.0f/x;
116  }
117  xSquared = x*x;
118  accumulator = 1.0f;
119  xPower = xSquared;
120  accumulator += 0.33288950512027f*xPower;
121  xPower *= xSquared;
122  accumulator += -0.08467922817644f*xPower;
123  xPower *= xSquared;
124  accumulator += 0.03252232640125f*xPower;
125  xPower *= xSquared;
126  accumulator += -0.00749305860992f*xPower;
127 
128  return offset + x/accumulator;
129 #endif
130 }
131 
132 inline float CLAM_atan2(float Imag, float Real)
133  {
134 #ifndef CLAM_OPTIMIZE
135  return (float) atan2((double)Imag, (double)Real);
136 #else
137  if(Real==0 && Imag==0) return 0.f;
138  register float accumulator, xPower, xSquared, offset, x;
139 
140  if (Imag > 0.0f)
141  {
142  if (Imag <= -Real)
143  {
144  offset = PI_;
145  x = Imag/Real;
146  }
147  else if (Imag > Real)
148  {
149  offset = PI_2;
150  x = -Real/Imag;
151  }
152  else
153  {
154  offset = 0.0f;
155  x = Imag/Real;
156  }
157  }
158  else
159  {
160  if (Imag >= Real)
161  {
162  offset = -PI_;
163  x = Imag/Real;
164  }
165  else if (Imag < -Real)
166  {
167  offset = -PI_2;
168  x = -Real/Imag;
169  }
170  else
171  {
172  offset = 0.0f;
173  x = Imag/Real;
174  }
175  }
176 
177  xSquared = x*x;
178  accumulator = 1.0f;
179  xPower = xSquared;
180  accumulator += 0.33288950512027f*xPower;
181  xPower *= xSquared;
182  accumulator += -0.08467922817644f*xPower;
183  xPower *= xSquared;
184  accumulator += 0.03252232640125f*xPower;
185  xPower *= xSquared;
186  accumulator += -0.00749305860992f*xPower;
187 
188  return offset + x/accumulator;
189 #endif
190 }
191 
192 inline float CLAM_exp2(register float x)
193 {
194 #ifndef CLAM_OPTIMIZE
195  return (float) exp(LN2*(double)x);
196 #else
197  if (x >= -127.0f)
198  {
199  register float accumulator, xPower;
200  register union {float f; long i;} xBits;
201 
202  xBits.i = (long)(x + FLOAT_OFFSET) - LONG_OFFSET; /* integer part */
203  x -= (float)(xBits.i); /* fractional part */
204 
205  accumulator = 1.0f + 0.69303212081966f*x;
206  xPower = x*x;
207  accumulator += 0.24137976293709f*xPower;
208  xPower *= x;
209  accumulator += 0.05203236900844f*xPower;
210  xPower *= x;
211  accumulator += 0.01355574723481f*xPower;
212 
213  xBits.i += 127; /* bias integer part */
214  xBits.i <<= 23; /* move biased int part into exponent bits */
215 
216  return accumulator * xBits.f;
217  }
218  else
219  {
220  return 0.0f;
221  }
222 #endif
223 }
224 
225 inline float CLAM_log2(register float x)
226 {
227 #ifndef CLAM_OPTIMIZE
228  return (float) (ONE_OVER_LN2*log((double)x));
229 #else
230  if (x > 5.877471754e-39f)
231  {
232  register float accumulator, xPower;
233  register long intPart;
234 
235  register union {float f; long i;} xBits;
236 
237  xBits.f = x;
238 
239  intPart = ((xBits.i)>>23); /* get biased exponent */
240  intPart -= 127; /* unbias it */
241 
242  x = (float)(xBits.i & 0x007FFFFF); /* mask off exponent leaving 0x800000*(mantissa - 1) */
243  x *= 1.192092895507812e-07f; /* divide by 0x800000 */
244 
245  accumulator = 1.44254494359510f*x;
246  xPower = x*x;
247  accumulator += -0.71814525675041f*xPower;
248  xPower *= x;
249  accumulator += 0.45754919692582f*xPower;
250  xPower *= x;
251  accumulator += -0.27790534462866f*xPower;
252  xPower *= x;
253  accumulator += 0.12179791068782f*xPower;
254  xPower *= x;
255  accumulator += -0.02584144982967f*xPower;
256 
257  return accumulator + (float)intPart;
258  }
259  else
260  {
261  return -HUGE_;
262  }
263 #endif
264 }
265 
266 inline float CLAM_pow(float x, float y)
267 {
268 #ifndef CLAM_OPTIMIZE
269  return (float) pow((double)x, (double)y);
270 #else
271  return CLAM_exp2(y*CLAM_log2(x));
272 #endif
273 }
274 
275 inline float CLAM_sqrt(register float x)
276  {
277 #ifndef CLAM_OPTIMIZE
278  return (float) sqrt((double)x);
279 #else
280  if (x > 5.877471754e-39f)
281  {
282  register float accumulator, xPower;
283  register long intPart;
284  register union {float f; long i;} xBits;
285 
286  xBits.f = x;
287 
288  intPart = ((xBits.i)>>23); /* get biased exponent */
289  intPart -= 127; /* unbias it */
290 
291  x = (float)(xBits.i & 0x007FFFFF); /* mask off exponent leaving 0x800000*(mantissa - 1) */
292  x *= 1.192092895507812e-07f; /* divide by 0x800000 */
293 
294  accumulator = 1.0f + 0.49959804148061f*x;
295  xPower = x*x;
296  accumulator += -0.12047308243453f*xPower;
297  xPower *= x;
298  accumulator += 0.04585425015501f*xPower;
299  xPower *= x;
300  accumulator += -0.01076564682800f*xPower;
301 
302  if (intPart & 0x00000001)
303  {
304  accumulator *= ROOT2; /* an odd input exponent means an extra sqrt(2) in the output */
305  }
306 
307  xBits.i = intPart >> 1; /* divide exponent by 2, lose LSB */
308  xBits.i += 127; /* rebias exponent */
309  xBits.i <<= 23; /* move biased exponent into exponent bits */
310 
311  return accumulator * xBits.f;
312  }
313  else
314  {
315  return 0.0f;
316  }
317 #endif
318  }
319 
320 inline float CLAM_log(register float x)
321 {
322 #ifndef CLAM_OPTIMIZE
323  return (float) log((double)x);
324 #else
325  return LN2*CLAM_log2(x);
326 #endif
327 }
328 
329 inline float CLAM_log10(register float x)
330 {
331 #ifndef CLAM_OPTIMIZE
332  return (float) log10((double)x);
333 #else
334  return LN2_OVER_LN10*CLAM_log2(x);
335 #endif
336 }
337 
338 inline float CLAM_20log10(register float x)
339 {
340 #ifndef CLAM_OPTIMIZE
341  return (float) 20*log10((double)x);
342 #else
344 #endif
345 }
346 
347 inline float CLAM_exp(register float x)
348 {
349 #ifndef CLAM_OPTIMIZE
350  return (float) exp((double)x);
351 #else
352  return CLAM_exp2(ONE_OVER_LN2*x);
353 #endif
354 }
355 
356 #if defined _MSC_VER && _MSC_VER < 1310 // MSVC++ 6
357 #undef min
358 #undef max
359  namespace std
360  {
361  template < typename T >
362  const T& max( const T& a, const T& b) {
363  return (a>=b)? a : b;
364  }
365  template < typename T >
366  const T& min(const T& a, const T& b) {
367  return (a<=b)? a : b;
368  }
369  } // namespace std
370 #endif // MSVC++ 6
371 
372 #if defined _MSC_VER // MSVC++7
373  namespace std
374  {
375  template <typename T>
376  bool isnan(T data)
377  {
378  return _isnan(data) == 1;
379  }
380  template <typename T>
381  bool isinf(T data)
382  {
383  return _isnan(data) == 1;
384  }
385  }
386 #endif // MSVC++ 7
387 
388 #ifndef __USE_ISOC99
389 #ifndef __APPLE__
390 inline double round(double _X)
391  {return (floor(_X+0.5)); }
392 inline float round(float _X)
393  {return (floorf(_X+0.5f)); }
394 #endif // __APPLE__
395 #endif // __USE_ISOC99
396 
397 
400 inline float log2lin( float x )
401 {
402 
403 // static double magic = 1.0 / (20.0 * log10(exp(1.0)))=0.1151292546497;
404 
405  return CLAM_exp( x * 0.1151292546497f );
406 
407 }
408 
415 {
416  return (((n - 1) & n) == 0);
417 }
418 
425 {
426  --n;
427 
428  n |= n >> 16;
429  n |= n >> 8;
430  n |= n >> 4;
431  n |= n >> 2;
432  n |= n >> 1;
433  ++n;
434 
435  return n;
436 }
437 
438 namespace CLAM
439 {
440 
441 /*Non member function, returns absolute value of class T*/
442 template <class T> inline T Abs(T value)
443 {
444  return ( value < 0 ) ? -value : value;
445 }
446 
447 /* DB */
448 
449 // Default scaling
450 #define CLAM_DB_SCALING 20
451 
452 inline double DB(double linData, int scaling=20)
453 {
454  return (scaling*CLAM_log10(linData));
455 }
456 
457 inline double Lin(double logData, int scaling=20 )
458 {
459  return (CLAM_pow(double(10),(logData/scaling)) );
460 }
461 
466 template<class T> inline
467 T CLAM_max(const T& x, const T& y)
468  {return (x < y ? y : x); }
469 
470 template<class T> inline
471 T CLAM_min(const T& x, const T& y)
472  {return (x > y ? y : x); }
473 }
474 
475 #endif // CLAM_Math.hxx
476