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