#ifndef COMPLEX_H
#define COMPLEX_H 1

#include <cmath>
#include <iostream>

class Complex {
private:
    double re;
    double im;
public:
    Complex(double r = 0., double i = 0.):
        re(r),
        im(i)
    {}
    
    double real() const { return re; }
    double imag() const { return im; }
    void real(double r) { re = r; }
    void imag(double i) { im = i; }
    
    Complex& operator+=(const Complex& z) {
        re += z.re;
        im += z.im;
        return *this;
    }
    
    Complex& operator-=(const Complex& z) {
        re -= z.re;
        im -= z.im;
        return *this;
    }
    
    Complex& operator*=(const Complex& z) {
        double r = re*z.re - im*z.im;
        double i = re*z.im + im*z.re;
        re = r;
        im = i;
        return *this;
    }
    
    Complex inverse() const {
        double abs2 = re*re + im*im;
        return Complex(
            re/abs2, -im/abs2
        );
    }
    
    Complex& operator/=(const Complex& z) {
        (*this) *= z.inverse();
        return *this;
    }    
};

inline Complex operator+(
    const Complex& z1, const Complex& z2
) {
    Complex res = z1;
    res += z2;
    return res;
}

inline Complex operator-(
    const Complex& z1, const Complex& z2
) {
    Complex res = z1;
    res -= z2;
    return res;
}

inline Complex operator*(
    const Complex& z1, const Complex& z2
) {
    Complex res = z1;
    res *= z2;
    return res;
}

inline Complex operator/(
    const Complex& z1, const Complex& z2
) {
    Complex res = z1;
    res /= z2;
    return res;
}

inline Complex conj(const Complex& z) {
    return Complex(z.real(), -z.imag());
}

inline double real(const Complex& z) { return z.real(); }
inline double imag(const Complex& z) { return z.imag(); }

inline double arg(const Complex& z) {
    return atan2(z.imag(), z.real());
}

inline double abs(const Complex& z) {
    return sqrt(
        z.real()*z.real() + z.imag()*z.imag()
    ); 
}

inline Complex polar(double magnitude, double angle) {
    return Complex(
        magnitude*cos(angle),
        magnitude*sin(angle)
    );
}

Complex rootN(const Complex& z, int n);
Complex unitRootN(int n);

inline std::ostream& operator<<(
    std::ostream& s, const Complex& z
) {
    s << z.real() << " + " << z.imag() << "*i";
    return s;
} 

inline std::istream& operator>>(
    std::istream& s, Complex& z
) {
    double r, i;
    if (s >> r >> i) {
        z.real(r); z.imag(i);
    }    
    return s;
} 

#endif