complex.h

00001 #include <iostream>
00002 #include <cmath>
00003 #include "random.h"
00004 #include "array.h"
00005 class Complex
00006 {
00007 
00008 public:
00009 
00010         Complex() : _real(0), _imag(0) {}
00011         explicit Complex( double r) : _real(r), _imag(0) {}
00012         Complex(double r, double i) : _real(r), _imag(i) {}
00013 
00014         Complex& operator+=(const double& d)
00015         {
00016                 _real += d;
00017                 return *this;
00018         }
00019 
00020         Complex& operator+=(const Complex& c)
00021         {
00022                 _real += c._real;
00023                 _imag += c._imag;
00024                 return *this;
00025         }
00026 
00027         Complex& operator-=(const double &d)
00028         {
00029                 _real -= d;
00030                 return *this;
00031         }
00032 
00033         Complex& operator-=(const Complex& c)
00034         {
00035                 _real -= c._real;
00036                 _imag -= c._imag;
00037                 return *this;
00038         }
00039 
00040         Complex& operator*=(const double& d)
00041         {
00042                 _real *= d;
00043                 _imag *= d;
00044                 return *this;
00045         }
00046 
00047         Complex& operator*=(const Complex& c)
00048         {
00049                 double re = _real;
00050                 double im = _imag;
00051                 _real = re * c._real - im * c._imag;
00052                 _imag = re * c._imag + im * c._real;
00053                 return *this;
00054         }
00055 
00056         Complex& operator/=(const double& d)
00057         {
00058                 _real /= d;
00059                 _imag /= d;
00060                 return *this;
00061         }
00062 
00063         Complex& operator/=(const Complex& c)
00064         {
00065                 double re = _real;
00066                 double im = _imag;
00067                 double d = c._real * c._real + c._imag * c._imag;
00068                 _real = (re * c._real + im * c._imag) / d;
00069                 _imag = (im * c._real - re * c._imag) / d;
00070                 return *this;
00071         }
00072 
00073         Complex Conj() const
00074         {
00075                 return Complex(_real, -_imag);
00076         }
00077 
00078         Complex TransToRadiusAngle() const
00079         {
00080                 if (_real == 0 && _imag == 0)
00081                 {
00082                         return Complex(0,0);
00083                 }
00084 
00085                 double r = sqrt(_real*_real + _imag*_imag);
00086                 double sinx = _imag / r;
00087                 double alpha = asin(sinx); // alpha will be in -pi/2 to pi/2
00088                 if (alpha >=0) // means in 1st or 2nd quardrant
00089                 {
00090                         if (_real < 0)
00091                         {
00092                                 alpha = PI - alpha;
00093                         }
00094                 }
00095                 else 
00096                 {
00097                         if (_real < 0)
00098                         {
00099                                 alpha = -alpha + PI;
00100                         }
00101                         else
00102                         {
00103                                 alpha = 2*PI + alpha;
00104                         }
00105                 }
00106 
00107                 return Complex(r, alpha);               
00108         }
00109 
00110         Array<Complex> Sqrt() const
00111         {
00112                 //
00113                 Array<Complex> roots;
00114 
00115                 if (_real == 0 && _imag == 0)
00116                 {                                 
00117                         roots.append(Complex(0,0));
00118                         return roots;
00119                 }
00120 
00121                 // transform to radius and angle form 
00122 
00123                 Complex ra = TransToRadiusAngle();
00124                 double rRoot = sqrt(ra.Real());
00125                 
00126 
00127                 for(int i = 0; i < 2; i++)
00128                 {
00129                         double angle = (ra.Imag() + 2*double(i)*PI)/2;
00130                         double real = rRoot * cos(angle);
00131                         double imag = rRoot * sin(angle);
00132                         roots.append(Complex(real,imag)); 
00133                 }
00134 
00135                 return roots;                           
00136         }
00137 
00138         Array<Complex> Root(int k) const
00139         {
00140                 Array<Complex> roots;
00141                 if (_real == 0 && _imag == 0)
00142                 {
00143                         roots.append(Complex(0,0));
00144                         return roots;
00145                 }
00146                 Complex ra = TransToRadiusAngle();
00147                 double rRoot = sqrt(ra.Real());
00148                 
00149 
00150                 for(int i = 0; i < k; i++)
00151                 {
00152                         double angle = (ra.Imag() + 2*double(i)*PI)/double(k);
00153                         double real = rRoot * cos(angle);
00154                         double imag = rRoot * sin(angle);
00155                         roots.append(Complex(real,imag)); 
00156                 }
00157 
00158                 return roots;                           
00159         }
00160         double Real() const { return _real; }
00161         double Imag() const { return _imag; }
00162         void Real(const double& re) { _real = re ; }
00163         void Imag(const double& im) { _imag = im ; }
00164         void Set(const double& re, const double& im){ _real = re; _imag = im ; }
00165         double Modsq() const { return _real*_real + _imag * _imag ; }
00166         double Mod() const { return sqrt(_real*_real + _imag * _imag); }
00167 
00168 private:
00169         double _real;
00170         double _imag;
00171 
00172 };
00173 
00174 inline Complex operator+(const Complex& c)
00175 {
00176         return Complex(c.Real(), c.Imag());
00177 }
00178 
00179 inline Complex operator-(const Complex& c)
00180 {
00181         return Complex(-c.Real(), -c.Imag());
00182 }
00183 
00184 inline Complex operator+(const Complex& c, const double& d)
00185 {
00186         return Complex(c.Real() + d, c.Imag());
00187 }
00188 
00189 inline Complex operator+(const double& d, const Complex& c)
00190 {
00191         return Complex(d + c.Real(), c.Imag());
00192 }
00193 
00194 inline Complex operator+(const Complex& c1, const Complex& c2)
00195 {
00196         return Complex(c1.Real() + c2.Real(), c1.Imag() + c2.Imag());
00197 }
00198 
00199 inline Complex operator-(const Complex& c, const double& d)
00200 {
00201         return Complex(c.Real() - d, c.Imag());
00202 }
00203 
00204 inline Complex operator-(const double& d, const Complex& c)
00205 {
00206         return Complex(d - c.Real(), -c.Imag());
00207 }
00208 
00209 inline Complex operator-(const Complex& c1, const Complex& c2)
00210 {
00211         return Complex(c1.Real() - c2.Real(), c1.Imag() - c2.Imag());
00212 }
00213 
00214 
00215 inline Complex operator*(const Complex& c, const double& d)
00216 {
00217         return Complex(c.Real() * d, c.Imag() * d);
00218 }
00219 
00220 inline Complex operator*(const double& d, const Complex& c)
00221 {
00222         return Complex(c.Real() * d, c.Imag() * d);
00223 }
00224 
00225 inline Complex operator*(const Complex& c1, const Complex& c2)
00226 {
00227         double real = c1.Real() * c2.Real() - c1.Imag() * c2.Imag();
00228         double imag = c1.Real() * c2.Imag() + c1.Imag() * c2.Real();
00229         return Complex(real, imag);
00230 }
00231 
00232 inline Complex operator/(const Complex& c, const double& d)
00233 {
00234         return Complex(c.Real() / d, c.Imag() / d);
00235 }
00236 
00237 inline Complex operator/(const double& d, const Complex& c)
00238 {
00239         double dd = c.Real() * c.Real() + c.Imag() * c.Imag();
00240         return Complex((d * c.Real())/dd, (-d * c.Imag())/dd);
00241 }
00242 
00243 inline Complex operator/(const Complex& c1, const Complex& c2)
00244 {
00245         double d = c2.Real() * c2.Real() + c2.Imag() * c2.Imag();
00246         double real = (c1.Real() * c2.Real() + c1.Imag() * c2.Imag()) / d;
00247         double imag = (c1.Imag() * c2.Real() - c1.Real() * c2.Imag()) / d;
00248         return Complex(real, imag);
00249 }
00250 
00251 inline double real(const Complex &c)
00252 {
00253         return c.Real();
00254 }
00255 
00256 inline double imag(const Complex &c)
00257 {
00258         return c.Imag();
00259 }
00260 
00261 inline double abs(const Complex &c)
00262 {
00263         return sqrt(c.Real() * c.Real() + c.Imag() * c.Imag());
00264 }
00265 
00266 inline double norm(const Complex &c)
00267 {
00268         return c.Real() * c.Real() + c.Imag() * c.Imag();
00269 }
00270 
00271 inline Complex conj(const Complex &c)
00272 {
00273         return Complex(c.Real(), -c.Imag());
00274 }
00275 
00276 inline ostream &operator<<(ostream &os, const Complex &c)
00277 {
00278         os<<c.Real()<<"+"<<c.Imag()<<"i";
00279         return os;
00280 } 

Generated on Sun Jun 7 11:55:21 2009 for Alchemy by  doxygen 1.5.1