Added an #ifdef to block out a PowerPC-specific header on Mac OS X.
[icculus/iodoom3.git] / 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__ */