source: projectionDesigner/tag/ProjectionDesigner_1.1.5/projdesigner/include/gmtl/Math.h

Last change on this file was 2, checked in by Torben Dannhauer, 14 years ago
File size: 15.0 KB
Line 
1/************************************************************** ggt-head beg
2 *
3 * GGT: Generic Graphics Toolkit
4 *
5 * Original Authors:
6 *   Allen Bierbaum
7 *
8 * -----------------------------------------------------------------
9 * File:          Math.h,v
10 * Date modified: 2005/06/23 21:13:28
11 * Version:       1.40
12 * -----------------------------------------------------------------
13 *
14 *********************************************************** ggt-head end */
15/*************************************************************** ggt-cpr beg
16*
17* GGT: The Generic Graphics Toolkit
18* Copyright (C) 2001,2002 Allen Bierbaum
19*
20* This library is free software; you can redistribute it and/or
21* modify it under the terms of the GNU Lesser General Public
22* License as published by the Free Software Foundation; either
23* version 2.1 of the License, or (at your option) any later version.
24*
25* This library is distributed in the hope that it will be useful,
26* but WITHOUT ANY WARRANTY; without even the implied warranty of
27* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
28* Lesser General Public License for more details.
29*
30* You should have received a copy of the GNU Lesser General Public
31* License along with this library; if not, write to the Free Software
32* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
33*
34 ************************************************************ ggt-cpr end */
35#ifndef _GMTL_MATH_H_
36#define _GMTL_MATH_H_
37
38#include <math.h>
39#include <stdlib.h>
40#include <gmtl/Defines.h>
41#include <gmtl/Util/Assert.h>
42
43namespace gmtl
44{
45
46/** Base class for Rotation orders
47 *  @ingroup Defines
48 * @see XYZ, ZYX, ZXY
49 */
50struct RotationOrderBase { enum { IS_ROTORDER = 1 }; };
51
52/** XYZ Rotation order
53 *  @ingroup Defines */
54struct XYZ : public RotationOrderBase { enum { ID = 0 }; };
55
56/** ZYX Rotation order
57 *  @ingroup Defines */
58struct ZYX : public RotationOrderBase { enum { ID = 1 }; };
59
60/** ZXY Rotation order
61 *  @ingroup Defines */
62struct ZXY : public RotationOrderBase { enum { ID = 2 }; };
63
64namespace Math
65{
66   /** @ingroup Math
67    *  @name Mathematical constants
68    *  @{
69    */
70   const float TWO_PI = 6.28318530717958647692f;
71   const float PI = 3.14159265358979323846f; //3.14159265358979323846264338327950288419716939937510;
72   const float PI_OVER_2 = 1.57079632679489661923f;
73   const float PI_OVER_4 = 0.78539816339744830962f;
74   /** @} */
75
76   /** @ingroup Math
77    *  @name C Math Abstraction
78    *  @{
79    */
80//----------------------------------------------------------------------------
81template <typename T>
82inline T abs( T iValue )
83{
84    return T( iValue >= ((T)0) ? iValue : -iValue );
85}
86
87inline float abs(float iValue)
88{  return fabsf(iValue); }
89inline double abs(double iValue)
90{  return fabs(iValue); }
91inline int abs(int iValue)
92{  return ::abs(iValue); }
93inline long abs(long iValue)
94{  return labs(iValue); }
95
96//----------------------------------------------------------------------------
97template <typename T>
98inline T ceil( T fValue );
99inline float ceil( float fValue )
100{
101#ifdef NO_CEILF
102   return float(::ceil(fValue));
103#else
104   return float( ::ceilf( fValue ) );
105#endif
106}
107inline double ceil( double fValue )
108{
109    return double( ::ceil( fValue ) );
110}
111//----------------------------------------------------------------------------
112template <typename T>
113inline T floor( T fValue ); // why do a floor of int?  shouldn't compile...
114inline float floor( float fValue )
115{
116#ifdef NO_FLOORF
117   return float(::floor(fValue));
118#else
119   return float( ::floorf( fValue ) );
120#endif
121}
122inline double floor( double fValue )
123{
124    return double( ::floor( fValue ) );
125}
126//----------------------------------------------------------------------------
127template <typename T>
128inline int sign( T iValue )
129{
130   if (iValue > T(0))
131   {
132      return 1;
133   }
134   else
135   {
136      if (iValue < T(0))
137      {
138         return -1;
139      }
140      else
141      {
142         return 0;
143      }
144   }
145}
146//----------------------------------------------------------------------------
147/**
148 * Clamps the given value down to zero if it is within epsilon of zero.
149 *
150 * @param value      the value to clamp
151 * @param eps        the epsilon tolerance or zero by default
152 *
153 * @return  zero if the value is close to 0, the value otherwise
154 */
155template <typename T>
156inline T zeroClamp( T value, T eps = T(0) )
157{
158   return ( (gmtl::Math::abs(value) <= eps) ? T(0) : value );
159}
160
161//----------------------------------------------------------------------------
162// don't allow non-float types, because their ret val wont be correct
163// i.e. with int, the int retval will be rounded up or down.
164//      we'd need a float retval to do it right, but we can't specialize by ret
165template <typename T>
166inline T aCos( T fValue );
167inline float aCos( float fValue )
168{
169    if ( -1.0f < fValue )
170    {
171        if ( fValue < 1.0f )
172        {
173#ifdef NO_ACOSF
174            return float(::acos(fValue));
175#else
176            return float( ::acosf( fValue ) );
177#endif
178        }
179        else
180            return 0.0f;
181    }
182    else
183    {
184        return (float)gmtl::Math::PI;
185    }
186}
187inline double aCos( double fValue )
188{
189    if ( -1.0 < fValue )
190    {
191        if ( fValue < 1.0 )
192            return double( ::acos( fValue ) );
193        else
194            return 0.0;
195    }
196    else
197    {
198        return (double)gmtl::Math::PI;
199    }
200}
201//----------------------------------------------------------------------------
202template <typename T>
203inline T aSin( T fValue );
204inline float aSin( float fValue )
205{
206    if ( -1.0f < fValue )
207    {
208        if ( fValue < 1.0f )
209        {
210#ifdef NO_ASINF
211            return float(::asin(fValue));
212#else
213            return float( ::asinf( fValue ) );
214#endif
215        }
216        else
217            return (float)-gmtl::Math::PI_OVER_2;
218    }
219    else
220    {
221        return (float)gmtl::Math::PI_OVER_2;
222    }
223}
224inline double aSin( double fValue )
225{
226    if ( -1.0 < fValue )
227    {
228        if ( fValue < 1.0 )
229            return double( ::asin( fValue ) );
230        else
231            return (double)-gmtl::Math::PI_OVER_2;
232    }
233    else
234    {
235        return (double)gmtl::Math::PI_OVER_2;
236    }
237}
238//----------------------------------------------------------------------------
239template <typename T>
240inline T aTan( T fValue );
241inline double aTan( double fValue )
242{
243    return ::atan( fValue );
244}
245inline float aTan( float fValue )
246{
247#ifdef NO_TANF
248   return float(::atan(fValue));
249#else
250   return float( ::atanf( fValue ) );
251#endif
252}
253//----------------------------------------------------------------------------
254template <typename T>
255inline T aTan2( T fY, T fX );
256inline float aTan2( float fY, float fX )
257{
258#ifdef NO_ATAN2F
259   return float(::atan2(fY, fX));
260#else
261   return float( ::atan2f( fY, fX ) );
262#endif
263}
264inline double aTan2( double fY, double fX )
265{
266    return double( ::atan2( fY, fX ) );
267}
268//----------------------------------------------------------------------------
269template <typename T>
270inline T cos( T fValue );
271inline float cos( float fValue )
272{
273#ifdef NO_COSF
274   return float(::cos(fValue));
275#else
276   return float( ::cosf( fValue ) );
277#endif
278}
279inline double cos( double fValue )
280{
281    return double( ::cos( fValue ) );
282}
283//----------------------------------------------------------------------------
284template <typename T>
285inline T exp( T fValue );
286inline float exp( float fValue )
287{
288#ifdef NO_EXPF
289   return float(::exp(fValue));
290#else
291   return float( ::expf( fValue ) );
292#endif
293}
294inline double exp( double fValue )
295{
296    return double( ::exp( fValue ) );
297}
298//----------------------------------------------------------------------------
299template <typename T>
300inline T log( T fValue );
301inline double log( double fValue )
302{
303    return double( ::log( fValue ) );
304}
305inline float log( float fValue )
306{
307#ifdef NO_LOGF
308   return float(::log(fValue));
309#else
310   return float( ::logf( fValue ) );
311#endif
312}
313//----------------------------------------------------------------------------
314inline double pow( double fBase, double fExponent)
315{
316    return double( ::pow( fBase, fExponent ) );
317}
318inline float pow( float fBase, float fExponent)
319{
320#ifdef NO_POWF
321   return float(::pow(fBase, fExponent));
322#else
323   return float( ::powf( fBase, fExponent ) );
324#endif
325}
326//----------------------------------------------------------------------------
327template <typename T>
328inline T sin( T fValue );
329inline double sin( double fValue )
330{
331    return double( ::sin( fValue ) );
332}
333inline float sin( float fValue )
334{
335#ifdef NO_SINF
336   return float(::sin(fValue));
337#else
338   return float( ::sinf( fValue ) );
339#endif
340}
341//----------------------------------------------------------------------------
342template <typename T>
343inline T tan( T fValue );
344inline double tan( double fValue )
345{
346    return double( ::tan( fValue ) );
347}
348inline float tan( float fValue )
349{
350#ifdef NO_TANF
351   return float(::tan(fValue));
352#else
353   return float( ::tanf( fValue ) );
354#endif
355}
356//----------------------------------------------------------------------------
357template <typename T>
358inline T sqr( T fValue )
359{
360    return T( fValue * fValue );
361}
362//----------------------------------------------------------------------------
363template <typename T>
364inline T sqrt( T fValue )
365{
366#ifdef NO_SQRTF
367   return T(::sqrt(((float)fValue)));
368#else
369   return T( ::sqrtf( ((float)fValue) ) );
370#endif
371}
372inline double sqrt( double fValue )
373{
374    return double( ::sqrt( fValue ) );
375}
376
377/** Fast inverse square root.
378 */
379inline float fastInvSqrt(float x)
380{
381   const float xhalf(0.5f*x);
382   long i = *(long*)&x;
383   i = 0x5f3759df - (i>>1);    // This hides a good amount of math
384   x = *(float*)&i;
385   x = x*(1.5f - xhalf*x*x);   // Repeat for more accuracy
386   return x;
387}
388
389inline float fastInvSqrt2(float x)
390{
391   const float xhalf(0.5f*x);
392   long i = *(long*)&x;
393   i = 0x5f3759df - (i>>1);    // This hides a good amount of math
394   x = *(float*)&i;
395   x = x*(1.5f - xhalf*x*x);   // Repeat for more accuracy
396   x = x*(1.5f - xhalf*x*x);
397   return x;
398}
399
400inline float fastInvSqrt3(float x)
401{
402   const float xhalf(0.5f*x);
403   long i = *(long*)&x;
404   i = 0x5f3759df - (i>>1);    // This hides a good amount of math
405   x = *(float*)&i;
406   x = x*(1.5f - xhalf*x*x);   // Repeat for more accuracy
407   x = x*(1.5f - xhalf*x*x);
408   x = x*(1.5f - xhalf*x*x);
409   return x;
410}
411
412
413//----------------------------------------------------------------------------
414/**
415 * Seeds the pseudorandom number generator with the given seed.
416 *
417 * @param seed  the seed for the pseudorandom number generator.
418 */
419inline void seedRandom(unsigned int seed)
420{
421   ::srand(seed);
422}
423
424/** get a random number between 0 and 1
425 * @post returns number between 0 and 1
426 */
427inline float unitRandom()
428{
429   return float(::rand())/float(RAND_MAX);
430}
431
432/** return a random number between x1 and x2
433 * RETURNS: random number between x1 and x2
434 */
435inline float rangeRandom( float x1, float x2 )
436{
437   float r = gmtl::Math::unitRandom();
438   float size = x2 - x1;
439   return float( r * size + x1 );
440}
441
442/*
443float SymmetricRandom ()
444{
445    return 2.0*float(rand())/float(RAND_MAX) - 1.0;
446}
447*/
448//----------------------------------------------------------------------------
449
450inline float deg2Rad( float fVal )
451{
452   return float( fVal * (float)(gmtl::Math::PI/180.0) );
453}
454inline double deg2Rad( double fVal )
455{
456   return double( fVal * (double)(gmtl::Math::PI/180.0) );
457}
458
459inline float rad2Deg( float fVal )
460{
461   return float( fVal * (float)(180.0/gmtl::Math::PI) );
462}
463inline double rad2Deg( double fVal )
464{
465   return double( fVal * (double)(180.0/gmtl::Math::PI) );
466}
467//----------------------------------------------------------------------------
468
469/** Is almost equal?
470 * test for equality within some tolerance...
471 * @PRE: tolerance must be >= 0
472 */
473template <class T>
474inline bool isEqual( const T& a, const T& b, const T& tolerance )
475{
476   gmtlASSERT( tolerance >= (T)0 );
477   return bool( gmtl::Math::abs( a - b ) <= tolerance );
478}
479//----------------------------------------------------------------------------
480/** cut off the digits after the decimal place */
481template <class T>
482inline T trunc( T val )
483{
484   return T( (val < ((T)0)) ? gmtl::Math::ceil( val ) : gmtl::Math::floor( val ) );
485}
486/** round to nearest integer */
487template <class T>
488inline T round( T p )
489{
490   return T( gmtl::Math::floor( p + (T)0.5 ) );
491}
492//----------------------------------------------------------------------------
493/** min returns the minimum of 2 values */
494template <class T>
495inline T Min( const T& x, const T& y )
496{
497   return ( x > y ) ? y : x;
498}
499/** min returns the minimum of 3 values */
500template <class T>
501inline T Min( const T& x, const T& y, const T& z )
502{
503   return Min( gmtl::Math::Min( x, y ), z );
504}
505/** min returns the minimum of 4 values */
506template <class T>
507inline T Min( const T& w, const T& x, const T& y, const T& z )
508{
509   return gmtl::Math::Min( gmtl::Math::Min( w, x ), gmtl::Math::Min( y, z ) );
510}
511
512/** max returns the maximum of 2 values */
513template <class T>
514inline T Max( const T& x, const T& y )
515{
516   return ( x > y ) ? x : y;
517}
518/** max returns the maximum of 3 values */
519template <class T>
520inline T Max( const T& x, const T& y, const T& z )
521{
522   return Max( gmtl::Math::Max( x, y ), z );
523}
524/** max returns the maximum of 4 values */
525template <class T>
526inline T Max( const T& w, const T& x, const T& y, const T& z )
527{
528   return gmtl::Math::Max( gmtl::Math::Max( w, x ), gmtl::Math::Max( y, z ) );
529}
530//----------------------------------------------------------------------------
531/** Compute the factorial.
532 * give - an object who's type has operator++, operator=, operator<=, and operator*= defined.
533 *        it should be a single valued scalar type such as an int, float, double etc....
534 * NOTE: This could be faster with a lookup table, but then wouldn't work templated : kevin
535 */
536template<class T>
537inline T factorial(T rhs)
538{
539   T lhs = (T)1;
540
541   for( T x = (T)1; x <= rhs; ++x )
542   {
543      lhs *= x;
544   }
545
546   return lhs;
547}
548/** @} */
549
550/**
551 * clamp "number" to a range between lo and hi
552 */
553template <class T>
554inline T clamp( T number, T lo, T hi )
555{
556   if (number > hi) number = hi;
557   else if (number < lo) number = lo;
558   return number;
559}
560
561/** @ingroup Interp
562 * @name Scalar type interpolation (for doubles, floats, etc...)
563 * @{
564 */
565
566/** Linear Interpolation between number [a] and [b].
567 *  lerp=0.0 returns a, lerp=1.0 returns b
568 *  @pre use double or float only...
569 */
570template <class T, typename U>
571inline void lerp( T& result, const U& lerp, const T& a, const T& b )
572{
573    T size = b - a;
574    result = ((U)a) + (((U)size) * lerp);
575}
576/** @} */
577
578/**
579 * Uses the quadratic formula to compute the 2 roots of the given 2nd degree
580 * polynomial in the form of Ax^2 + Bx + C.
581 *
582 * @param r1   set to the first root
583 * @param r2   set to the second root
584 * @param a    the coefficient to x^2
585 * @param b    the coefficient to x^1
586 * @param c    the coefficient to x^0
587 *
588 * @return  true if both r1 and r2 are real; false otherwise
589 */
590template <class T>
591inline bool quadraticFormula(T& r1, T& r2, const T& a, const T& b, const T& c)
592{
593   const T q = b*b - T(4)*a*c;
594
595   // the result has real roots
596   if (q >= 0)
597   {
598      const T sq = gmtl::Math::sqrt(q);
599      const T d = T(1) / (T(2) * a);
600      r1 = (-b + sq) * d;
601      r2 = (-b - sq) * d;
602      return true;
603   }
604   // the result has complex roots
605   else
606   {
607      return false;
608   }
609}
610
611} // end namespace Math
612} // end namespace gmtl
613
614#endif
Note: See TracBrowser for help on using the repository browser.