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);
00088 if (alpha >=0)
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
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 }