//-------------------------------------------------- // vs_complex.h // /*! \brief Template class for complex numbers \class complex complex.h this routines are based on volya_complex.h (Volya at MSU) Version 1.0.1, he did \verbatim -- Some modification for speed -- Some protections from overflow and underflow -- Solving of equations of order 2,3,4: Solve2, Solve3, Solve4; \endverbatim this routine extends Volyas routine by using templates for more flexibility. Even it is tried to get as close as possible to Volyas routine, it is not the same. Where a scalar x = {int, long int, float, double, long double} occured, those spezializations are implemented . But only where it is useful eg. complex += x but not x += complex (integers are skipped) */ // // Carsten Schwarz (C.Schwarz@gsi.de) // GSI // Plankstr. 1 // 64291 Darmstadt // Germany // // 21.7.99 // //--------------------------------------------------- // $Id: complex.hxx,v 1.2 2003/05/01 13:11:48 carsten Exp $ // cvs log at end of file #ifndef __COMPLEX_CXX__ #define __COMPLEX_CXX__ // endif at end of file #ifndef small_epsilon #define small_epsilon 1e-12 #endif #ifndef PI #define PI 3.141592653589793238462643 #endif #include #include #include #include template class complex { public: complex() { re=0.0; im=0.0; } complex(scalar r,scalar i) : re(r),im(i){} public: scalar re; scalar im; }; //!specialization of complex numbers template template <> class complex { public: complex(float r=0.0, float i=0.0) : re(r), im(i) {} explicit complex(const complex& a); explicit complex(const complex& a); public: float re,im; }; //-------- //!specialization of complex numbers template template <> class complex { public: complex(double r=0.0, double i=0.0) : re(r), im(i) {} complex (const complex& a); explicit complex (const complex& a); public: double re,im; }; //-------- //!specialization of complex numbers template template <> class complex { public: complex(long double r=0.0, long double i=0.0) : re(r), im(i) {} complex (const complex& a); complex (const complex& a); public: long double re,im; }; //-------------------------------------------------------------- // += -= *= /= //-------------------------------------------------------------- // complex += complex template inline complex& operator += (complex& a, const complex& b) { a.re += b.re; a.im += b.im; return a; } // scalar += complex template inline scalar& operator += (float &a, const complex& b) { assert(fabs(b.im) < small_epsilon); a+=b.re; return a; } template inline scalar& operator += (double &a, const complex& b) { assert(fabs(b.im) < small_epsilon); a+=b.re; return a; } template inline scalar& operator += (long double &a, const complex& b) { assert(fabs(b.im) < small_epsilon); a+=b.re; return a; } // complex += scalar template inline complex& operator += (complex& a, const int &b) { a.re += b; return a; } template inline complex& operator += (complex& a, const long int &b) { a.re += b; return a; } template inline complex& operator += (complex& a, const float &b) { a.re += b; return a; } template inline complex& operator += (complex& a, const double &b) { a.re += b; return a; } template inline complex& operator += (complex& a, const long double &b) { a.re += b; return a; } // complex -= complex template inline complex& operator -= (complex& a, const complex& b) { a.re -= b.re; a.im -= b.im; return a; } // scalar -= complex template inline scalar& operator -= (float &a, const complex& b) { assert(fabs(b.im) < small_epsilon); a-=b.re; return a; } template inline scalar& operator -= (double &a, const complex& b) { assert(fabs(b.im) < small_epsilon); a-=b.re; return a; } template inline scalar& operator -= (long double &a, const complex& b) { assert(fabs(b.im) < small_epsilon); a-=b.re; return a; } // complex -= scalar template inline complex& operator -= (complex& a, const int &b) { a.re -= b; return a; } template inline complex& operator -= (complex& a, const long int &b) { a.re -= b; return a; } template inline complex& operator -= (complex& a, const float &b) { a.re -= b; return a; } template inline complex& operator -= (complex& a, const double &b) { a.re -= b; return a; } template inline complex& operator -= (complex& a, const long double &b) { a.re -= b; return a; } // complex *= complex template inline complex& operator *= (complex& a, const complex& b) { scalar t1=a.re*b.re, t2=a.im*b.im; a.im = ((a.re+a.im)*(b.re+b.im)-t1-t2); a.re = (t1 - t2); return a; } // scalar *= complex template inline scalar& operator *= (float &a, const complex& b) { assert(fabs(b.im) < small_epsilon); a *= b.re; return a; } template inline scalar& operator *= (double &a, const complex& b) { assert(fabs(b.im) < small_epsilon); a *= b.re; return a; } template inline scalar& operator *= (long double &a, const complex& b) { assert(fabs(b.im) < small_epsilon); a *= b.re; return a; } // complex *= scalar template inline complex& operator *= (complex& a, const int &b) { a.re *= b; a.im *= b; return a; } template inline complex& operator *= (complex& a, const long int &b) { a.re *= b; a.im *= b; return a; } template inline complex& operator *= (complex& a, const float &b) { a.re *= b; a.im *= b; return a; } template inline complex& operator *= (complex& a, const double &b) { a.re *= b; a.im *= b; return a; } template inline complex& operator *= (complex& a, const long double &b) { a.re *= b; a.im *= b; return a; } // complex /= complex template inline complex& operator /= (complex& a, const complex& b) { scalar t1, t2, temp; if (fabs(b.re) >= fabs(b.im)) { t1 = b.im / b.re; t2 = b.re + b.im * t1; temp = (a.re + a.im * t1) / t2; a.im = (a.im - a.re * t1) / t2; a.re = temp; } else { t1 = b.re / b.im; t2 = b.re * t1 + b.im; temp = (a.re * t1 + a.im) / t2; a.im = (a.im * t1 - a.re) / t2; a.re = temp; } return a; } // scalar /= complex template inline scalar& operator /= (float &a, const complex& b) { assert(fabs(b.im) < small_epsilon); a /= b.re; return a; } template inline scalar& operator /= (double &a, const complex& b) { assert(fabs(b.im) < small_epsilon); a /= b.re; return a; } template inline scalar& operator /= (long double &a, const complex& b) { assert(fabs(b.im) < small_epsilon); a /= b.re; return a; } // complex /= scalar template inline complex& operator /= (complex& a, const int &b) { a.re /= b; a.im /= b; return a; } template inline complex& operator /= (complex& a, const long int &b) { a.re /= b; a.im /= b; return a; } template inline complex& operator /= (complex& a, const float &b) { a.re /= b; a.im /= b; return a; } template inline complex& operator /= (complex& a, const double &b) { a.re /= b; a.im /= b; return a; } template inline complex& operator /= (complex& a, const long double &b) { a.re /= b; a.im /= b; return a; } //-------------------------------------------------------------- // unary operators ! + - //-------------------------------------------------------------- template inline complex operator! (const complex& a) //complex conj. {return complex(a.re,-a.im);} template inline complex operator+ (const complex& a) //complex conj. {return complex(a.re,a.im);} template inline complex operator- (const complex& a) //complex conj. {return complex(-a.re,-a.im);} //--------------------------------------------------------------- // binary operators + - * / template inline complex operator+ (const complex &a, const complex &b) {return complex(a.re+b.re, a.im+b.im);} //.. template inline complex operator+ (const complex &a, const int &b) {return complex(a.re+b, a.im);} template inline complex operator+ (const complex &a, const long int &b) {return complex(a.re+b, a.im);} template inline complex operator+ (const complex &a, const float &b) {return complex(a.re+b, a.im);} template inline complex operator+ (const complex &a, const double &b) {return complex(a.re+b, a.im);} template inline complex operator+ (const complex &a, const long double &b) {return complex(a.re+b, a.im);} //.. template inline complex operator+ (const int &a, const complex &b) {return complex(a+b.re, b.im);} template inline complex operator+ (const long int &a, const complex &b) {return complex(a+b.re, b.im);} template inline complex operator+ (const float &a, const complex &b) {return complex(a+b.re, b.im);} template inline complex operator+ (const double &a, const complex &b) {return complex(a+b.re, b.im);} template inline complex operator+ (const long double &a, const complex &b) {return complex(a+b.re, b.im);} //--------- template inline complex operator- (const complex &a, const complex &b) {return complex(a.re-b.re, a.im-b.im);} //.. template inline complex operator- (const complex &a, const int &b) {return complex(a.re-b, a.im);} template inline complex operator- (const complex &a, const long int &b) {return complex(a.re-b, a.im);} template inline complex operator- (const complex &a, const float &b) {return complex(a.re-b, a.im);} template inline complex operator- (const complex &a, const double &b) {return complex(a.re-b, a.im);} template inline complex operator- (const complex &a, const long double &b) {return complex(a.re-b, a.im);} //.. template inline complex operator- (const int &a, const complex &b) {return complex(a-b.re, -b.im);} template inline complex operator- (const long int &a, const complex &b) {return complex(a-b.re, -b.im);} template inline complex operator- (const float &a, const complex &b) {return complex(a-b.re, -b.im);} template inline complex operator- (const double &a, const complex &b) {return complex(a-b.re, -b.im);} template inline complex operator- (const long double &a, const complex &b) {return complex(a-b.re, -b.im);} //--------- template inline complex operator* (const complex& a, const complex& b) { scalar t1 = a.re * b.re; scalar t2 = a.im * b.im; return complex(t1-t2, (a.re+a.im)* (b.re+b.im)- (t1+t2) ); // why this trick? imag = a.r*b.i + b.r*a.i } //... template inline complex operator* (const complex& a, const int &b) {return complex(a.re*b, a.im*b);} template inline complex operator* (const complex& a, const long int &b) {return complex(a.re*b, a.im*b);} template inline complex operator* (const complex& a, const float &b) {return complex(a.re*b, a.im*b);} template inline complex operator* (const complex& a, const double &b) {return complex(a.re*b, a.im*b);} template inline complex operator* (const complex& a, const long double &b) {return complex(a.re*b, a.im*b);} //... template inline complex operator* (const int &a, const complex& b) {return complex(a*b.re, a*b.im);} template inline complex operator* (const long int &a, const complex& b) {return complex(a*b.re, a*b.im);} template inline complex operator* (const float &a, const complex& b) {return complex(a*b.re, a*b.im);} template inline complex operator* (const double &a, const complex& b) {return complex(a*b.re, a*b.im);} template inline complex operator* (const long double &a, const complex& b) {return complex(a*b.re, a*b.im);} //------------ template inline complex operator/ (const complex& a, const complex& b) { scalar t1, t2; if (fabs(b.re) >= fabs(b.im)) { t1 = b.im / b.re; t2 = b.re + b.im * t1; return complex((a.re+a.im*t1)/t2, (a.im-a.re*t1)/t2); } else { t1 = b.re / b.im; t2 = b.re*t1 + b.im; return complex((a.re*t1 + a.im) /t2, (a.im*t1 - a.re) /t2); } } //.. template inline complex operator/ (const int &a, const complex& b) {return complex(a,0)/b;} template inline complex operator/ (const long int &a, const complex& b) {return complex(a,0)/b;} template inline complex operator/ (const float &a, const complex& b) {return complex(a,0)/b;} template inline complex operator/ (const double &a, const complex& b) {return complex(a,0)/b;} template inline complex operator/ (const long double &a, const complex& b) {return complex(a,0)/b;} //.. template inline complex operator/ (const complex& a, const int &b) {return complex(a.re/b,a.im/b);} template inline complex operator/ (const complex& a, const long int &b) {return complex(a.re/b,a.im/b);} template inline complex operator/ (const complex& a, const float &b) {return complex(a.re/b,a.im/b);} template inline complex operator/ (const complex& a, const double &b) {return complex(a.re/b,a.im/b);} template inline complex operator/ (const complex& a, const long double &b) {return complex(a.re/b,a.im/b);} //------------------------------------------------------------------ // stream operators //------------------------------------------------------------------ template std::ostream& operator << (std::ostream &stream, const complex &a) { stream<<"("< std::istream& operator >> (std::istream &stream, complex &a) { stream>>a.re>>a.im; return stream; } //--------------------------------------------------------------------- // bool operators == != //--------------------------------------------------------------------- template inline bool operator == (const complex& a, const complex& b) { return (a.re==b.re) && (a.im==b.im); } //... template inline bool operator == (const int &a, const complex& b) { return (a==b.re) && (b.im inline bool operator == (const long int &a, const complex& b) { return (a==b.re) && (b.im inline bool operator == (const float &a, const complex& b) { return (a==b.re) && (b.im inline bool operator == (const double &a, const complex& b) { return (a==b.re) && (b.im inline bool operator == (const long double &a, const complex& b) { return (a==b.re) && (b.im inline bool operator == (const complex& a, const int &b) { return (a.re==b) && (a.im inline bool operator == (const complex& a, const long int &b) { return (a.re==b) && (a.im inline bool operator == (const complex& a, const float &b) { return (a.re==b) && (a.im inline bool operator == (const complex& a, const double &b) { return (a.re==b) && (a.im inline bool operator == (const complex& a, const long double &b) { return (a.re==b) && (a.im inline bool operator != (const complex& a, const complex& b) { return (a.re!=b.re) || (a.im!=b.im); } //... template inline bool operator != (const int &a, const complex& b) { return (a!=b.re) || (b.im>=small_epsilon); } template inline bool operator != (const long int &a, const complex& b) { return (a!=b.re) || (b.im>=small_epsilon); } template inline bool operator != (const float &a, const complex& b) { return (a!=b.re) || (b.im>=small_epsilon); } template inline bool operator != (const double &a, const complex& b) { return (a!=b.re) || (b.im>=small_epsilon); } template inline bool operator != (const long double &a, const complex& b) { return (a!=b.re) || (b.im>=small_epsilon); } //... template inline bool operator != (const complex& a, const int &b) { return (a.re!=b) || (a.im>=small_epsilon); } template inline bool operator != (const complex& a, const long int &b) { return (a.re!=b) || (a.im>=small_epsilon); } template inline bool operator != (const complex& a, const float &b) { return (a.re!=b) || (a.im>=small_epsilon); } template inline bool operator != (const complex& a, const double &b) { return (a.re!=b) || (a.im>=small_epsilon); } template inline bool operator != (const complex& a, const long double &b) { return (a.re!=b) || (a.im>=small_epsilon); } //--------------------------------------------------------------------- // bool operators >= <= //--------------------------------------------------------------------- template inline bool operator >= (const complex& a, const complex& b) { return (abs(a) >= abs(b)); } //... template inline bool operator >= (const int &a, const complex& b) { return (a >= abs(b)); } template inline bool operator >= (const long int &a, const complex& b) { return (a >= abs(b)); } template inline bool operator >= (const float &a, const complex& b) { return (a >= abs(b)); } template inline bool operator >= (const double &a, const complex& b) { return (a >= abs(b)); } template inline bool operator >= (const long double &a, const complex& b) { return (a >= abs(b)); } //... template inline bool operator >= (const complex& a, const int &b) { return (abs(a) >= b); } template inline bool operator >= (const complex& a, const long int &b) { return (abs(a) >= b); } template inline bool operator >= (const complex& a, const float &b) { return (abs(a) >= b); } template inline bool operator >= (const complex& a, const double &b) { return (abs(a) >= b); } template inline bool operator >= (const complex& a, const long double &b) { return (abs(a) >= b); } //... template inline bool operator <= (const complex& a, const complex