00001 #ifndef __CMathTools__
00002 #define __CMathTools__
00003
00004 #include "../Basics/CCountedObject.hpp"
00005 #include <algorithm>
00006 using namespace std;
00007
00008 #ifdef WIN32
00009 #define IMAN 0
00010 #else
00011 #ifdef __i386
00012 #define IMAN 0
00013 #else
00014 #define IMAN 1
00015 #endif
00016 #endif
00017
00018 #ifdef WIN32
00019 #define exponent_position 1
00020 #else
00021 #define exponent_position 0
00022 #endif
00023
00024
00025 static float shift23=(1<<23);
00026 static float PowBodge=0.33971f;
00027
00028
00029
00030 namespace Exponent
00031 {
00032 namespace MathTools
00033 {
00052 class CMathTools
00053 {
00054 public:
00055
00056
00057
00058 const static float CMATH_HALF_PI_FLOAT;
00059 const static float CMATH_PI_FLOAT;
00060 const static float CMATH_2PI_FLOAT;
00061 const static float CMATH_4PI_FLOAT;
00063
00064
00065 const static double CMATH_HALF_PI_DOUBLE;
00066 const static double CMATH_PI_DOUBLE;
00067 const static double CMATH_2PI_DOUBLE;
00068 const static double CMATH_4PI_DOUBLE;
00070
00071
00072 const static double CMATH_LN2_DOUBLE;
00073 const static double CMATH_LN2_INV_DOUBLE;
00075
00076
00077 const static float CMATH_LN2_FLOAT;
00078 const static float CMATH_LN2_INV_FLOAT;
00080
00081
00082 const static double CMATH_ONE_THIRD_DOUBLE;
00083 const static double CMATH_TWO_THIRDS_DOUBLE;
00085
00086
00087 const static float CMATH_ONE_THIRD_FLOAT;
00088 const static float CMATH_TWO_THIRDS_FLOAT;
00090
00091
00098 static FORCEINLINE double minimum(const double a, const double b)
00099 {
00100 return min(a, b);
00101 }
00102
00109 static FORCEINLINE float minimum(const float a, const float b)
00110 {
00111 return min(a, b);
00112 }
00113
00120 static FORCEINLINE long minimum(const long a, const long b)
00121 {
00122 return min(a, b);
00123 }
00124
00131 static FORCEINLINE double fastMinimum(const double a, const double b)
00132 {
00133 return 0.5 * (a + b - fabs(a - b));
00134 }
00135
00142 static FORCEINLINE float fastMinimum(const float a, const float b)
00143 {
00144 return 0.5f * (a + b - fabsf(a - b));
00145 }
00146
00147
00148
00155 static FORCEINLINE double maximum(const double a, const double b)
00156 {
00157 return max(a, b);
00158 }
00159
00166 static FORCEINLINE float maximum(const float a, const float b)
00167 {
00168 return max(a, b);
00169 }
00170
00177 static FORCEINLINE long maximum(const long a, const long b)
00178 {
00179 return max(a, b);
00180 }
00181
00188 static FORCEINLINE double fastMaximum(const double a, const double b)
00189 {
00190 return 0.5 * (a + b + fabs(a - b));
00191 }
00192
00199 static FORCEINLINE float fastMaximum(const float a, const float b)
00200 {
00201 return 0.5f * (a + b + fabsf(a - b));
00202 }
00203
00204
00205
00211 static FORCEINLINE double fastModulus(const double x)
00212 {
00213 return x - floor(x);
00214 }
00215
00221 static FORCEINLINE float fastModulus(const float x)
00222 {
00223 return x - floorf(x);
00224 }
00225
00226
00227
00234 static FORCEINLINE double fastAtan(const double x)
00235 {
00236 return (x / (1.0 + 0.28 * (x * x)));
00237 }
00238
00245 static FORCEINLINE float fastAtan(const float x)
00246 {
00247 return (x / (1.f + 0.28f * (x * x)));
00248 }
00249
00250
00251
00259 static FORCEINLINE double clamp(const double min, const double max, const double x)
00260 {
00261 return 0.5 * (fabs(x - min) + (min + max) - fabs(x - max));
00262 }
00263
00269 static FORCEINLINE double clamp(const double x)
00270 {
00271 return 0.5 * (fabs(x) + 1.0 - fabs(x - 1.0));
00272 }
00273
00281 static FORCEINLINE float clamp(const float min, const float max, const float x)
00282 {
00283 return 0.5f * (fabsf(x - min) + (min + max) - fabsf(x - max));
00284 }
00285
00286
00287
00293 static FORCEINLINE double cube(const double x)
00294 {
00295 return x * x * x;
00296 }
00297
00303 static FORCEINLINE float cube(const float x)
00304 {
00305 return x * x * x;
00306 }
00307
00308
00309
00315 static FORCEINLINE double square(const double x)
00316 {
00317 return x * x;
00318 }
00319
00325 static FORCEINLINE float square(const float x)
00326 {
00327 return x * x;
00328 }
00329
00330
00331
00338 static FORCEINLINE long fastDoubleToLong(double x)
00339 {
00340 #ifdef WIN32
00341 x = x + 103079215104;
00342 #else
00343 #ifndef __i386
00344 x = x + 103079215104.0;
00345 #else
00346 return (long)x;
00347 #endif
00348 #endif
00349 return ((long*)&x)[IMAN] >> 16;
00350 }
00351
00359 static FORCEINLINE long fastDoubleToLongFloor(const double x)
00360 {
00361 #ifdef WIN32
00362 long returnValue;
00363 __asm
00364 {
00365 fld x
00366 fistp returnValue
00367 }
00368 return returnValue;
00369 #else
00370 return (long)floor(x);
00371 #endif
00372 }
00373
00374
00375
00382 static FORCEINLINE double iZero(const double y)
00383 {
00384 double s = 1.0;
00385 double ds = 1.0;
00386 double d = 0.0;
00387 do
00388 {
00389 d = d + 2.0;
00390 ds = ds * (y * y) / (d * d);
00391 s = s + ds;
00392 }while(ds > 1e-7 * s);
00393 return s;
00394 }
00395
00402 static FORCEINLINE double zeroethOrderBessel(const double x)
00403 {
00404 double bessel = 1.0;
00405 double term = 0.0;
00406 long i = 1;
00407
00408 do
00409 {
00410
00411 term = pow(0.5 * x, (double)i) / factorial((double)i);
00412
00413
00414 bessel += (term * term);
00415
00416
00417 i++;
00418 }while (term > 0.000001 * bessel);
00419
00420
00421 return bessel;
00422 }
00423
00429 static FORCEINLINE double factorial(const double x)
00430 {
00431
00432 double outputValue = 1.0;
00433 double value = x;
00434
00435
00436 while (value > 1)
00437 {
00438 outputValue *= value--;
00439 }
00440
00441
00442 return outputValue;
00443 }
00444
00445
00446
00449
00450
00458 static FORCEINLINE int fastSign(const float x)
00459 {
00460 #ifdef WIN32
00461 return 1 + (((*(int *) &x) >> 31) << 1);
00462 #else
00463 if (x > 0)
00464 {
00465 return 1;
00466 }
00467 else if (x < 0)
00468 {
00469 return -1;
00470 }
00471 return 0;
00472 #endif
00473 }
00474
00482 static FORCEINLINE int fastSign(const double x)
00483 {
00484 #ifdef WIN32
00485 return fastSign((float)x);
00486 #else
00487 if (x > 0)
00488 {
00489 return 1;
00490 }
00491 else if (x < 0)
00492 {
00493 return -1;
00494 }
00495 return 0;
00496 #endif
00497 }
00498
00499
00500
00506 static FORCEINLINE double fastLog(const double x)
00507 {
00508 #ifdef WIN32
00509 return (double)fastLog((float)x);
00510 #else
00511 return logf(x);
00512 #endif
00513 }
00514
00520 static FORCEINLINE float fastLog(const float x)
00521 {
00522 #ifdef WIN32
00523 return (fastLog2(x) * 0.69314718f);
00524 #else
00525 return log(x);
00526 #endif
00527 }
00528
00529
00530
00536 static FORCEINLINE double fastAbsolute(const double x)
00537 {
00538 #ifdef WIN32
00539 return (double)fastAbsolute((float)x);
00540 #else
00541 return fabs(x);
00542 #endif
00543 }
00544
00550 static FORCEINLINE float fastAbsolute(const float x)
00551 {
00552 #ifdef WIN32
00553 int i = ((*(int*)&x)&0x7fffffff);
00554 return (*(float*)&i);
00555 #else
00556 return fabsf(x);
00557 #endif
00558 }
00559
00560
00561
00562
00569 static FORCEINLINE double fastSquareRoot(const double x)
00570 {
00571 #ifdef WIN32
00572 double y, z, temp;
00573 unsigned long *tempPtr = ((unsigned long *)&temp)+1;
00574
00575 temp = x;
00576 *tempPtr = (0xbfcd4600 - *tempPtr) >> 1;
00577 y = temp;
00578 z = x * 0.5;
00579 y *= 1.5 - y * y * z;
00580 y *= 1.5 - y * y * z;
00581 y *= 1.5 - y * y * z;
00582 y *= 1.5 - y * y * z;
00583 return y * x;
00584 #else
00585 return sqrt(x);
00586 #endif
00587 }
00588
00595 static FORCEINLINE float fastSquareRoot(const float x)
00596 {
00597 #ifdef WIN32
00598 return (float)fastSquareRoot((double)x);
00599 #else
00600 return sqrtf(x);
00601 #endif
00602 }
00603
00604
00611 static FORCEINLINE float fastSquareRootVersion2(const float x)
00612 {
00613 #ifdef WIN32
00614 float result;
00615 _asm
00616 {
00617 mov eax, x
00618 sub eax, 0x3f800000
00619 sar eax, 1
00620 add eax, 0x3f800000
00621 mov result, eax
00622 }
00623 return result;
00624 #else
00625 return sqrtf(x);
00626 #endif
00627 }
00628
00629
00630
00637 static FORCEINLINE double fastExponential(const double x)
00638 {
00639 return myPow2(x * CMATH_LN2_INV_DOUBLE);
00640 }
00641
00642
00643 static FORCEINLINE float myPow2(float i)
00644 {
00645
00646 float x;
00647 float y=i-floorf(i);
00648 y=(y-y*y)*PowBodge;
00649
00650 x=i+127-y;
00651 x*= shift23;
00652 *(int*)&x=(int)x;
00653 return x;
00654 }
00655
00656 static FORCEINLINE double myPow2(double x)
00657 {
00658
00659 long *p = (long*)(&x) + exponent_position;
00660
00661 const long tx = fastDoubleToLongFloor(x);
00662
00663 x -= static_cast<double>(tx);
00664
00665
00666 x = 1.0 + x * (0.6 + x * 0.3);
00667
00668 *p += (tx<<20);
00669 return x;
00670 }
00671
00677 static FORCEINLINE double fastExp2(const double x)
00678 {
00679 #ifdef WIN32
00680 int e;
00681 double ret;
00682
00683 if (x >= 0)
00684 {
00685 e = int (x);
00686 ret = x - (e - 1);
00687 ((*(1 + (int *) &ret)) &= ~(2047 << 20)) += (e + 1023) << 20;
00688 }
00689 else
00690 {
00691 e = int (x + 1023);
00692 ret = x - (e - 1024);
00693 ((*(1 + (int *) &ret)) &= ~(2047 << 20)) += e << 20;
00694 }
00695 return (ret);
00696 #else
00697 return exp2(x);
00698 #endif
00699 }
00700
00706 static FORCEINLINE float fastExp2(const float x)
00707 {
00708 #ifdef WIN32
00709 return (float)fastExp2((double)x);
00710 #else
00711 return exp2f(x);
00712 #endif
00713 }
00714
00723 static FORCEINLINE double fastExp(const double x)
00724 {
00725
00726
00727
00728
00729
00730
00731
00732
00733
00734
00735
00736
00737
00738
00739
00740
00741
00742
00743
00744
00745
00746
00747 union
00748 {
00749 double d;
00750 #ifdef WIN32
00751 struct { int j, i; } n;
00752 #else
00753 struct { int i, j; } n;
00754 #endif
00755 }eco;
00756
00757 eco.n.i = (int)(1512775.3951951856938358403823062 * x) + (1072632447);
00758 eco.n.j = 0;
00759
00760 return eco.d;
00761 }
00762
00763
00764
00770 static FORCEINLINE double fastPower4(double x)
00771 {
00772 #ifdef WIN32
00773 return (double)fastPower4((float)x);
00774 #else
00775 return x * x * x * x;
00776 #endif
00777 }
00778
00784 static FORCEINLINE float fastPower4(float x)
00785 {
00786 #ifdef WIN32
00787 long *longPointer = (long *)(&x);
00788 long myLong = *longPointer;
00789 myLong -= 0x3F800000l;
00790 myLong <<= 2;
00791 myLong += 0x3F800000l;
00792 *longPointer = myLong;
00793 return x;
00794 #else
00795 return x * x * x * x;
00796 #endif
00797 }
00798
00806 static FORCEINLINE double fastPower(const double b, const double e)
00807 {
00808 return exp(e * log(b));
00809 }
00810
00811
00812
00818 static FORCEINLINE double fastSin(const double x)
00819 {
00820 const double xSquared = x * x;
00821 double returnValue = 0.00761;
00822 returnValue *= xSquared;
00823 returnValue -= 0.16605;
00824 returnValue *= xSquared;
00825 returnValue += 1.0;
00826 returnValue *= x;
00827 return returnValue;
00828 }
00829
00830
00831
00837 static FORCEINLINE double tanhApproximation(const double x)
00838 {
00839 double a = fabs(x);
00840 a*=(6+a*(3+a));
00841 return ((x<0)?-1:1)*a/(a+12);
00842
00843
00844
00845 }
00846
00852 static FORCEINLINE double fastTanh(const double x)
00853 {
00854
00855 return 1.0 - 2.0 / (fastExp(2.0 * x) + 1.0);
00856 }
00857
00863 static FORCEINLINE float fastTanh(const float x)
00864 {
00865 return 1.f - 2.f / (expf(2.f * x) + 1.f);
00866 }
00867
00868 #ifdef WIN32
00869
00875 static FORCEINLINE double fastTanh2(const double x)
00876 {
00877 double absolute = fabs(x);
00878 absolute *= (6 + absolute * (3 + absolute));
00879 return fastSign(x) * absolute / (absolute + 12);
00880 }
00881
00887 static FORCEINLINE float fastTanh2(const float x)
00888 {
00889 float absolute = fabsf(x);
00890 absolute *= (6 + absolute * (3 + absolute));
00891 return fastSign(x) * absolute / (absolute + 12);
00892 }
00893 #endif
00894
00895
00896
00897 protected:
00898
00899
00900 #ifdef WIN32
00901
00906 static FORCEINLINE float fastLog2(float x)
00907 {
00908 int * const exp_ptr = reinterpret_cast <int *> (&x);
00909 int val = *exp_ptr;
00910 const int log_2 = ((val >> 23) & 255) - 128;
00911 val &= ~(255 << 23);
00912 val += 127 << 23;
00913 *exp_ptr = val;
00914
00915 return (x + log_2);
00916 }
00917 #endif
00918
00921 };
00922 }
00923 }
00924 #endif // End of CMathTools.hpp