neo/idlib/math/Complex.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_COMPLEX_H__
  30 #define __MATH_COMPLEX_H__
  31
  32 /*
  33 ===============================================================================
  34
  35   Complex number
  36
  37 ===============================================================================
  38 */
  39
  40 class idComplex {
  41 public:
  42         float                           r;              // real part
  43         float                           i;              // imaginary part
  44
  45                                                 idComplex( void );
  46                                                 idComplex( const float r, const float i );
  47
  48         void                            Set( const float r, const float i );
  49         void                            Zero( void );
  50
  51         float                           operator[]( int index ) const;
  52         float &                         operator[]( int index );
  53
  54         idComplex                       operator-() const;
  55         idComplex &                     operator=( const idComplex &a );
  56
  57         idComplex                       operator*( const idComplex &a ) const;
  58         idComplex                       operator/( const idComplex &a ) const;
  59         idComplex                       operator+( const idComplex &a ) const;
  60         idComplex                       operator-( const idComplex &a ) const;
  61
  62         idComplex &                     operator*=( const idComplex &a );
  63         idComplex &                     operator/=( const idComplex &a );
  64         idComplex &                     operator+=( const idComplex &a );
  65         idComplex &                     operator-=( const idComplex &a );
  66
  67         idComplex                       operator*( const float a ) const;
  68         idComplex                       operator/( const float a ) const;
  69         idComplex                       operator+( const float a ) const;
  70         idComplex                       operator-( const float a ) const;
  71
  72         idComplex &                     operator*=( const float a );
  73         idComplex &                     operator/=( const float a );
  74         idComplex &                     operator+=( const float a );
  75         idComplex &                     operator-=( const float a );
  76
  77         friend idComplex        operator*( const float a, const idComplex &b );
  78         friend idComplex        operator/( const float a, const idComplex &b );
  79         friend idComplex        operator+( const float a, const idComplex &b );
  80         friend idComplex        operator-( const float a, const idComplex &b );
  81
  82         bool                            Compare( const idComplex &a ) const;                                            // exact compare, no epsilon
  83         bool                            Compare( const idComplex &a, const float epsilon ) const;       // compare with epsilon
  84         bool                            operator==(     const idComplex &a ) const;                                             // exact compare, no epsilon
  85         bool                            operator!=(     const idComplex &a ) const;                                             // exact compare, no epsilon
  86
  87         idComplex                       Reciprocal( void ) const;
  88         idComplex                       Sqrt( void ) const;
  89         float                           Abs( void ) const;
  90
  91         int                                     GetDimension( void ) const;
  92
  93         const float *           ToFloatPtr( void ) const;
  94         float *                         ToFloatPtr( void );
  95         const char *            ToString( int precision = 2 ) const;
  96 };
  97
  98 extern idComplex complex_origin;
  99 #define complex_zero complex_origin
 100
 101 ID_INLINE idComplex::idComplex( void ) {
 102 }
 103
 104 ID_INLINE idComplex::idComplex( const float r, const float i ) {
 105         this->r = r;
 106         this->i = i;
 107 }
 108
 109 ID_INLINE void idComplex::Set( const float r, const float i ) {
 110         this->r = r;
 111         this->i = i;
 112 }
 113
 114 ID_INLINE void idComplex::Zero( void ) {
 115         r = i = 0.0f;
 116 }
 117
 118 ID_INLINE float idComplex::operator[]( int index ) const {
 119         assert( index >= 0 && index < 2 );
 120         return ( &r )[ index ];
 121 }
 122
 123 ID_INLINE float& idComplex::operator[]( int index ) {
 124         assert( index >= 0 && index < 2 );
 125         return ( &r )[ index ];
 126 }
 127
 128 ID_INLINE idComplex idComplex::operator-() const {
 129         return idComplex( -r, -i );
 130 }
 131
 132 ID_INLINE idComplex &idComplex::operator=( const idComplex &a ) {
 133         r = a.r;
 134         i = a.i;
 135         return *this;
 136 }
 137
 138 ID_INLINE idComplex idComplex::operator*( const idComplex &a ) const {
 139         return idComplex( r * a.r - i * a.i, i * a.r + r * a.i );
 140 }
 141
 142 ID_INLINE idComplex idComplex::operator/( const idComplex &a ) const {
 143         float s, t;
 144         if ( idMath::Fabs( a.r ) >= idMath::Fabs( a.i ) ) {
 145                 s = a.i / a.r;
 146                 t = 1.0f / ( a.r + s * a.i );
 147                 return idComplex( ( r + s * i ) * t, ( i - s * r ) * t );
 148         } else {
 149                 s = a.r / a.i;
 150                 t = 1.0f / ( s * a.r + a.i );
 151                 return idComplex( ( r * s + i ) * t, ( i * s - r ) * t );
 152         }
 153 }
 154
 155 ID_INLINE idComplex idComplex::operator+( const idComplex &a ) const {
 156         return idComplex( r + a.r, i + a.i );
 157 }
 158
 159 ID_INLINE idComplex idComplex::operator-( const idComplex &a ) const {
 160         return idComplex( r - a.r, i - a.i );
 161 }
 162
 163 ID_INLINE idComplex &idComplex::operator*=( const idComplex &a ) {
 164         *this = idComplex( r * a.r - i * a.i, i * a.r + r * a.i );
 165         return *this;
 166 }
 167
 168 ID_INLINE idComplex &idComplex::operator/=( const idComplex &a ) {
 169         float s, t;
 170         if ( idMath::Fabs( a.r ) >= idMath::Fabs( a.i ) ) {
 171                 s = a.i / a.r;
 172                 t = 1.0f / ( a.r + s * a.i );
 173                 *this = idComplex( ( r + s * i ) * t, ( i - s * r ) * t );
 174         } else {
 175                 s = a.r / a.i;
 176                 t = 1.0f / ( s * a.r + a.i );
 177                 *this = idComplex( ( r * s + i ) * t, ( i * s - r ) * t );
 178         }
 179         return *this;
 180 }
 181
 182 ID_INLINE idComplex &idComplex::operator+=( const idComplex &a ) {
 183         r += a.r;
 184         i += a.i;
 185         return *this;
 186 }
 187
 188 ID_INLINE idComplex &idComplex::operator-=( const idComplex &a ) {
 189         r -= a.r;
 190         i -= a.i;
 191         return *this;
 192 }
 193
 194 ID_INLINE idComplex idComplex::operator*( const float a ) const {
 195         return idComplex( r * a, i * a );
 196 }
 197
 198 ID_INLINE idComplex idComplex::operator/( const float a ) const {
 199         float s = 1.0f / a;
 200         return idComplex( r * s, i * s );
 201 }
 202
 203 ID_INLINE idComplex idComplex::operator+( const float a ) const {
 204         return idComplex( r + a, i );
 205 }
 206
 207 ID_INLINE idComplex idComplex::operator-( const float a ) const {
 208         return idComplex( r - a, i );
 209 }
 210
 211 ID_INLINE idComplex &idComplex::operator*=( const float a ) {
 212         r *= a;
 213         i *= a;
 214         return *this;
 215 }
 216
 217 ID_INLINE idComplex &idComplex::operator/=( const float a ) {
 218         float s = 1.0f / a;
 219         r *= s;
 220         i *= s;
 221         return *this;
 222 }
 223
 224 ID_INLINE idComplex &idComplex::operator+=( const float a ) {
 225         r += a;
 226         return *this;
 227 }
 228
 229 ID_INLINE idComplex &idComplex::operator-=( const float a ) {
 230         r -= a;
 231         return *this;
 232 }
 233
 234 ID_INLINE idComplex operator*( const float a, const idComplex &b ) {
 235         return idComplex( a * b.r, a * b.i );
 236 }
 237
 238 ID_INLINE idComplex operator/( const float a, const idComplex &b ) {
 239         float s, t;
 240         if ( idMath::Fabs( b.r ) >= idMath::Fabs( b.i ) ) {
 241                 s = b.i / b.r;
 242                 t = a / ( b.r + s * b.i );
 243                 return idComplex( t, - s * t );
 244         } else {
 245                 s = b.r / b.i;
 246                 t = a / ( s * b.r + b.i );
 247                 return idComplex( s * t, - t );
 248         }
 249 }
 250
 251 ID_INLINE idComplex operator+( const float a, const idComplex &b ) {
 252         return idComplex( a + b.r, b.i );
 253 }
 254
 255 ID_INLINE idComplex operator-( const float a, const idComplex &b ) {
 256         return idComplex( a - b.r, -b.i );
 257 }
 258
 259 ID_INLINE idComplex idComplex::Reciprocal( void ) const {
 260         float s, t;
 261         if ( idMath::Fabs( r ) >= idMath::Fabs( i ) ) {
 262                 s = i / r;
 263                 t = 1.0f / ( r + s * i );
 264                 return idComplex( t, - s * t );
 265         } else {
 266                 s = r / i;
 267                 t = 1.0f / ( s * r + i );
 268                 return idComplex( s * t, - t );
 269         }
 270 }
 271
 272 ID_INLINE idComplex idComplex::Sqrt( void ) const {
 273         float x, y, w;
 274
 275         if ( r == 0.0f && i == 0.0f ) {
 276                 return idComplex( 0.0f, 0.0f );
 277         }
 278         x = idMath::Fabs( r );
 279         y = idMath::Fabs( i );
 280         if ( x >= y ) {
 281                 w = y / x;
 282                 w = idMath::Sqrt( x ) * idMath::Sqrt( 0.5f * ( 1.0f + idMath::Sqrt( 1.0f + w * w ) ) );
 283         } else {
 284                 w = x / y;
 285                 w = idMath::Sqrt( y ) * idMath::Sqrt( 0.5f * ( w + idMath::Sqrt( 1.0f + w * w ) ) );
 286         }
 287         if ( w == 0.0f ) {
 288                 return idComplex( 0.0f, 0.0f );
 289         }
 290         if ( r >= 0.0f ) {
 291                 return idComplex( w, 0.5f * i / w );
 292         } else {
 293                 return idComplex( 0.5f * y / w, ( i >= 0.0f ) ? w : -w );
 294         }
 295 }
 296
 297 ID_INLINE float idComplex::Abs( void ) const {
 298         float x, y, t;
 299         x = idMath::Fabs( r );
 300         y = idMath::Fabs( i );
 301         if ( x == 0.0f ) {
 302                 return y;
 303         } else if ( y == 0.0f ) {
 304                 return x;
 305         } else if ( x > y ) {
 306                 t = y / x;
 307                 return x * idMath::Sqrt( 1.0f + t * t );
 308         } else {
 309                 t = x / y;
 310                 return y * idMath::Sqrt( 1.0f + t * t );
 311         }
 312 }
 313
 314 ID_INLINE bool idComplex::Compare( const idComplex &a ) const {
 315         return ( ( r == a.r ) && ( i == a.i ) );
 316 }
 317
 318 ID_INLINE bool idComplex::Compare( const idComplex &a, const float epsilon ) const {
 319         if ( idMath::Fabs( r - a.r ) > epsilon ) {
 320                 return false;
 321         }
 322         if ( idMath::Fabs( i - a.i ) > epsilon ) {
 323                 return false;
 324         }
 325         return true;
 326 }
 327
 328 ID_INLINE bool idComplex::operator==( const idComplex &a ) const {
 329         return Compare( a );
 330 }
 331
 332 ID_INLINE bool idComplex::operator!=( const idComplex &a ) const {
 333         return !Compare( a );
 334 }
 335
 336 ID_INLINE int idComplex::GetDimension( void ) const {
 337         return 2;
 338 }
 339
 340 ID_INLINE const float *idComplex::ToFloatPtr( void ) const {
 341         return &r;
 342 }
 343
 344 ID_INLINE float *idComplex::ToFloatPtr( void ) {
 345         return &r;
 346 }
 347
 348 #endif /* !__MATH_COMPLEX_H__ */