neo/idlib/math/Math.h

   1 /*
   2 ===========================================================================
   3
   4 Doom 3 GPL Source Code
   5 Copyright (C) 1999-2011 id Software LLC, a ZeniMax Media company.
   6
   7 This file is part of the Doom 3 GPL Source Code (?Doom 3 Source Code?).
   8
   9 Doom 3 Source Code is free software: you can redistribute it and/or modify
  10 it under the terms of the GNU General Public License as published by
  11 the Free Software Foundation, either version 3 of the License, or
  12 (at your option) any later version.
  13
  14 Doom 3 Source Code is distributed in the hope that it will be useful,
  15 but WITHOUT ANY WARRANTY; without even the implied warranty of
  16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  17 GNU General Public License for more details.
  18
  19 You should have received a copy of the GNU General Public License
  20 along with Doom 3 Source Code.  If not, see <http://www.gnu.org/licenses/>.
  21
  22 In addition, the Doom 3 Source Code is also subject to certain additional terms. You should have received a copy of these additional terms immediately following the terms and conditions of the GNU General Public License which accompanied the Doom 3 Source Code.  If not, please request a copy in writing from id Software at the address below.
  23
  24 If you have questions concerning this license or the applicable additional terms, you may contact in writing id Software LLC, c/o ZeniMax Media Inc., Suite 120, Rockville, Maryland 20850 USA.
  25
  26 ===========================================================================
  27 */
  28
  29 #ifndef __MATH_MATH_H__
  30 #define __MATH_MATH_H__
  31
  32 #if defined(MACOS_X) && defined(__powerpc__)
  33 // for square root estimate instruction
  34 #include <ppc_intrinsics.h>
  35 // for FLT_MIN
  36 #include <float.h>
  37 #endif
  38 /*
  39 ===============================================================================
  40
  41   Math
  42
  43 ===============================================================================
  44 */
  45
  46 #ifdef INFINITY
  47 #undef INFINITY
  48 #endif
  49
  50 #ifdef FLT_EPSILON
  51 #undef FLT_EPSILON
  52 #endif
  53
  54 #define DEG2RAD(a)                              ( (a) * idMath::M_DEG2RAD )
  55 #define RAD2DEG(a)                              ( (a) * idMath::M_RAD2DEG )
  56
  57 #define SEC2MS(t)                               ( idMath::FtoiFast( (t) * idMath::M_SEC2MS ) )
  58 #define MS2SEC(t)                               ( (t) * idMath::M_MS2SEC )
  59
  60 #define ANGLE2SHORT(x)                  ( idMath::FtoiFast( (x) * 65536.0f / 360.0f ) & 65535 )
  61 #define SHORT2ANGLE(x)                  ( (x) * ( 360.0f / 65536.0f ) )
  62
  63 #define ANGLE2BYTE(x)                   ( idMath::FtoiFast( (x) * 256.0f / 360.0f ) & 255 )
  64 #define BYTE2ANGLE(x)                   ( (x) * ( 360.0f / 256.0f ) )
  65
  66 #define FLOATSIGNBITSET(f)              ((*(const unsigned long *)&(f)) >> 31)
  67 #define FLOATSIGNBITNOTSET(f)   ((~(*(const unsigned long *)&(f))) >> 31)
  68 #define FLOATNOTZERO(f)                 ((*(const unsigned long *)&(f)) & ~(1<<31) )
  69 #define INTSIGNBITSET(i)                (((const unsigned long)(i)) >> 31)
  70 #define INTSIGNBITNOTSET(i)             ((~((const unsigned long)(i))) >> 31)
  71
  72 #define FLOAT_IS_NAN(x)                 (((*(const unsigned long *)&x) & 0x7f800000) == 0x7f800000)
  73 #define FLOAT_IS_INF(x)                 (((*(const unsigned long *)&x) & 0x7fffffff) == 0x7f800000)
  74 #define FLOAT_IS_IND(x)                 ((*(const unsigned long *)&x) == 0xffc00000)
  75 #define FLOAT_IS_DENORMAL(x)    (((*(const unsigned long *)&x) & 0x7f800000) == 0x00000000 && \
  76                                                                  ((*(const unsigned long *)&x) & 0x007fffff) != 0x00000000 )
  77
  78 #define IEEE_FLT_MANTISSA_BITS  23
  79 #define IEEE_FLT_EXPONENT_BITS  8
  80 #define IEEE_FLT_EXPONENT_BIAS  127
  81 #define IEEE_FLT_SIGN_BIT               31
  82
  83 #define IEEE_DBL_MANTISSA_BITS  52
  84 #define IEEE_DBL_EXPONENT_BITS  11
  85 #define IEEE_DBL_EXPONENT_BIAS  1023
  86 #define IEEE_DBL_SIGN_BIT               63
  87
  88 #define IEEE_DBLE_MANTISSA_BITS 63
  89 #define IEEE_DBLE_EXPONENT_BITS 15
  90 #define IEEE_DBLE_EXPONENT_BIAS 0
  91 #define IEEE_DBLE_SIGN_BIT              79
  92
  93 template<class T> ID_INLINE int MaxIndex( T x, T y ) { return  ( x > y ) ? 0 : 1; }
  94 template<class T> ID_INLINE int MinIndex( T x, T y ) { return ( x < y ) ? 0 : 1; }
  95
  96 template<class T> ID_INLINE T   Max3( T x, T y, T z ) { return ( x > y ) ? ( ( x > z ) ? x : z ) : ( ( y > z ) ? y : z ); }
  97 template<class T> ID_INLINE T   Min3( T x, T y, T z ) { return ( x < y ) ? ( ( x < z ) ? x : z ) : ( ( y < z ) ? y : z ); }
  98 template<class T> ID_INLINE int Max3Index( T x, T y, T z ) { return ( x > y ) ? ( ( x > z ) ? 0 : 2 ) : ( ( y > z ) ? 1 : 2 ); }
  99 template<class T> ID_INLINE int Min3Index( T x, T y, T z ) { return ( x < y ) ? ( ( x < z ) ? 0 : 2 ) : ( ( y < z ) ? 1 : 2 ); }
 100
 101 template<class T> ID_INLINE T   Sign( T f ) { return ( f > 0 ) ? 1 : ( ( f < 0 ) ? -1 : 0 ); }
 102 template<class T> ID_INLINE T   Square( T x ) { return x * x; }
 103 template<class T> ID_INLINE T   Cube( T x ) { return x * x * x; }
 104
 105
 106 class idMath {
 107 public:
 108
 109         static void                                     Init( void );
 110
 111         static float                            RSqrt( float x );                       // reciprocal square root, returns huge number when x == 0.0
 112
 113         static float                            InvSqrt( float x );                     // inverse square root with 32 bits precision, returns huge number when x == 0.0
 114         static float                            InvSqrt16( float x );           // inverse square root with 16 bits precision, returns huge number when x == 0.0
 115         static double                           InvSqrt64( float x );           // inverse square root with 64 bits precision, returns huge number when x == 0.0
 116
 117         static float                            Sqrt( float x );                        // square root with 32 bits precision
 118         static float                            Sqrt16( float x );                      // square root with 16 bits precision
 119         static double                           Sqrt64( float x );                      // square root with 64 bits precision
 120
 121         static float                            Sin( float a );                         // sine with 32 bits precision
 122         static float                            Sin16( float a );                       // sine with 16 bits precision, maximum absolute error is 2.3082e-09
 123         static double                           Sin64( float a );                       // sine with 64 bits precision
 124
 125         static float                            Cos( float a );                         // cosine with 32 bits precision
 126         static float                            Cos16( float a );                       // cosine with 16 bits precision, maximum absolute error is 2.3082e-09
 127         static double                           Cos64( float a );                       // cosine with 64 bits precision
 128
 129         static void                                     SinCos( float a, float &s, float &c );          // sine and cosine with 32 bits precision
 130         static void                                     SinCos16( float a, float &s, float &c );        // sine and cosine with 16 bits precision
 131         static void                                     SinCos64( float a, double &s, double &c );      // sine and cosine with 64 bits precision
 132
 133         static float                            Tan( float a );                         // tangent with 32 bits precision
 134         static float                            Tan16( float a );                       // tangent with 16 bits precision, maximum absolute error is 1.8897e-08
 135         static double                           Tan64( float a );                       // tangent with 64 bits precision
 136
 137         static float                            ASin( float a );                        // arc sine with 32 bits precision, input is clamped to [-1, 1] to avoid a silent NaN
 138         static float                            ASin16( float a );                      // arc sine with 16 bits precision, maximum absolute error is 6.7626e-05
 139         static double                           ASin64( float a );                      // arc sine with 64 bits precision
 140
 141         static float                            ACos( float a );                        // arc cosine with 32 bits precision, input is clamped to [-1, 1] to avoid a silent NaN
 142         static float                            ACos16( float a );                      // arc cosine with 16 bits precision, maximum absolute error is 6.7626e-05
 143         static double                           ACos64( float a );                      // arc cosine with 64 bits precision
 144
 145         static float                            ATan( float a );                        // arc tangent with 32 bits precision
 146         static float                            ATan16( float a );                      // arc tangent with 16 bits precision, maximum absolute error is 1.3593e-08
 147         static double                           ATan64( float a );                      // arc tangent with 64 bits precision
 148
 149         static float                            ATan( float y, float x );       // arc tangent with 32 bits precision
 150         static float                            ATan16( float y, float x );     // arc tangent with 16 bits precision, maximum absolute error is 1.3593e-08
 151         static double                           ATan64( float y, float x );     // arc tangent with 64 bits precision
 152
 153         static float                            Pow( float x, float y );        // x raised to the power y with 32 bits precision
 154         static float                            Pow16( float x, float y );      // x raised to the power y with 16 bits precision
 155         static double                           Pow64( float x, float y );      // x raised to the power y with 64 bits precision
 156
 157         static float                            Exp( float f );                         // e raised to the power f with 32 bits precision
 158         static float                            Exp16( float f );                       // e raised to the power f with 16 bits precision
 159         static double                           Exp64( float f );                       // e raised to the power f with 64 bits precision
 160
 161         static float                            Log( float f );                         // natural logarithm with 32 bits precision
 162         static float                            Log16( float f );                       // natural logarithm with 16 bits precision
 163         static double                           Log64( float f );                       // natural logarithm with 64 bits precision
 164
 165         static int                                      IPow( int x, int y );           // integral x raised to the power y
 166         static int                                      ILog2( float f );                       // integral base-2 logarithm of the floating point value
 167         static int                                      ILog2( int i );                         // integral base-2 logarithm of the integer value
 168
 169         static int                                      BitsForFloat( float f );        // minumum number of bits required to represent ceil( f )
 170         static int                                      BitsForInteger( int i );        // minumum number of bits required to represent i
 171         static int                                      MaskForFloatSign( float f );// returns 0x00000000 if x >= 0.0f and returns 0xFFFFFFFF if x <= -0.0f
 172         static int                                      MaskForIntegerSign( int i );// returns 0x00000000 if x >= 0 and returns 0xFFFFFFFF if x < 0
 173         static int                                      FloorPowerOfTwo( int x );       // round x down to the nearest power of 2
 174         static int                                      CeilPowerOfTwo( int x );        // round x up to the nearest power of 2
 175         static bool                                     IsPowerOfTwo( int x );          // returns true if x is a power of 2
 176         static int                                      BitCount( int x );                      // returns the number of 1 bits in x
 177         static int                                      BitReverse( int x );            // returns the bit reverse of x
 178
 179         static int                                      Abs( int x );                           // returns the absolute value of the integer value (for reference only)
 180         static float                            Fabs( float f );                        // returns the absolute value of the floating point value
 181         static float                            Floor( float f );                       // returns the largest integer that is less than or equal to the given value
 182         static float                            Ceil( float f );                        // returns the smallest integer that is greater than or equal to the given value
 183         static float                            Rint( float f );                        // returns the nearest integer
 184         static int                                      Ftoi( float f );                        // float to int conversion
 185         static int                                      FtoiFast( float f );            // fast float to int conversion but uses current FPU round mode (default round nearest)
 186         static unsigned long            Ftol( float f );                        // float to long conversion
 187         static unsigned long            FtolFast( float );                      // fast float to long conversion but uses current FPU round mode (default round nearest)
 188
 189         static signed char                      ClampChar( int i );
 190         static signed short                     ClampShort( int i );
 191         static int                                      ClampInt( int min, int max, int value );
 192         static float                            ClampFloat( float min, float max, float value );
 193
 194         static float                            AngleNormalize360( float angle );
 195         static float                            AngleNormalize180( float angle );
 196         static float                            AngleDelta( float angle1, float angle2 );
 197
 198         static int                                      FloatToBits( float f, int exponentBits, int mantissaBits );
 199         static float                            BitsToFloat( int i, int exponentBits, int mantissaBits );
 200
 201         static int                                      FloatHash( const float *array, const int numFloats );
 202
 203         static const float                      PI;                                                     // pi
 204         static const float                      TWO_PI;                                         // pi * 2
 205         static const float                      HALF_PI;                                        // pi / 2
 206         static const float                      ONEFOURTH_PI;                           // pi / 4
 207         static const float                      E;                                                      // e
 208         static const float                      SQRT_TWO;                                       // sqrt( 2 )
 209         static const float                      SQRT_THREE;                                     // sqrt( 3 )
 210         static const float                      SQRT_1OVER2;                            // sqrt( 1 / 2 )
 211         static const float                      SQRT_1OVER3;                            // sqrt( 1 / 3 )
 212         static const float                      M_DEG2RAD;                                      // degrees to radians multiplier
 213         static const float                      M_RAD2DEG;                                      // radians to degrees multiplier
 214         static const float                      M_SEC2MS;                                       // seconds to milliseconds multiplier
 215         static const float                      M_MS2SEC;                                       // milliseconds to seconds multiplier
 216         static const float                      INFINITY;                                       // huge number which should be larger than any valid number used
 217         static const float                      FLT_EPSILON;                            // smallest positive number such that 1.0+FLT_EPSILON != 1.0
 218
 219 private:
 220         enum {
 221                 LOOKUP_BITS                             = 8,
 222                 EXP_POS                                 = 23,
 223                 EXP_BIAS                                = 127,
 224                 LOOKUP_POS                              = (EXP_POS-LOOKUP_BITS),
 225                 SEED_POS                                = (EXP_POS-8),
 226                 SQRT_TABLE_SIZE                 = (2<<LOOKUP_BITS),
 227                 LOOKUP_MASK                             = (SQRT_TABLE_SIZE-1)
 228         };
 229
 230         union _flint {
 231                 dword                                   i;
 232                 float                                   f;
 233         };
 234
 235         static dword                            iSqrt[SQRT_TABLE_SIZE];
 236         static bool                                     initialized;
 237 };
 238
 239 ID_INLINE float idMath::RSqrt( float x ) {
 240
 241         long i;
 242         float y, r;
 243
 244         y = x * 0.5f;
 245         i = *reinterpret_cast<long *>( &x );
 246         i = 0x5f3759df - ( i >> 1 );
 247         r = *reinterpret_cast<float *>( &i );
 248         r = r * ( 1.5f - r * r * y );
 249         return r;
 250 }
 251
 252 ID_INLINE float idMath::InvSqrt16( float x ) {
 253
 254         dword a = ((union _flint*)(&x))->i;
 255         union _flint seed;
 256
 257         assert( initialized );
 258
 259         double y = x * 0.5f;
 260         seed.i = (( ( (3*EXP_BIAS-1) - ( (a >> EXP_POS) & 0xFF) ) >> 1)<<EXP_POS) | iSqrt[(a >> (EXP_POS-LOOKUP_BITS)) & LOOKUP_MASK];
 261         double r = seed.f;
 262         r = r * ( 1.5f - r * r * y );
 263         return (float) r;
 264 }
 265
 266 ID_INLINE float idMath::InvSqrt( float x ) {
 267
 268         dword a = ((union _flint*)(&x))->i;
 269         union _flint seed;
 270
 271         assert( initialized );
 272
 273         double y = x * 0.5f;
 274         seed.i = (( ( (3*EXP_BIAS-1) - ( (a >> EXP_POS) & 0xFF) ) >> 1)<<EXP_POS) | iSqrt[(a >> (EXP_POS-LOOKUP_BITS)) & LOOKUP_MASK];
 275         double r = seed.f;
 276         r = r * ( 1.5f - r * r * y );
 277         r = r * ( 1.5f - r * r * y );
 278         return (float) r;
 279 }
 280
 281 ID_INLINE double idMath::InvSqrt64( float x ) {
 282         dword a = ((union _flint*)(&x))->i;
 283         union _flint seed;
 284
 285         assert( initialized );
 286
 287         double y = x * 0.5f;
 288         seed.i = (( ( (3*EXP_BIAS-1) - ( (a >> EXP_POS) & 0xFF) ) >> 1)<<EXP_POS) | iSqrt[(a >> (EXP_POS-LOOKUP_BITS)) & LOOKUP_MASK];
 289         double r = seed.f;
 290         r = r * ( 1.5f - r * r * y );
 291         r = r * ( 1.5f - r * r * y );
 292         r = r * ( 1.5f - r * r * y );
 293         return r;
 294 }
 295
 296 ID_INLINE float idMath::Sqrt16( float x ) {
 297         return x * InvSqrt16( x );
 298 }
 299
 300 ID_INLINE float idMath::Sqrt( float x ) {
 301         return x * InvSqrt( x );
 302 }
 303
 304 ID_INLINE double idMath::Sqrt64( float x ) {
 305         return x * InvSqrt64( x );
 306 }
 307
 308 ID_INLINE float idMath::Sin( float a ) {
 309         return sinf( a );
 310 }
 311
 312 ID_INLINE float idMath::Sin16( float a ) {
 313         float s;
 314
 315         if ( ( a < 0.0f ) || ( a >= TWO_PI ) ) {
 316                 a -= floorf( a / TWO_PI ) * TWO_PI;
 317         }
 318 #if 1
 319         if ( a < PI ) {
 320                 if ( a > HALF_PI ) {
 321                         a = PI - a;
 322                 }
 323         } else {
 324                 if ( a > PI + HALF_PI ) {
 325                         a = a - TWO_PI;
 326                 } else {
 327                         a = PI - a;
 328                 }
 329         }
 330 #else
 331         a = PI - a;
 332         if ( fabs( a ) >= HALF_PI ) {
 333                 a = ( ( a < 0.0f ) ? -PI : PI ) - a;
 334         }
 335 #endif
 336         s = a * a;
 337         return a * ( ( ( ( ( -2.39e-08f * s + 2.7526e-06f ) * s - 1.98409e-04f ) * s + 8.3333315e-03f ) * s - 1.666666664e-01f ) * s + 1.0f );
 338 }
 339
 340 ID_INLINE double idMath::Sin64( float a ) {
 341         return sin( a );
 342 }
 343
 344 ID_INLINE float idMath::Cos( float a ) {
 345         return cosf( a );
 346 }
 347
 348 ID_INLINE float idMath::Cos16( float a ) {
 349         float s, d;
 350
 351         if ( ( a < 0.0f ) || ( a >= TWO_PI ) ) {
 352                 a -= floorf( a / TWO_PI ) * TWO_PI;
 353         }
 354 #if 1
 355         if ( a < PI ) {
 356                 if ( a > HALF_PI ) {
 357                         a = PI - a;
 358                         d = -1.0f;
 359                 } else {
 360                         d = 1.0f;
 361                 }
 362         } else {
 363                 if ( a > PI + HALF_PI ) {
 364                         a = a - TWO_PI;
 365                         d = 1.0f;
 366                 } else {
 367                         a = PI - a;
 368                         d = -1.0f;
 369                 }
 370         }
 371 #else
 372         a = PI - a;
 373         if ( fabs( a ) >= HALF_PI ) {
 374                 a = ( ( a < 0.0f ) ? -PI : PI ) - a;
 375                 d = 1.0f;
 376         } else {
 377                 d = -1.0f;
 378         }
 379 #endif
 380         s = a * a;
 381         return d * ( ( ( ( ( -2.605e-07f * s + 2.47609e-05f ) * s - 1.3888397e-03f ) * s + 4.16666418e-02f ) * s - 4.999999963e-01f ) * s + 1.0f );
 382 }
 383
 384 ID_INLINE double idMath::Cos64( float a ) {
 385         return cos( a );
 386 }
 387
 388 ID_INLINE void idMath::SinCos( float a, float &s, float &c ) {
 389 #ifdef _WIN32
 390         _asm {
 391                 fld             a
 392                 fsincos
 393                 mov             ecx, c
 394                 mov             edx, s
 395                 fstp    dword ptr [ecx]
 396                 fstp    dword ptr [edx]
 397         }
 398 #else
 399         s = sinf( a );
 400         c = cosf( a );
 401 #endif
 402 }
 403
 404 ID_INLINE void idMath::SinCos16( float a, float &s, float &c ) {
 405         float t, d;
 406
 407         if ( ( a < 0.0f ) || ( a >= idMath::TWO_PI ) ) {
 408                 a -= floorf( a / idMath::TWO_PI ) * idMath::TWO_PI;
 409         }
 410 #if 1
 411         if ( a < PI ) {
 412                 if ( a > HALF_PI ) {
 413                         a = PI - a;
 414                         d = -1.0f;
 415                 } else {
 416                         d = 1.0f;
 417                 }
 418         } else {
 419                 if ( a > PI + HALF_PI ) {
 420                         a = a - TWO_PI;
 421                         d = 1.0f;
 422                 } else {
 423                         a = PI - a;
 424                         d = -1.0f;
 425                 }
 426         }
 427 #else
 428         a = PI - a;
 429         if ( fabs( a ) >= HALF_PI ) {
 430                 a = ( ( a < 0.0f ) ? -PI : PI ) - a;
 431                 d = 1.0f;
 432         } else {
 433                 d = -1.0f;
 434         }
 435 #endif
 436         t = a * a;
 437         s = a * ( ( ( ( ( -2.39e-08f * t + 2.7526e-06f ) * t - 1.98409e-04f ) * t + 8.3333315e-03f ) * t - 1.666666664e-01f ) * t + 1.0f );
 438         c = d * ( ( ( ( ( -2.605e-07f * t + 2.47609e-05f ) * t - 1.3888397e-03f ) * t + 4.16666418e-02f ) * t - 4.999999963e-01f ) * t + 1.0f );
 439 }
 440
 441 ID_INLINE void idMath::SinCos64( float a, double &s, double &c ) {
 442 #ifdef _WIN32
 443         _asm {
 444                 fld             a
 445                 fsincos
 446                 mov             ecx, c
 447                 mov             edx, s
 448                 fstp    qword ptr [ecx]
 449                 fstp    qword ptr [edx]
 450         }
 451 #else
 452         s = sin( a );
 453         c = cos( a );
 454 #endif
 455 }
 456
 457 ID_INLINE float idMath::Tan( float a ) {
 458         return tanf( a );
 459 }
 460
 461 ID_INLINE float idMath::Tan16( float a ) {
 462         float s;
 463         bool reciprocal;
 464
 465         if ( ( a < 0.0f ) || ( a >= PI ) ) {
 466                 a -= floorf( a / PI ) * PI;
 467         }
 468 #if 1
 469         if ( a < HALF_PI ) {
 470                 if ( a > ONEFOURTH_PI ) {
 471                         a = HALF_PI - a;
 472                         reciprocal = true;
 473                 } else {
 474                         reciprocal = false;
 475                 }
 476         } else {
 477                 if ( a > HALF_PI + ONEFOURTH_PI ) {
 478                         a = a - PI;
 479                         reciprocal = false;
 480                 } else {
 481                         a = HALF_PI - a;
 482                         reciprocal = true;
 483                 }
 484         }
 485 #else
 486         a = HALF_PI - a;
 487         if ( fabs( a ) >= ONEFOURTH_PI ) {
 488                 a = ( ( a < 0.0f ) ? -HALF_PI : HALF_PI ) - a;
 489                 reciprocal = false;
 490         } else {
 491                 reciprocal = true;
 492         }
 493 #endif
 494         s = a * a;
 495         s = a * ( ( ( ( ( ( 9.5168091e-03f * s + 2.900525e-03f ) * s + 2.45650893e-02f ) * s + 5.33740603e-02f ) * s + 1.333923995e-01f ) * s + 3.333314036e-01f ) * s + 1.0f );
 496         if ( reciprocal ) {
 497                 return 1.0f / s;
 498         } else {
 499                 return s;
 500         }
 501 }
 502
 503 ID_INLINE double idMath::Tan64( float a ) {
 504         return tan( a );
 505 }
 506
 507 ID_INLINE float idMath::ASin( float a ) {
 508         if ( a <= -1.0f ) {
 509                 return -HALF_PI;
 510         }
 511         if ( a >= 1.0f ) {
 512                 return HALF_PI;
 513         }
 514         return asinf( a );
 515 }
 516
 517 ID_INLINE float idMath::ASin16( float a ) {
 518         if ( FLOATSIGNBITSET( a ) ) {
 519                 if ( a <= -1.0f ) {
 520                         return -HALF_PI;
 521                 }
 522                 a = fabs( a );
 523                 return ( ( ( -0.0187293f * a + 0.0742610f ) * a - 0.2121144f ) * a + 1.5707288f ) * sqrt( 1.0f - a ) - HALF_PI;
 524         } else {
 525                 if ( a >= 1.0f ) {
 526                         return HALF_PI;
 527                 }
 528                 return HALF_PI - ( ( ( -0.0187293f * a + 0.0742610f ) * a - 0.2121144f ) * a + 1.5707288f ) * sqrt( 1.0f - a );
 529         }
 530 }
 531
 532 ID_INLINE double idMath::ASin64( float a ) {
 533         if ( a <= -1.0f ) {
 534                 return -HALF_PI;
 535         }
 536         if ( a >= 1.0f ) {
 537                 return HALF_PI;
 538         }
 539         return asin( a );
 540 }
 541
 542 ID_INLINE float idMath::ACos( float a ) {
 543         if ( a <= -1.0f ) {
 544                 return PI;
 545         }
 546         if ( a >= 1.0f ) {
 547                 return 0.0f;
 548         }
 549         return acosf( a );
 550 }
 551
 552 ID_INLINE float idMath::ACos16( float a ) {
 553         if ( FLOATSIGNBITSET( a ) ) {
 554                 if ( a <= -1.0f ) {
 555                         return PI;
 556                 }
 557                 a = fabs( a );
 558                 return PI - ( ( ( -0.0187293f * a + 0.0742610f ) * a - 0.2121144f ) * a + 1.5707288f ) * sqrt( 1.0f - a );
 559         } else {
 560                 if ( a >= 1.0f ) {
 561                         return 0.0f;
 562                 }
 563                 return ( ( ( -0.0187293f * a + 0.0742610f ) * a - 0.2121144f ) * a + 1.5707288f ) * sqrt( 1.0f - a );
 564         }
 565 }
 566
 567 ID_INLINE double idMath::ACos64( float a ) {
 568         if ( a <= -1.0f ) {
 569                 return PI;
 570         }
 571         if ( a >= 1.0f ) {
 572                 return 0.0f;
 573         }
 574         return acos( a );
 575 }
 576
 577 ID_INLINE float idMath::ATan( float a ) {
 578         return atanf( a );
 579 }
 580
 581 ID_INLINE float idMath::ATan16( float a ) {
 582         float s;
 583
 584         if ( fabs( a ) > 1.0f ) {
 585                 a = 1.0f / a;
 586                 s = a * a;
 587                 s = - ( ( ( ( ( ( ( ( ( 0.0028662257f * s - 0.0161657367f ) * s + 0.0429096138f ) * s - 0.0752896400f )
 588                                 * s + 0.1065626393f ) * s - 0.1420889944f ) * s + 0.1999355085f ) * s - 0.3333314528f ) * s ) + 1.0f ) * a;
 589                 if ( FLOATSIGNBITSET( a ) ) {
 590                         return s - HALF_PI;
 591                 } else {
 592                         return s + HALF_PI;
 593                 }
 594         } else {
 595                 s = a * a;
 596                 return ( ( ( ( ( ( ( ( ( 0.0028662257f * s - 0.0161657367f ) * s + 0.0429096138f ) * s - 0.0752896400f )
 597                         * s + 0.1065626393f ) * s - 0.1420889944f ) * s + 0.1999355085f ) * s - 0.3333314528f ) * s ) + 1.0f ) * a;
 598         }
 599 }
 600
 601 ID_INLINE double idMath::ATan64( float a ) {
 602         return atan( a );
 603 }
 604
 605 ID_INLINE float idMath::ATan( float y, float x ) {
 606         return atan2f( y, x );
 607 }
 608
 609 ID_INLINE float idMath::ATan16( float y, float x ) {
 610         float a, s;
 611
 612         if ( fabs( y ) > fabs( x ) ) {
 613                 a = x / y;
 614                 s = a * a;
 615                 s = - ( ( ( ( ( ( ( ( ( 0.0028662257f * s - 0.0161657367f ) * s + 0.0429096138f ) * s - 0.0752896400f )
 616                                 * s + 0.1065626393f ) * s - 0.1420889944f ) * s + 0.1999355085f ) * s - 0.3333314528f ) * s ) + 1.0f ) * a;
 617                 if ( FLOATSIGNBITSET( a ) ) {
 618                         return s - HALF_PI;
 619                 } else {
 620                         return s + HALF_PI;
 621                 }
 622         } else {
 623                 a = y / x;
 624                 s = a * a;
 625                 return ( ( ( ( ( ( ( ( ( 0.0028662257f * s - 0.0161657367f ) * s + 0.0429096138f ) * s - 0.0752896400f )
 626                         * s + 0.1065626393f ) * s - 0.1420889944f ) * s + 0.1999355085f ) * s - 0.3333314528f ) * s ) + 1.0f ) * a;
 627         }
 628 }
 629
 630 ID_INLINE double idMath::ATan64( float y, float x ) {
 631         return atan2( y, x );
 632 }
 633
 634 ID_INLINE float idMath::Pow( float x, float y ) {
 635         return powf( x, y );
 636 }
 637
 638 ID_INLINE float idMath::Pow16( float x, float y ) {
 639         return Exp16( y * Log16( x ) );
 640 }
 641
 642 ID_INLINE double idMath::Pow64( float x, float y ) {
 643         return pow( x, y );
 644 }
 645
 646 ID_INLINE float idMath::Exp( float f ) {
 647         return expf( f );
 648 }
 649
 650 ID_INLINE float idMath::Exp16( float f ) {
 651         int i, s, e, m, exponent;
 652         float x, x2, y, p, q;
 653
 654         x = f * 1.44269504088896340f;           // multiply with ( 1 / log( 2 ) )
 655 #if 1
 656         i = *reinterpret_cast<int *>(&x);
 657         s = ( i >> IEEE_FLT_SIGN_BIT );
 658         e = ( ( i >> IEEE_FLT_MANTISSA_BITS ) & ( ( 1 << IEEE_FLT_EXPONENT_BITS ) - 1 ) ) - IEEE_FLT_EXPONENT_BIAS;
 659         m = ( i & ( ( 1 << IEEE_FLT_MANTISSA_BITS ) - 1 ) ) | ( 1 << IEEE_FLT_MANTISSA_BITS );
 660         i = ( ( m >> ( IEEE_FLT_MANTISSA_BITS - e ) ) & ~( e >> 31 ) ) ^ s;
 661 #else
 662         i = (int) x;
 663         if ( x < 0.0f ) {
 664                 i--;
 665         }
 666 #endif
 667         exponent = ( i + IEEE_FLT_EXPONENT_BIAS ) << IEEE_FLT_MANTISSA_BITS;
 668         y = *reinterpret_cast<float *>(&exponent);
 669         x -= (float) i;
 670         if ( x >= 0.5f ) {
 671                 x -= 0.5f;
 672                 y *= 1.4142135623730950488f;    // multiply with sqrt( 2 )
 673         }
 674         x2 = x * x;
 675         p = x * ( 7.2152891511493f + x2 * 0.0576900723731f );
 676         q = 20.8189237930062f + x2;
 677         x = y * ( q + p ) / ( q - p );
 678         return x;
 679 }
 680
 681 ID_INLINE double idMath::Exp64( float f ) {
 682         return exp( f );
 683 }
 684
 685 ID_INLINE float idMath::Log( float f ) {
 686         return logf( f );
 687 }
 688
 689 ID_INLINE float idMath::Log16( float f ) {
 690         int i, exponent;
 691         float y, y2;
 692
 693         i = *reinterpret_cast<int *>(&f);
 694         exponent = ( ( i >> IEEE_FLT_MANTISSA_BITS ) & ( ( 1 << IEEE_FLT_EXPONENT_BITS ) - 1 ) ) - IEEE_FLT_EXPONENT_BIAS;
 695         i -= ( exponent + 1 ) << IEEE_FLT_MANTISSA_BITS;        // get value in the range [.5, 1>
 696         y = *reinterpret_cast<float *>(&i);
 697         y *= 1.4142135623730950488f;                                            // multiply with sqrt( 2 )
 698         y = ( y - 1.0f ) / ( y + 1.0f );
 699         y2 = y * y;
 700         y = y * ( 2.000000000046727f + y2 * ( 0.666666635059382f + y2 * ( 0.4000059794795f + y2 * ( 0.28525381498f + y2 * 0.2376245609f ) ) ) );
 701         y += 0.693147180559945f * ( (float)exponent + 0.5f );
 702         return y;
 703 }
 704
 705 ID_INLINE double idMath::Log64( float f ) {
 706         return log( f );
 707 }
 708
 709 ID_INLINE int idMath::IPow( int x, int y ) {
 710         int r; for( r = x; y > 1; y-- ) { r *= x; } return r;
 711 }
 712
 713 ID_INLINE int idMath::ILog2( float f ) {
 714         return ( ( (*reinterpret_cast<int *>(&f)) >> IEEE_FLT_MANTISSA_BITS ) & ( ( 1 << IEEE_FLT_EXPONENT_BITS ) - 1 ) ) - IEEE_FLT_EXPONENT_BIAS;
 715 }
 716
 717 ID_INLINE int idMath::ILog2( int i ) {
 718         return ILog2( (float)i );
 719 }
 720
 721 ID_INLINE int idMath::BitsForFloat( float f ) {
 722         return ILog2( f ) + 1;
 723 }
 724
 725 ID_INLINE int idMath::BitsForInteger( int i ) {
 726         return ILog2( (float)i ) + 1;
 727 }
 728
 729 ID_INLINE int idMath::MaskForFloatSign( float f ) {
 730         return ( (*reinterpret_cast<int *>(&f)) >> 31 );
 731 }
 732
 733 ID_INLINE int idMath::MaskForIntegerSign( int i ) {
 734         return ( i >> 31 );
 735 }
 736
 737 ID_INLINE int idMath::FloorPowerOfTwo( int x ) {
 738         return CeilPowerOfTwo( x ) >> 1;
 739 }
 740
 741 ID_INLINE int idMath::CeilPowerOfTwo( int x ) {
 742         x--;
 743         x |= x >> 1;
 744         x |= x >> 2;
 745         x |= x >> 4;
 746         x |= x >> 8;
 747         x |= x >> 16;
 748         x++;
 749         return x;
 750 }
 751
 752 ID_INLINE bool idMath::IsPowerOfTwo( int x ) {
 753         return ( x & ( x - 1 ) ) == 0 && x > 0;
 754 }
 755
 756 ID_INLINE int idMath::BitCount( int x ) {
 757         x -= ( ( x >> 1 ) & 0x55555555 );
 758         x = ( ( ( x >> 2 ) & 0x33333333 ) + ( x & 0x33333333 ) );
 759         x = ( ( ( x >> 4 ) + x ) & 0x0f0f0f0f );
 760         x += ( x >> 8 );
 761         return ( ( x + ( x >> 16 ) ) & 0x0000003f );
 762 }
 763
 764 ID_INLINE int idMath::BitReverse( int x ) {
 765         x = ( ( ( x >> 1 ) & 0x55555555 ) | ( ( x & 0x55555555 ) << 1 ) );
 766         x = ( ( ( x >> 2 ) & 0x33333333 ) | ( ( x & 0x33333333 ) << 2 ) );
 767         x = ( ( ( x >> 4 ) & 0x0f0f0f0f ) | ( ( x & 0x0f0f0f0f ) << 4 ) );
 768         x = ( ( ( x >> 8 ) & 0x00ff00ff ) | ( ( x & 0x00ff00ff ) << 8 ) );
 769         return ( ( x >> 16 ) | ( x << 16 ) );
 770 }
 771
 772 ID_INLINE int idMath::Abs( int x ) {
 773    int y = x >> 31;
 774    return ( ( x ^ y ) - y );
 775 }
 776
 777 ID_INLINE float idMath::Fabs( float f ) {
 778         int tmp = *reinterpret_cast<int *>( &f );
 779         tmp &= 0x7FFFFFFF;
 780         return *reinterpret_cast<float *>( &tmp );
 781 }
 782
 783 ID_INLINE float idMath::Floor( float f ) {
 784         return floorf( f );
 785 }
 786
 787 ID_INLINE float idMath::Ceil( float f ) {
 788         return ceilf( f );
 789 }
 790
 791 ID_INLINE float idMath::Rint( float f ) {
 792         return floorf( f + 0.5f );
 793 }
 794
 795 ID_INLINE int idMath::Ftoi( float f ) {
 796         return (int) f;
 797 }
 798
 799 ID_INLINE int idMath::FtoiFast( float f ) {
 800 #ifdef _WIN32
 801         int i;
 802         __asm fld               f
 803         __asm fistp             i               // use default rouding mode (round nearest)
 804         return i;
 805 #elif 0                                         // round chop (C/C++ standard)
 806         int i, s, e, m, shift;
 807         i = *reinterpret_cast<int *>(&f);
 808         s = i >> IEEE_FLT_SIGN_BIT;
 809         e = ( ( i >> IEEE_FLT_MANTISSA_BITS ) & ( ( 1 << IEEE_FLT_EXPONENT_BITS ) - 1 ) ) - IEEE_FLT_EXPONENT_BIAS;
 810         m = ( i & ( ( 1 << IEEE_FLT_MANTISSA_BITS ) - 1 ) ) | ( 1 << IEEE_FLT_MANTISSA_BITS );
 811         shift = e - IEEE_FLT_MANTISSA_BITS;
 812         return ( ( ( ( m >> -shift ) | ( m << shift ) ) & ~( e >> 31 ) ) ^ s ) - s;
 813 //#elif defined( __i386__ )
 814 #elif 0
 815         int i = 0;
 816         __asm__ __volatile__ (
 817                                                   "fld %1\n" \
 818                                                   "fistp %0\n" \
 819                                                   : "=m" (i) \
 820                                                   : "m" (f) );
 821         return i;
 822 #else
 823         return (int) f;
 824 #endif
 825 }
 826
 827 ID_INLINE unsigned long idMath::Ftol( float f ) {
 828         return (unsigned long) f;
 829 }
 830
 831 ID_INLINE unsigned long idMath::FtolFast( float f ) {
 832 #ifdef _WIN32
 833         // FIXME: this overflows on 31bits still .. same as FtoiFast
 834         unsigned long i;
 835         __asm fld               f
 836         __asm fistp             i               // use default rouding mode (round nearest)
 837         return i;
 838 #elif 0                                         // round chop (C/C++ standard)
 839         int i, s, e, m, shift;
 840         i = *reinterpret_cast<int *>(&f);
 841         s = i >> IEEE_FLT_SIGN_BIT;
 842         e = ( ( i >> IEEE_FLT_MANTISSA_BITS ) & ( ( 1 << IEEE_FLT_EXPONENT_BITS ) - 1 ) ) - IEEE_FLT_EXPONENT_BIAS;
 843         m = ( i & ( ( 1 << IEEE_FLT_MANTISSA_BITS ) - 1 ) ) | ( 1 << IEEE_FLT_MANTISSA_BITS );
 844         shift = e - IEEE_FLT_MANTISSA_BITS;
 845         return ( ( ( ( m >> -shift ) | ( m << shift ) ) & ~( e >> 31 ) ) ^ s ) - s;
 846 //#elif defined( __i386__ )
 847 #elif 0
 848         // for some reason, on gcc I need to make sure i == 0 before performing a fistp
 849         int i = 0;
 850         __asm__ __volatile__ (
 851                                                   "fld %1\n" \
 852                                                   "fistp %0\n" \
 853                                                   : "=m" (i) \
 854                                                   : "m" (f) );
 855         return i;
 856 #else
 857         return (unsigned long) f;
 858 #endif
 859 }
 860
 861 ID_INLINE signed char idMath::ClampChar( int i ) {
 862         if ( i < -128 ) {
 863                 return -128;
 864         }
 865         if ( i > 127 ) {
 866                 return 127;
 867         }
 868         return i;
 869 }
 870
 871 ID_INLINE signed short idMath::ClampShort( int i ) {
 872         if ( i < -32768 ) {
 873                 return -32768;
 874         }
 875         if ( i > 32767 ) {
 876                 return 32767;
 877         }
 878         return i;
 879 }
 880
 881 ID_INLINE int idMath::ClampInt( int min, int max, int value ) {
 882         if ( value < min ) {
 883                 return min;
 884         }
 885         if ( value > max ) {
 886                 return max;
 887         }
 888         return value;
 889 }
 890
 891 ID_INLINE float idMath::ClampFloat( float min, float max, float value ) {
 892         if ( value < min ) {
 893                 return min;
 894         }
 895         if ( value > max ) {
 896                 return max;
 897         }
 898         return value;
 899 }
 900
 901 ID_INLINE float idMath::AngleNormalize360( float angle ) {
 902         if ( ( angle >= 360.0f ) || ( angle < 0.0f ) ) {
 903                 angle -= floor( angle / 360.0f ) * 360.0f;
 904         }
 905         return angle;
 906 }
 907
 908 ID_INLINE float idMath::AngleNormalize180( float angle ) {
 909         angle = AngleNormalize360( angle );
 910         if ( angle > 180.0f ) {
 911                 angle -= 360.0f;
 912         }
 913         return angle;
 914 }
 915
 916 ID_INLINE float idMath::AngleDelta( float angle1, float angle2 ) {
 917         return AngleNormalize180( angle1 - angle2 );
 918 }
 919
 920 ID_INLINE int idMath::FloatHash( const float *array, const int numFloats ) {
 921         int i, hash = 0;
 922         const int *ptr;
 923
 924         ptr = reinterpret_cast<const int *>( array );
 925         for ( i = 0; i < numFloats; i++ ) {
 926                 hash ^= ptr[i];
 927         }
 928         return hash;
 929 }
 930
 931 #endif /* !__MATH_MATH_H__ */