CLAM-Development
1.1
|
00001 /* 00002 * Copyright (c) 2004 MUSIC TECHNOLOGY GROUP (MTG) 00003 * UNIVERSITAT POMPEU FABRA 00004 * 00005 * 00006 * This program is free software; you can redistribute it and/or modify 00007 * it under the terms of the GNU General Public License as published by 00008 * the Free Software Foundation; either version 2 of the License, or 00009 * (at your option) any later version. 00010 * 00011 * This program is distributed in the hope that it will be useful, 00012 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00013 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00014 * GNU General Public License for more details. 00015 * 00016 * You should have received a copy of the GNU General Public License 00017 * along with this program; if not, write to the Free Software 00018 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00019 * 00020 */ 00021 00022 #ifndef __CLAM_MATH__ 00023 #define __CLAM_MATH__ 00024 00025 #include <cmath> 00026 #include "DataTypes.hxx" 00027 #include "FastRounding.hxx" 00028 00029 00030 //The following constants are defined also in OSDefines but only for windows and using the #define preprocessor 00031 //directive. It is much better to use const float declarations 00032 const float PI_ = 3.1415926535897932384626433832795028841972; /* pi */ 00033 const float ONE_OVER_PI = (0.3183098861837906661338147750939f); 00034 const float TWOPI = (6.2831853071795864769252867665590057683943f); /* 2*pi */ 00035 const float ONE_OVER_TWOPI = (0.15915494309189535682609381638f); 00036 const float PI_2 = (1.5707963267948966192313216916397514420986f); /* pi/2 */ 00037 const float TWO_OVER_PI = (0.636619772367581332267629550188f); 00038 const float LN2 = (0.6931471805599453094172321214581765680755f); /* ln(2) */ 00039 const float ONE_OVER_LN2 = (1.44269504088896333066907387547f); 00040 const float LN10 = (2.3025850929940456840179914546843642076011f); /* ln(10) */ 00041 const float ONE_OVER_LN10 = (0.43429448190325177635683940025f); 00042 const float LN2_OVER_LN10 = LN2*ONE_OVER_LN10; 00043 const float TIMES20LN2_OVER_LN10 = 20*LN2_OVER_LN10; 00044 const long LONG_OFFSET = 4096L; 00045 const float FLOAT_OFFSET = 4096.0; 00046 const float HUGE_ = 1.0e8; 00047 const float ROOT2 = (1.4142135623730950488016887242096980785697f); /* sqrt(2) */ 00048 00050 inline float CLAM_sin(register float x) 00051 { 00052 #ifndef CLAM_OPTIMIZE 00053 return (float) sin((double)x); 00054 #else 00055 x *= ONE_OVER_PI; 00056 register float accumulator, xPower, xSquared; 00057 register long evenIntPart = ((long)(0.5f*x + 1024.5) - 1024)<<1; 00058 x -= (float)evenIntPart; 00059 xSquared = x*x; 00060 accumulator = 3.14159265358979f*x; 00061 xPower = xSquared*x; 00062 accumulator += -5.16731953364340f*xPower; 00063 xPower *= xSquared; 00064 accumulator += 2.54620566822659f*xPower; 00065 xPower *= xSquared; 00066 accumulator += -0.586027023087261f*xPower; 00067 xPower *= xSquared; 00068 accumulator += 0.06554823491427f*xPower; 00069 return accumulator; 00070 #endif 00071 } 00072 00073 inline float CLAM_cos(register float x) 00074 { 00075 #ifndef CLAM_OPTIMIZE 00076 return (float) cos((double)x); 00077 #else 00078 x *= ONE_OVER_PI; 00079 register float accumulator, xPower, xSquared; 00080 00081 register long evenIntPart = ((long)(0.5f*x + 1024.5f) - 1024)<<1; 00082 x -= (float)evenIntPart; 00083 00084 xSquared = x*x; 00085 accumulator = 1.57079632679490f*x; /* series for sin(PI/2*x) */ 00086 xPower = xSquared*x; 00087 accumulator += -0.64596406188166f*xPower; 00088 xPower *= xSquared; 00089 accumulator += 0.07969158490912f*xPower; 00090 xPower *= xSquared; 00091 accumulator += -0.00467687997706f*xPower; 00092 xPower *= xSquared; 00093 accumulator += 0.00015303015470f*xPower; 00094 return 1.0f - 2.0f*accumulator*accumulator; /* cos(w) = 1 - 2*(sin(w/2))^2 */ 00095 #endif 00096 } 00097 00098 inline float CLAM_atan(register float x) 00099 { 00100 #ifndef CLAM_OPTIMIZE 00101 return (float) atan((double)x); 00102 #else 00103 register float accumulator, xPower, xSquared, offset; 00104 00105 offset = 0.0f; 00106 00107 if (x < -1.0f) 00108 { 00109 offset = -PI_2; 00110 x = -1.0f/x; 00111 } 00112 else if (x > 1.0f) 00113 { 00114 offset = PI_2; 00115 x = -1.0f/x; 00116 } 00117 xSquared = x*x; 00118 accumulator = 1.0f; 00119 xPower = xSquared; 00120 accumulator += 0.33288950512027f*xPower; 00121 xPower *= xSquared; 00122 accumulator += -0.08467922817644f*xPower; 00123 xPower *= xSquared; 00124 accumulator += 0.03252232640125f*xPower; 00125 xPower *= xSquared; 00126 accumulator += -0.00749305860992f*xPower; 00127 00128 return offset + x/accumulator; 00129 #endif 00130 } 00131 00132 inline float CLAM_atan2(float Imag, float Real) 00133 { 00134 #ifndef CLAM_OPTIMIZE 00135 return (float) atan2((double)Imag, (double)Real); 00136 #else 00137 if(Real==0 && Imag==0) return 0.f; 00138 register float accumulator, xPower, xSquared, offset, x; 00139 00140 if (Imag > 0.0f) 00141 { 00142 if (Imag <= -Real) 00143 { 00144 offset = PI_; 00145 x = Imag/Real; 00146 } 00147 else if (Imag > Real) 00148 { 00149 offset = PI_2; 00150 x = -Real/Imag; 00151 } 00152 else 00153 { 00154 offset = 0.0f; 00155 x = Imag/Real; 00156 } 00157 } 00158 else 00159 { 00160 if (Imag >= Real) 00161 { 00162 offset = -PI_; 00163 x = Imag/Real; 00164 } 00165 else if (Imag < -Real) 00166 { 00167 offset = -PI_2; 00168 x = -Real/Imag; 00169 } 00170 else 00171 { 00172 offset = 0.0f; 00173 x = Imag/Real; 00174 } 00175 } 00176 00177 xSquared = x*x; 00178 accumulator = 1.0f; 00179 xPower = xSquared; 00180 accumulator += 0.33288950512027f*xPower; 00181 xPower *= xSquared; 00182 accumulator += -0.08467922817644f*xPower; 00183 xPower *= xSquared; 00184 accumulator += 0.03252232640125f*xPower; 00185 xPower *= xSquared; 00186 accumulator += -0.00749305860992f*xPower; 00187 00188 return offset + x/accumulator; 00189 #endif 00190 } 00191 00192 inline float CLAM_exp2(register float x) 00193 { 00194 #ifndef CLAM_OPTIMIZE 00195 return (float) exp(LN2*(double)x); 00196 #else 00197 if (x >= -127.0f) 00198 { 00199 register float accumulator, xPower; 00200 register union {float f; long i;} xBits; 00201 00202 xBits.i = (long)(x + FLOAT_OFFSET) - LONG_OFFSET; /* integer part */ 00203 x -= (float)(xBits.i); /* fractional part */ 00204 00205 accumulator = 1.0f + 0.69303212081966f*x; 00206 xPower = x*x; 00207 accumulator += 0.24137976293709f*xPower; 00208 xPower *= x; 00209 accumulator += 0.05203236900844f*xPower; 00210 xPower *= x; 00211 accumulator += 0.01355574723481f*xPower; 00212 00213 xBits.i += 127; /* bias integer part */ 00214 xBits.i <<= 23; /* move biased int part into exponent bits */ 00215 00216 return accumulator * xBits.f; 00217 } 00218 else 00219 { 00220 return 0.0f; 00221 } 00222 #endif 00223 } 00224 00225 inline float CLAM_log2(register float x) 00226 { 00227 #ifndef CLAM_OPTIMIZE 00228 return (float) (ONE_OVER_LN2*log((double)x)); 00229 #else 00230 if (x > 5.877471754e-39f) 00231 { 00232 register float accumulator, xPower; 00233 register long intPart; 00234 00235 register union {float f; long i;} xBits; 00236 00237 xBits.f = x; 00238 00239 intPart = ((xBits.i)>>23); /* get biased exponent */ 00240 intPart -= 127; /* unbias it */ 00241 00242 x = (float)(xBits.i & 0x007FFFFF); /* mask off exponent leaving 0x800000*(mantissa - 1) */ 00243 x *= 1.192092895507812e-07f; /* divide by 0x800000 */ 00244 00245 accumulator = 1.44254494359510f*x; 00246 xPower = x*x; 00247 accumulator += -0.71814525675041f*xPower; 00248 xPower *= x; 00249 accumulator += 0.45754919692582f*xPower; 00250 xPower *= x; 00251 accumulator += -0.27790534462866f*xPower; 00252 xPower *= x; 00253 accumulator += 0.12179791068782f*xPower; 00254 xPower *= x; 00255 accumulator += -0.02584144982967f*xPower; 00256 00257 return accumulator + (float)intPart; 00258 } 00259 else 00260 { 00261 return -HUGE_; 00262 } 00263 #endif 00264 } 00265 00266 inline float CLAM_pow(float x, float y) 00267 { 00268 #ifndef CLAM_OPTIMIZE 00269 return (float) pow((double)x, (double)y); 00270 #else 00271 return CLAM_exp2(y*CLAM_log2(x)); 00272 #endif 00273 } 00274 00275 inline float CLAM_sqrt(register float x) 00276 { 00277 #ifndef CLAM_OPTIMIZE 00278 return (float) sqrt((double)x); 00279 #else 00280 if (x > 5.877471754e-39f) 00281 { 00282 register float accumulator, xPower; 00283 register long intPart; 00284 register union {float f; long i;} xBits; 00285 00286 xBits.f = x; 00287 00288 intPart = ((xBits.i)>>23); /* get biased exponent */ 00289 intPart -= 127; /* unbias it */ 00290 00291 x = (float)(xBits.i & 0x007FFFFF); /* mask off exponent leaving 0x800000*(mantissa - 1) */ 00292 x *= 1.192092895507812e-07f; /* divide by 0x800000 */ 00293 00294 accumulator = 1.0f + 0.49959804148061f*x; 00295 xPower = x*x; 00296 accumulator += -0.12047308243453f*xPower; 00297 xPower *= x; 00298 accumulator += 0.04585425015501f*xPower; 00299 xPower *= x; 00300 accumulator += -0.01076564682800f*xPower; 00301 00302 if (intPart & 0x00000001) 00303 { 00304 accumulator *= ROOT2; /* an odd input exponent means an extra sqrt(2) in the output */ 00305 } 00306 00307 xBits.i = intPart >> 1; /* divide exponent by 2, lose LSB */ 00308 xBits.i += 127; /* rebias exponent */ 00309 xBits.i <<= 23; /* move biased exponent into exponent bits */ 00310 00311 return accumulator * xBits.f; 00312 } 00313 else 00314 { 00315 return 0.0f; 00316 } 00317 #endif 00318 } 00319 00320 inline float CLAM_log(register float x) 00321 { 00322 #ifndef CLAM_OPTIMIZE 00323 return (float) log((double)x); 00324 #else 00325 return LN2*CLAM_log2(x); 00326 #endif 00327 } 00328 00329 inline float CLAM_log10(register float x) 00330 { 00331 #ifndef CLAM_OPTIMIZE 00332 return (float) log10((double)x); 00333 #else 00334 return LN2_OVER_LN10*CLAM_log2(x); 00335 #endif 00336 } 00337 00338 inline float CLAM_20log10(register float x) 00339 { 00340 #ifndef CLAM_OPTIMIZE 00341 return (float) 20*log10((double)x); 00342 #else 00343 return TIMES20LN2_OVER_LN10*CLAM_log2(x); 00344 #endif 00345 } 00346 00347 inline float CLAM_exp(register float x) 00348 { 00349 #ifndef CLAM_OPTIMIZE 00350 return (float) exp((double)x); 00351 #else 00352 return CLAM_exp2(ONE_OVER_LN2*x); 00353 #endif 00354 } 00355 00356 #if defined _MSC_VER && _MSC_VER < 1310 // MSVC++ 6 00357 #undef min 00358 #undef max 00359 namespace std 00360 { 00361 template < typename T > 00362 const T& max( const T& a, const T& b) { 00363 return (a>=b)? a : b; 00364 } 00365 template < typename T > 00366 const T& min(const T& a, const T& b) { 00367 return (a<=b)? a : b; 00368 } 00369 } // namespace std 00370 #endif // MSVC++ 6 00371 00372 #if defined _MSC_VER // MSVC++7 00373 namespace std 00374 { 00375 template <typename T> 00376 bool isnan(T data) 00377 { 00378 return _isnan(data) == 1; 00379 } 00380 template <typename T> 00381 bool isinf(T data) 00382 { 00383 return _isnan(data) == 1; 00384 } 00385 } 00386 #endif // MSVC++ 7 00387 00388 #ifndef __USE_ISOC99 00389 #ifndef __APPLE__ 00390 inline double round(double _X) 00391 {return (floor(_X+0.5)); } 00392 inline float round(float _X) 00393 {return (floorf(_X+0.5f)); } 00394 #endif // __APPLE__ 00395 #endif // __USE_ISOC99 00396 00397 00400 inline float log2lin( float x ) 00401 { 00402 00403 // static double magic = 1.0 / (20.0 * log10(exp(1.0)))=0.1151292546497; 00404 00405 return CLAM_exp( x * 0.1151292546497f ); 00406 00407 } 00408 00414 inline bool isPowerOfTwo( CLAM::TUInt32 n) 00415 { 00416 return (((n - 1) & n) == 0); 00417 } 00418 00424 inline CLAM::TUInt32 nextPowerOfTwo( CLAM::TUInt32 n) 00425 { 00426 --n; 00427 00428 n |= n >> 16; 00429 n |= n >> 8; 00430 n |= n >> 4; 00431 n |= n >> 2; 00432 n |= n >> 1; 00433 ++n; 00434 00435 return n; 00436 } 00437 00438 namespace CLAM 00439 { 00440 00441 /*Non member function, returns absolute value of class T*/ 00442 template <class T> inline T Abs(T value) 00443 { 00444 return ( value < 0 ) ? -value : value; 00445 } 00446 00447 /* DB */ 00448 00449 // Default scaling 00450 #define CLAM_DB_SCALING 20 00451 00452 inline double DB(double linData, int scaling=20) 00453 { 00454 return (scaling*CLAM_log10(linData)); 00455 } 00456 00457 inline double Lin(double logData, int scaling=20 ) 00458 { 00459 return (CLAM_pow(double(10),(logData/scaling)) ); 00460 } 00461 00466 template<class T> inline 00467 T CLAM_max(const T& x, const T& y) 00468 {return (x < y ? y : x); } 00469 00470 template<class T> inline 00471 T CLAM_min(const T& x, const T& y) 00472 {return (x > y ? y : x); } 00473 } 00474 00475 #endif // CLAM_Math.hxx 00476