C:/temp/src/j2k/Deprecated/Complex.old/COMPLEX.CPP Source File

00001 #ifndef __J2K__JComplex_CPP__
00002 #define __J2K__JComplex_CPP__
00003 
00004 #include <j2k/Fred/Math/Complex.hpp>
00005 
00006 /////////////////////////////////////////////////////////////////////////////
00007 /// Scaling operations                                                    ///
00008 /////////////////////////////////////////////////////////////////////////////
00009 inline void JComplex::operator*=(const CplxDouble r) {
00010   Real = Real * r;
00011   Imag = Imag * r;
00012 }
00013 
00014 inline void JComplex::operator/=(const CplxDouble r) {
00015   Real = Real / r;
00016   Imag = Imag / r;
00017 }
00018 
00019 /////////////////////////////////////////////////////////////////////////////
00020 /// Multiplication and Division operations                                ///
00021 /////////////////////////////////////////////////////////////////////////////
00022 JComplex& operator*=( JComplex& x, const JComplex& y ) {
00023    // Reminder: (a+bj)(c+dj) = (ac + bcj + adj + bdjj)
00024    // Result:   (ac - bd) + ( bc + ad )j
00025    // jj = -1, by definition.
00026 
00027    CplxNumber Temp_Re = (x.Real * y.Real) - (x.Imag * y.Imag);
00028    CplxNumber Temp_Im = (x.Imag * y.Real) + (x.Real * y.Imag);
00029 
00030    x.Real = Temp_Re;
00031    x.Imag = Temp_Im;
00032 
00033    return (x);
00034 }
00035 
00036 JComplex& operator/=( JComplex& x, const JComplex& y ) {
00037    // Reminder: (a+bj)/(c+dj) = (a+bj)(c-dj) / [ (c+dj)(c-dj) ]
00038    // = (ac + bcj - adj - bdjj) / (cc - ddjj)
00039    // = [ (ac + bd) + ( bc - ad )j ] / [ cc + dd ]
00040    // jj = -1, by definition.
00041 
00042    // Doesn't handle Divide by Zero, but STL does it.
00043 
00044    CplxNumber Temp_delta = (y.Real * y.Real) + (y.Imag * y.Imag);
00045    CplxNumber Temp_Re    = (x.Real * y.Real) + (x.Imag * y.Imag);
00046    CplxNumber Temp_Im    = (x.Imag * y.Real) - (x.Real * y.Imag);
00047 
00048    x.Real = Temp_Re / Temp_delta;
00049    x.Imag = Temp_Im / Temp_delta;
00050 
00051    return (x);
00052 }
00053 
00054 static void JComplex::operator*=(const JComplex& c) {
00055   // Don't affect either until both are done, do new real to temp
00056   const CplxDouble tmpReal = (Real * c.Real) - (Imag * c.Imag);
00057   Imag = (Real * c.Imag) + (Imag * c.Real);
00058   Real = tmpReal;
00059 }
00060 
00061 JComplex operator*(const JComplex& left, const JComplex& right) {
00062   return JComplex (
00063    (left.Real * right.Real) - (left.Imag * right.Imag),
00064    (left.Real * right.Imag) + (left.Imag * right.Real)
00065   );
00066 }
00067 
00068 void JComplex::operator/=(const JComplex& c) {
00069  const CplxDouble Tmp1 = c.Real / c.Imag;
00070  const CplxDouble Tmp2 = c.Imag + (Tmp1 * c.Real);
00071 
00072  // Don't affect either until both are done, do new real to temp
00073  const CplxDouble TmpReal = (Real * Tmp1 + Imag) / Tmp2;
00074  Imag = (Imag * Tmp1 - Real) / Tmp2;
00075  Real = TmpReal;
00076 }
00077 
00078 JComplex operator/(const JComplex& left, const JComplex& right) {
00079  const CplxDouble Tmp1 = right.Real / right.Imag;
00080  const CplxDouble Tmp2 = right.Imag + (Tmp1 * right.Real);
00081 
00082  return JComplex(
00083   ( (left.Real * Tmp1 + left.Imag) / Tmp2 ),
00084   ( (left.Imag * Tmp1 - left.Real) / Tmp2 )
00085  );
00086 }
00087 
00088 /////////////////////////////////////////////////////////////////////////////
00089 /// Utilities                                                             ///
00090 /////////////////////////////////////////////////////////////////////////////
00091    CplxNumber Real(const CplxNumber& x) {                                  
00092       return (Real = x);                                                   
00093    }                                                                       
00094                                                                            
00095    CplxNumber Imag(const CplxNumber& y) {                                  
00096       return (Imag = y);                                                   
00097    }                                                                       
00098                                                                            
00099    CplxNumber Real() const {                                               
00100       return (Real);                                                       
00101    }                                                                       
00102                                                                            
00103    CplxNumber Imag() const {                                               
00104       return (Imag);                                                       
00105    }                                                                       
00106                                                                            
00107    // Arithmetic Operations
00108    JComplex& operator=( const CplxNumber& x ) {                            
00109       Real = x;                                                            
00110       Imag = 0;                                                            
00111       return (*this);                                                      
00112    }                                                                       
00113                                                                            
00114    JComplex& operator=( const JComplex& c ) {
00115       Real = c.Real;                                                       
00116       Imag = c.Imag;                                                       
00117       return (*this);                                                      
00118    }                                                                       
00119                                                                            
00120    JComplex& operator+=(const CplxNumber& x ) {                            
00121       Real = Real + x;                                                     
00122       return (*this);                                                      
00123    }                                                                       
00124                                                                            
00125    JComplex& operator+=( JComplex& x, const JComplex& y ) {                
00126      x.Real = x.Real + y.Real;
00127      x.Imag = y.Imag + y.Imag;
00128      return (x);                                                           
00129    }                                                                       
00130                                                                            
00131    JComplex& operator-=(const CplxNumber& x) {                             
00132       Real = Real - x;                                                     
00133       return (*this);                                                      
00134    }                                                                       
00135                                                                            
00136    JComplex& operator-=( JComplex& x, const JComplex& y ) {                
00137      x.Real = x.Real - y.Real;
00138      x.Imag = y.Imag - y.Imag;
00139      return (x);                                                           
00140    }                                                                       
00141                                                                            
00142    JComplex& operator*=(const CplxNumber& x) {                             
00143       Real = Real * x;                                                     
00144       Imag = Imag * x;                                                     
00145       return (*this);                                                      
00146    }                                                                       
00147                                                                            
00148    JComplex& operator/=(const CplxNumber& x) {                             
00149       Real = Real / x;                                                     
00150       Imag = Imag / x;                                                     
00151       return (*this);                                                      
00152    }                                                                       
00153                                                                            
00154 protected:                                                                 
00155    CplxNumber Real;                                                        
00156    CplxNumber Imag;                                                        
00157                                                                            
00158 public:                                                                    
00159                                                                            
00160   static BOOL Compare( const JComplex& l, const JComplex& r ) {            
00161     return ( (l.Real == r.Real) && (l.Imag == r.Imag) );                   
00162   }                                                                        
00163                                                                            
00164   inline friend BOOL operator==( const JComplex& l, const JComplex& r ) {  
00165     return Compare( l, r );                                                
00166   }                                                                        
00167                                                                            
00168   inline friend BOOL operator!=( const Number& l, const Number& r ) {      
00169     return !Compare( l, r );                                               
00170   }                                                                        
00171                                                                            
00172                                                                            
00173 
00174    JComplex& operator*=( JComplex& x, const JComplex& y ) {
00175      // Reminder: (a+bj)(c+dj) = (ac + bcj + adj + bdjj)
00176      // Result:   (ac - bd) + ( bc + ad )j
00177      // jj = -1, by definition.
00178 
00179      CplxNumber Temp_Re = (x.Real * y.Real) - (x.Imag * y.Imag);
00180      CplxNumber Temp_Im = (x.Imag * y.Real) + (x.Real * y.Imag);
00181 
00182      x.Real = Temp_Re;
00183      x.Imag = Temp_Im;
00184 
00185      return (x);
00186    }
00187 
00188    JComplex& operator/=( JComplex& x, const JComplex& y ) {
00189      // Reminder: (a+bj)/(c+dj) = (a+bj)(c-dj) / [ (c+dj)(c-dj) ]
00190      // = (ac + bcj - adj - bdjj) / (cc - ddjj)
00191      // = [ (ac + bd) + ( bc - ad )j ] / [ cc + dd ]
00192      // jj = -1, by definition.
00193 
00194      // Doesn't handle Divide by Zero, but STL does it.
00195 
00196      CplxNumber Temp_delta = (y.Real * y.Real) + (y.Imag * y.Imag);
00197 
00198      CplxNumber Temp_Re = (x.Real * y.Real) + (x.Imag * y.Imag);
00199      CplxNumber Temp_Im = (x.Imag * y.Real) - (x.Real * y.Imag);
00200 
00201      x.Real = Temp_Re / Temp_delta;
00202      x.Imag = Temp_Im / Temp_delta;
00203 
00204      return (x);
00205    }
00206 
00207    static CplxNumber _Infv(CplxNumber) {                                   
00208      return ( _Inf._D );                                                   
00209    }                                                                       
00210                                                                            
00211    static bool _Isinf(CplxNumber x) {                                      
00212      double temp = (double)x;                                              
00213      return (_Dtest(&temp) == _INFCODE);                                   
00214    }                                                                       
00215                                                                            
00216    static bool _Isnan(CplxNumber x) {                                      
00217      double temp = (double)x;                                              
00218      return (_Dtest(&temp) == _NANCODE);                                   
00219    }                                                                       
00220                                                                            
00221    static CplxNumber _Nanv(CplxNumber) {                                   
00222      return (_Nan._D);                                                     
00223    }                                                                       
00224 
00225 /* STL content
00226 // It use Complex double as temp to verify for NaN, Inf
00227 // Might use Double for that.
00228 {
00229    if ( _Isnan( Temp_re ) || _Isnan( Temp_im ) ) {
00230      x.Real( _Nanv( Temp_re )), x.Imag( x.Real );
00231    } else if (
00232    (Temp_im < 0 ? -Temp_im : +Temp_im) < (Temp_re < 0 ? -Temp_re : +Temp_re)
00233    ) {
00234      CplxNumber _Wr = Temp_im / Temp_re;
00235      CplxNumber _Wd = Temp_re + _Wr * Temp_im;
00236 
00237      if (_Isnan(_Wd) || _Wd == 0) {
00238        x.Real(_Nanv(Temp_re)), x.Imag(x.Real());
00239      } else {
00240        CplxNumber temp = (x.Real() + x.Imag() * _Wr) / _Wd;
00241        x.Imag((x.Imag() - x.Real() * _Wr) / _Wd);
00242        x.Real(temp);
00243      }
00244 
00245    } else if (Temp_im == 0) {
00246       x.Real(_Nanv(Temp_re)), x.Imag(x.Real());
00247    } else {
00248       CplxNumber _Wr = Temp_re / Temp_im;
00249       CplxNumber _Wd = Temp_im + _Wr * Temp_re;
00250 
00251       if (_Isnan(_Wd) || _Wd == 0) {
00252          x.Real(_Nanv(Temp_re)), x.Imag(x.Real());
00253       } else {
00254          CplxNumber temp = (x.Real() * _Wr + x.Imag()) / _Wd;
00255          x.Imag((x.Imag() * _Wr - x.Real()) / _Wd);
00256          x.Real(temp);
00257       }
00258    }
00259    return (x);
00260 }
00261 */
00262 
00263 /*
00264 
00265 // TEMPLATE FUNCTION operator>>
00266 template<class Exponent, class _Tr, class _U> inline
00267    basic_istream<Exponent, _Tr>&  operator>>(
00268       basic_istream<Exponent, _Tr>& Im, JComplex& x)
00269 {
00270    Exponent _Ch;
00271    long double Real, Imag;
00272    if (Im >> _Ch && _Ch != '(')
00273       Im.putback(_Ch), Im >> Real, Imag = 0;
00274    else if (Im >> Real >> _Ch && _Ch != ',')
00275       if (_Ch == ')')
00276          Imag = 0;
00277       else
00278          Im.putback(_Ch), Im.setstate(ios_base::failbit);
00279    else if (Im >> Imag >> _Ch && _Ch != ')')
00280          Im.putback(_Ch), Im.setstate(ios_base::failbit);
00281    if (!Im.fail())
00282       x = JComplex((_U)Real, (_U)Imag);
00283    return (Im);
00284 }
00285 
00286 // TEMPLATE FUNCTION operator<<
00287 template<class Exponent, class _Tr, class _U> inline
00288    basic_ostream<Exponent, _Tr>&  operator<<(
00289       basic_ostream<Exponent, _Tr>& _O, const JComplex& x)
00290    {basic_ostringstream<Exponent, _Tr, allocator<Exponent> > _S;
00291    _S.flags(_O.flags());
00292    _S.imbue(_O.getloc());
00293    _S.precision(_O.precision());
00294    _S << '(' << Real(x) << ',' << Imag(x) << ')';
00295    return (_O << _S.str().c_str()); }
00296  #include <xcomplex>
00297 */
00298 
00299 #endif