00001 /*************************************************************************** 00002 ** 00003 ** This file is part of QGpCoreTools. 00004 ** 00005 ** This library is free software; you can redistribute it and/or 00006 ** modify it under the terms of the GNU Lesser General Public 00007 ** License as published by the Free Software Foundation; either 00008 ** version 2.1 of the License, or (at your option) any later version. 00009 ** 00010 ** This file is distributed in the hope that it will be useful, but WITHOUT 00011 ** ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 00012 ** FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 00013 ** License for more details. 00014 ** 00015 ** You should have received a copy of the GNU Lesser General Public 00016 ** License along with this library; if not, write to the Free Software 00017 ** Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 00018 ** 00019 ** See http://www.geopsy.org for more information. 00020 ** 00021 ** Created : 2004-01-21 00022 ** Authors : 00023 ** Marc Wathelet 00024 ** Marc Wathelet (ULg, Liège, Belgium) 00025 ** Marc Wathelet (LGIT, Grenoble, France) 00026 ** 00027 ***************************************************************************/ 00028 00029 #ifndef COMPLEX_H 00030 #define COMPLEX_H 00031 00032 #include <math.h> 00033 #include <QtCore> 00034 00035 #include "QGpCoreToolsDLLExport.h" 00036 00037 namespace QGpCoreTools { 00038 00039 class StringSection; 00040 00041 // Trying to use SSE extension 00042 typedef double v2df __attribute__ ((vector_size (16))); 00043 00044 class QGPCORETOOLS_EXPORT Complex 00045 { 00046 public: 00047 inline Complex(const Complex& c); 00048 inline Complex(double re=0, double im=0); 00049 ~Complex() {} 00050 00051 inline bool operator==(const Complex& c) const; 00052 00053 inline Complex& operator=(const Complex& c); 00054 inline Complex& operator+=(const Complex& c); 00055 inline Complex& operator-=(const Complex& c); 00056 inline Complex& operator*=(const Complex& c); 00057 inline Complex& operator/=(const Complex& c); 00058 00059 inline Complex& operator=(const double& d); 00060 inline Complex& operator+=(const double& d); 00061 inline Complex& operator-=(const double& d); 00062 inline Complex& operator*=(const double& d); 00063 inline Complex& operator/=(const double& d); 00064 00065 inline Complex operator+(const Complex& c) const; 00066 inline Complex operator-(const Complex& c) const; 00067 inline Complex operator*(const Complex& c) const; 00068 inline Complex operator/(const Complex& c) const; 00069 00070 inline Complex operator+(const double& d) const; 00071 inline Complex operator-(const double& d) const; 00072 inline Complex operator*(const double& d) const; 00073 inline Complex operator/(const double& d) const; 00074 00075 bool isNull() const {return _data.values.re==0 && _data.values.im==0;} 00076 double re() const {return _data.values.re;} 00077 double im() const {return _data.values.im;} 00078 void set(double re, double im) {_data.values.re=re; _data.values.im=im;} 00079 void setRe(double re) {_data.values.re=re;} 00080 void setIm(double im) {_data.values.im=im;} 00081 void setExp(double fac, double arg) 00082 { 00083 _data.values.re=fac*::cos(arg); 00084 _data.values.im=fac*::sin(arg); 00085 } 00086 void setUnitExp(double arg) 00087 { 00088 _data.values.re=::cos(arg); 00089 _data.values.im=::sin(arg); 00090 } 00091 double abs() const {return ::sqrt(abs2());} 00092 double abs2() const {return _data.values.re * _data.values.re + _data.values.im * _data.values.im;} 00093 double phase() const {return atan2(_data.values.im, _data.values.re);} 00094 void setAbs(double fac) {setExp(fac, phase());} 00095 void setPhase(double arg) {setExp(abs(), arg);} 00096 00097 bool fromString(const StringSection& str); 00098 QString toString(int precision=6) const; 00099 00100 static const Complex null; 00101 private: 00102 union { 00103 v2df vector; 00104 struct { 00105 double re, im; 00106 } values; 00107 } _data; 00108 }; 00109 00110 inline Complex operator+(const double& d, const Complex& c); 00111 inline Complex operator-(const double& d, const Complex& c); 00112 inline Complex operator*(const double& d, const Complex& c); 00113 inline Complex operator/(const double& d, const Complex& c); 00114 00115 inline Complex sqrt(const Complex& c); 00116 inline Complex cos(const Complex& c); 00117 inline Complex sin(const Complex& c); 00118 inline Complex tan(const Complex& c); 00119 inline Complex asin(const Complex& c); 00120 inline Complex log(const Complex& c); 00121 inline Complex exp(const Complex& c); 00122 inline Complex conjugate(const Complex& c); 00123 inline Complex inverse(const Complex& c); 00124 inline double abs(const Complex& c); 00125 inline double abs2(const Complex& c); 00126 00127 QGPCORETOOLS_EXPORT QTextStream& operator<< (QTextStream& s, const Complex& c); 00128 QGPCORETOOLS_EXPORT QTextStream& operator>> (QTextStream& s, Complex& c); 00129 QGPCORETOOLS_EXPORT QDataStream& operator<< (QDataStream& s, const Complex& c); 00130 QGPCORETOOLS_EXPORT QDataStream& operator>> (QDataStream& s, Complex& c); 00131 00132 inline Complex::Complex(const Complex& c) 00133 { 00134 #ifdef QT_HAVE_SSE 00135 _data.vector=c._data.vector; 00136 #else 00137 _data.values.re=c._data.values.re; 00138 _data.values.im=c._data.values.im; 00139 #endif 00140 } 00141 00142 inline Complex::Complex(double re, double im) 00143 { 00144 _data.values.re=re; 00145 _data.values.im=im; 00146 } 00147 00148 inline bool Complex::operator==(const Complex& c) const 00149 { 00150 return _data.values.re==c._data.values.re && _data.values.im==c._data.values.im; 00151 } 00152 00153 inline Complex& Complex::operator=(const double& d) 00154 { 00155 _data.values.re=d; 00156 _data.values.im=0; 00157 return *this; 00158 } 00159 00160 inline Complex& Complex::operator=(const Complex& c) 00161 { 00162 #ifdef QT_HAVE_SSE 00163 _data.vector=c._data.vector; 00164 #else 00165 _data.values.re=c._data.values.re; 00166 _data.values.im=c._data.values.im; 00167 #endif 00168 return *this; 00169 } 00170 00171 inline Complex& Complex::operator+=(const double& d) 00172 { 00173 _data.values.re += d; 00174 return *this; 00175 } 00176 00177 inline Complex& Complex::operator+=(const Complex& c) 00178 { 00179 #ifdef QT_HAVE_SSE 00180 _data.vector+=c._data.vector; 00181 #else 00182 _data.values.re+=c._data.values.re; 00183 _data.values.im+=c._data.values.im; 00184 #endif 00185 return *this; 00186 } 00187 00188 inline Complex& Complex::operator-=(const double& d) 00189 { 00190 _data.values.re -= d; 00191 return *this; 00192 } 00193 00194 inline Complex& Complex::operator-=(const Complex& c) 00195 { 00196 #ifdef QT_HAVE_SSE 00197 _data.vector -= c._data.vector; 00198 #else 00199 _data.values.re-=c._data.values.re; 00200 _data.values.im-=c._data.values.im; 00201 #endif 00202 return *this; 00203 } 00204 00205 inline Complex& Complex::operator*=(const double& d) 00206 { 00207 #if defined(QT_HAVE_3DNOW) && defined(QT_HAVE_SSE) 00208 v2df factor=__builtin_ia32_movddup(d); 00209 __builtin_ia32_mulpd(_data.vector, factor); 00210 #else 00211 _data.values.re *= d; 00212 _data.values.im *= d; 00213 #endif 00214 return *this; 00215 } 00216 00217 inline Complex& Complex::operator*=(const Complex& c) 00218 { 00219 #if defined(QT_HAVE_3DNOW) && defined(QT_HAVE_SSE) 00220 v2df thisIm=__builtin_ia32_shufpd(_data.vector, _data.vector,3); // Im part of this in both 00221 v2df thisRe=__builtin_ia32_shufpd(_data.vector, _data.vector,0); // Re part of this in both 00222 // Next two lines not accepted by gcc ?? even with -msse3 option 00223 //v2df thisIm=__builtin_ia32_movddup(_data.values.im); // Im part of this in both 00224 //v2df thisRe=__builtin_ia32_movddup(_data.values.re); // Re part of this in both 00225 v2df c_flipped=__builtin_ia32_shufpd(c._data.vector,c._data.vector,1); // Swap re and im parts of c 00226 _data.vector=__builtin_ia32_addsubpd(__builtin_ia32_mulpd(thisRe, c._data.vector), 00227 __builtin_ia32_mulpd(thisIm, c_flipped)); 00228 #elif defined(QT_HAVE_SSE2) 00229 v2df thisIm=__builtin_ia32_shufpd(_data.vector, _data.vector,3); // Im part of this in both 00230 v2df thisRe=__builtin_ia32_shufpd(_data.vector, _data.vector,0); // Re part of this in both 00231 v2df c_flipped=__builtin_ia32_shufpd(c._data.vector,c._data.vector,1); // Swap re and im parts of c 00232 static const union { // (signbit,0) 00233 int i[4]; v2df v; 00234 } signbitlow={{0,0x80000000,0,0}}; 00235 thisIm=__builtin_ia32_xorpd(thisIm, signbitlow.v); // Change sign of low 00236 // Multiply and add: 00237 _data.vector=__builtin_ia32_mulpd(thisRe, c._data.vector) + 00238 __builtin_ia32_mulpd(thisIm, c_flipped); 00239 #else 00240 double tmpre=_data.values.re * c._data.values.re - _data.values.im * c._data.values.im; 00241 _data.values.im=_data.values.re * c._data.values.im + _data.values.im * c._data.values.re; 00242 _data.values.re=tmpre; 00243 #endif 00244 return *this; 00245 } 00246 00247 inline Complex& Complex::operator/=(const double& d) 00248 { 00249 double invd=1.0/d; 00250 _data.values.re*=invd; 00251 _data.values.im*=invd; 00252 return *this; 00253 } 00254 00255 inline Complex& Complex::operator/=(const Complex& c) 00256 { 00257 operator*=(conjugate(c)); 00258 operator/=(c.abs2()); 00259 return *this; 00260 } 00261 00262 inline Complex Complex::operator+(const double& d) const 00263 { 00264 return Complex(_data.values.re+d, _data.values.im); 00265 } 00266 00267 inline Complex operator+(const double& d, const Complex& c) 00268 { 00269 return c+d; 00270 } 00271 00272 inline Complex Complex::operator+(const Complex& c) const 00273 { 00274 Complex r; 00275 #ifdef QT_HAVE_SSE 00276 r._data.vector=_data.vector+c._data.vector; 00277 #else 00278 r._data.values.re=_data.values.re+c._data.values.re; 00279 r._data.values.im=_data.values.im+c._data.values.im; 00280 #endif 00281 return r; 00282 } 00283 00284 inline Complex Complex::operator-(const double& d) const 00285 { 00286 return Complex(_data.values.re-d, _data.values.im); 00287 } 00288 00289 inline Complex operator-(const double& d, const Complex& c) 00290 { 00291 return Complex(d-c.re(), -c.im()); 00292 } 00293 00294 inline Complex Complex::operator-(const Complex& c) const 00295 { 00296 Complex r; 00297 #ifdef QT_HAVE_SSE 00298 r._data.vector=_data.vector-c._data.vector; 00299 #else 00300 r._data.values.re=_data.values.re-c._data.values.re; 00301 r._data.values.im=_data.values.im-c._data.values.im; 00302 #endif 00303 return r; 00304 } 00305 00306 inline Complex Complex::operator*(const double& d) const 00307 { 00308 return Complex(_data.values.re*d, _data.values.im*d); 00309 } 00310 00311 inline Complex operator*(const double& d, const Complex& c) 00312 { 00313 return Complex(d*c.re(), d*c.im()); 00314 } 00315 00316 inline Complex Complex::operator*(const Complex& c) const 00317 { 00318 Complex r(*this); 00319 r*=c; 00320 return r; 00321 } 00322 00323 inline Complex Complex::operator/(const double& d) const 00324 { 00325 return Complex(_data.values.re/d, _data.values.im/d); 00326 } 00327 00328 inline Complex operator/(const double& d, const Complex& c) 00329 { 00330 double m2=c.abs2(); 00331 Complex r(d); 00332 r*=conjugate(c); 00333 r/=m2; 00334 return r; 00335 } 00336 00337 inline Complex Complex::operator/(const Complex& c) const 00338 { 00339 double m2=c.abs2(); 00340 Complex r(*this); 00341 r*=conjugate(c); 00342 r/=m2; 00343 return r; 00344 } 00345 00346 inline Complex sqrt(const Complex& c) 00347 { 00348 double r2=::sqrt(c.abs()); 00349 double phi2=c.phase()/2; 00350 return Complex(r2*::cos(phi2), r2*::sin(phi2)); 00351 } 00352 00353 inline Complex inverse(const Complex& c) 00354 { 00355 double m2=1.0/c.abs2(); 00356 return Complex(c.re() * m2, -m2 * c.im()); 00357 } 00358 00359 inline Complex exp(const Complex& c) 00360 { 00361 double r=::exp(c.re()); 00362 return Complex(r*::cos(c.im()), r*::sin(c.im())); 00363 } 00364 00365 inline Complex log(const Complex& c) 00366 { 00367 return Complex(::log(c.abs()), c.phase()); 00368 } 00369 00373 inline Complex asin(const Complex& c) 00374 { 00375 return Complex(0.0,-1.0)*log(sqrt(Complex(1.0)-c*c)+Complex(0.0,1.0)*c); 00376 } 00377 00378 inline Complex cos(const Complex& c) 00379 { 00380 return Complex( ::cos(c.re()) *::cosh(c.im()), -::sin(c.re()) * sinh(c.im())); 00381 } 00382 00383 inline Complex sin(const Complex& c) 00384 { 00385 return Complex( ::sin(c.re()) *::cosh(c.im()), ::cos(c.re()) * sinh(c.im())); 00386 } 00387 00388 inline Complex tan(const Complex& c) 00389 { 00390 return sin(c)/cos(c); 00391 } 00392 00393 inline Complex conjugate(const Complex& c) 00394 { 00395 return Complex(c.re(),-c.im()); 00396 } 00397 00398 double abs(const Complex& c) 00399 { 00400 return c.abs(); 00401 } 00402 00403 double abs2(const Complex& c) 00404 { 00405 return c.abs2(); 00406 } 00407 00408 } // namespace QGpCoreTools 00409 00410 #endif // COMPLEX_H
