#ifndef __MY_DCOMPLEX__ #define __MY_DCOMPLEX__ #include #include #include class dcomplex{ private: double re, im; public: dcomplex (double r = 0, double i = 0): re (r), im (i) { } dcomplex& operator+= (const double r); dcomplex& operator-= (const double r); dcomplex& operator*= (const double r); dcomplex& operator/= (const double r); dcomplex& operator+= (const dcomplex& r); dcomplex& operator-= (const dcomplex& r); dcomplex& operator*= (const dcomplex& r); dcomplex& operator/= (const dcomplex& r); double& real() { return re; } double& imag() { return im; } double real() const { return re; } double imag() const { return im; } // dcomplex& conj() { im = -im; return *this;} dcomplex conj() const {return dcomplex(re,-im);} void Set(double r, double i){re=r; im=i;} void Add(double r, double i){re+=r, im+=i;}; private: friend dcomplex operator + (const dcomplex& x, double y); friend dcomplex operator + (double x, const dcomplex& y); friend dcomplex operator + (const dcomplex& x, const dcomplex& y); friend dcomplex operator - (const dcomplex& x, double y); friend dcomplex operator - (double x, const dcomplex& y); friend dcomplex operator - (const dcomplex& x, const dcomplex& y); friend dcomplex operator * (const dcomplex& x, double y); friend dcomplex operator * (double x, const dcomplex& y); friend dcomplex operator * (const dcomplex& x, const dcomplex& y); friend dcomplex operator / (const dcomplex& x, double y); friend dcomplex operator / (double x, const dcomplex& y); friend dcomplex operator / (const dcomplex& x, const dcomplex& y); friend bool operator == (const dcomplex& x, double y); friend bool operator == (double x, const dcomplex& y); friend bool operator == (const dcomplex& x, const dcomplex& y); friend bool operator != (const dcomplex& x, double y); friend bool operator != (double x, const dcomplex& y); friend bool operator != (const dcomplex& x, const dcomplex& y); friend dcomplex conj (const dcomplex& r); friend std::ostream& operator<< (std::ostream& stream, const dcomplex& r); friend std::istream& operator>> (std::istream& stream, dcomplex& r); friend double real(const dcomplex& r) { return r.re;} friend double imag(const dcomplex& r) { return r.im;} friend double& real(dcomplex& r) { return r.re;} friend double& imag(dcomplex& r) { return r.im;} friend dcomplex sqrt(const dcomplex& z); friend dcomplex sqrt_(const dcomplex& z); friend dcomplex sqrt(const dcomplex& z1, const dcomplex x2); friend bool isnan(const dcomplex& z) { return ::isnan(z.re) || ::isnan(z.im);} friend double arg(const dcomplex& x) { return atan2(x.imag(), x.real());} friend double norm (const dcomplex& r) { return r.re*r.re+r.im*r.im;} friend double abs(const dcomplex& r) { return sqrt(norm(r));} friend dcomplex log(const dcomplex& x); friend dcomplex operator- (const dcomplex& r); friend double real_product(const dcomplex& z1, const dcomplex& z2) {return z1.re*z2.re-z1.im*z2.im;} friend double imag_product(const dcomplex& z1, const dcomplex& z2) {return z1.re*z2.im+z1.im*z2.re;} friend dcomplex coth(const dcomplex& z); friend dcomplex atan(const dcomplex& z); }; template inline T sqr(const T& x) { return x*x; } inline dcomplex& dcomplex::operator+= (const double r) { re += r; return *this; } inline dcomplex& dcomplex::operator-= (const double r) { re -= r; return *this; } inline dcomplex& dcomplex::operator*= (const double r) { re *= r; im *= r; return *this; } inline dcomplex& dcomplex::operator/= (const double r) { re /= r; im /= r; return *this; } inline dcomplex& dcomplex::operator+= (const dcomplex& r) { re += r.re; im += r.im; return *this; } inline dcomplex& dcomplex::operator-= (const dcomplex& r) { re -= r.re; im -= r.im; return *this; } inline dcomplex& dcomplex::operator*= (const dcomplex& r) { double a = re * r.re - im * r.im; im = re * r.im + im * r.re; re = a; return *this; } inline dcomplex& dcomplex::operator/= (const dcomplex& r) { double norm_ = norm(r); double a = re * r.re + im * r.im; im = im * r.re - re * r.im; re = a/norm_; im /= norm_; return *this; } inline dcomplex operator + (const dcomplex& x, double y) { return dcomplex(x.re + y, x.im); } inline dcomplex operator + (double x, const dcomplex& y) { return dcomplex(x + y.re, y.im); } inline dcomplex operator + (const dcomplex& x, const dcomplex& y) { return dcomplex(x.re + y.re, x.im + y.im); } inline dcomplex operator - (const dcomplex& x, double y) { return dcomplex(x.re - y, x.im); } inline dcomplex operator - (double x, const dcomplex& y) { return dcomplex(x - y.re, -y.im); } inline dcomplex operator - (const dcomplex& x, const dcomplex& y) { return dcomplex(x.re - y.re, x.im - y.im); } inline dcomplex operator * (const dcomplex& x, double y) { return dcomplex(x.re * y, x.im * y); } inline dcomplex operator * (double x, const dcomplex& y) { return dcomplex(x * y.re, x * y.im); } inline dcomplex operator * (const dcomplex& x, const dcomplex& y) { return dcomplex(x.re * y.re - x.im * y.im, x.re * y.im + x.im * y.re); } inline dcomplex operator / (const dcomplex& x, double y) { return dcomplex(x.re/y, x.im/y); } inline dcomplex operator / (double x, const dcomplex& y) { double xnorm = x/norm(y); return dcomplex(xnorm * y.re, -xnorm * y.im); } inline dcomplex operator / (const dcomplex& x, const dcomplex& y) { double norm_ = norm(y); return dcomplex((x.re * y.re + x.im * y.im)/norm_, (x.im * y.re - x.re * y.im)/norm_); } inline bool operator == (const dcomplex& x, double y) { return (x.re==y) && (x.im==0); } inline bool operator == (double x, const dcomplex& y) { return (x==y.re) && (y.im==0); } inline bool operator == (const dcomplex& x, const dcomplex& y) { return (x.re==y.re) && (x.im==y.im); } inline bool operator != (const dcomplex& x, double y) { return !(operator==(x,y)); } inline bool operator != (double x, const dcomplex& y) { return !(operator==(x,y)); } inline bool operator != (const dcomplex& x, const dcomplex& y) { return !(operator==(x,y)); } inline dcomplex operator- (const dcomplex& r) { return dcomplex(-r.re,-r.im); } inline dcomplex conj (const dcomplex& r) { return dcomplex(r.re, -r.im); } inline dcomplex log(const dcomplex& x) { return dcomplex(log(norm(x))*0.5, arg(x)); } inline dcomplex sqrt(const dcomplex& z1, const dcomplex z2) { double fi1 = atan2(z1.im, z1.re); double fi2 = atan2(z2.im, z2.re); double fi = 0.5*(fi1+fi2); double r = sqrt(sqrt((sqr(z1.re)+sqr(z1.im))*(sqr(z2.re)+sqr(z2.im)))); return dcomplex(r*cos(fi),r*sin(fi)); } inline dcomplex sqrt(const dcomplex& x) { double r = abs (x); double nr, ni; if (r == 0.0) nr = ni = r; else if (real (x) > 0) { nr = sqrt (0.5 * (r + real (x))); ni = imag (x) / nr / 2; } else { ni = sqrt (0.5 * (r - real (x))); if (imag (x) < 0) ni = - ni; nr = imag (x) / ni / 2; } return dcomplex (nr, ni); } inline dcomplex sqrt_(const dcomplex& x) { double r = sqrt(abs(x)); double fi = atan2(x.im, x.re); return dcomplex(r*cos(0.5*fi), r*sin(0.5*fi)); } inline dcomplex coth(const dcomplex& z) { double ex = exp(z.real()), emx = exp(-z.real()); double cy = cos(z.imag()), sy = sin(z.imag()); double u1 = ex+emx, u2 = ex-emx; double im = 1./(sqr(u2*cy)+sqr(u1*sy)); return dcomplex (u1*u2*im, cy*sy*(sqr(u2)-sqr(u1))*im); } inline dcomplex atan(const dcomplex& z) { double x = z.real(), y = z.imag(), x2 = x*x, y2=y*y; double den = 1/(x2+sqr(y+1.)); double re = 0.5*atan2(2*x*den,(1-x2-y2)*den); double im = -0.25*log(den*(x2+sqr(y-1.))); return dcomplex(re,im); } inline std::ostream& operator<< (std::ostream& stream, const dcomplex& r){ int width = stream.width(); stream << std::setw(width)<< r.re << " " << std::setw(width) << r.im << " "; return stream; } inline std::istream& operator>> (std::istream& stream, dcomplex& r){ stream >> r.re; stream >> r.im; return stream; } #endif