// complex.cc
#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;
#define PI 3.14159265358979323846
#define E  2.71828182845904523536

template <class T>
T absval(T x) { return (x >= 0) ? x : -x;}

template <class REAL, class FLOAT>
class complex {
  REAL r; REAL i;
public:
  complex() {};
  complex(REAL re, REAL im = 0) { r=re; i=im;}
  // Addition:
  complex operator+=(const complex& y) {    // Member-Funktion
    r += y.r; i += y.i; return *this; }
  complex operator+=(REAL y) { r += y; return *this; }
  friend complex operator+(const complex& x, const complex& y) {
    return complex(x.r+y.r, x.i+y.i);   }
  friend complex operator+(const complex& x, REAL y) {
    return complex(x.r+y, x.i);   }
  friend complex operator+(REAL x, const complex& y) {
    return complex(x+y.r, y.i);   }
  // Subtraktion:
  complex operator-=(const complex& y) {    // Member-Funktion
    r -= y.r; i -= y.i; return *this; }
  complex operator-=(REAL y) { r -= y; return *this; }
  friend complex operator-(const complex& x, const complex& y) {
    return complex(x.r-y.r, x.i-y.i);   }
  friend complex operator-(const complex& x, REAL y) {
    return complex(x.r-y, x.i);   }
  friend complex operator-(REAL x, const complex& y) {
    return complex(x-y.r, -y.i);   }
  friend complex operator-(const complex& x) {
    return complex(-x.r, -x.i);   }
  // Multiplikation:
  complex operator*=(const complex& y) {
    REAL rr = r*y.r-i*y.i; i = r*y.i+i*y.r; r = rr; return *this; }
  complex operator*=(REAL y) { r *= y; i *= y; return *this; }
  friend complex operator*(const complex& x, const complex& y) {
    return complex(x.r*y.r-x.i*y.i, x.r*y.i+x.i*y.r);   }
  friend complex operator*(const complex& x, REAL y) {
    return complex(x.r*y, x.i*y);   }
  friend complex operator*(REAL x, const complex& y) {
    return complex(x*y.r, x*y.i);   }
  // Division:
  complex operator/=(const complex& y) {
    REAL a = y.r*y.r + y.i*y.i; REAL rr = (r*y.r+i*y.i)/a;
    i = (i*y.r-r*y.i)/a; r = rr; return *this; }
  complex operator/=(REAL y) { r /= y; i /= y; return *this; }
  friend complex operator/(const complex& x, const complex& y) {
    REAL a = y.r*y.r + y.i*y.i;
    return complex((x.r*y.r+x.i*y.i)/a, (x.i*y.r-x.r*y.i)/a);   }
  friend complex operator/(const complex& x, REAL y) {
    return complex(x.r/y, x.i/y);   }
  friend complex operator/(REAL x, const complex& y) {
    REAL a = y.r*y.r + y.i*y.i; x/=a;
    return complex(x*y.r, -x*y.i); }
  friend complex inv(complex z) {
    REAL a = z.r*z.r + z.i*z.i;
    return complex(z.r/a,-z.i/a);
  }
  // Absolutbetrag, Argument:
  friend FLOAT abs(const complex& x) {
    return sqrt(x.r*x.r+x.i*x.i); }
  friend FLOAT arg(const complex& x) {
    if (x.r == 0) return (x.i == 0) ? 0 : ( x.i > 0 ? 0.5*PI : 1.5*PI);
    if (x.r <  0) return atan(x.i/x.r) + PI;
    // x.r > 0 :
    return  (x.i >= 0) ? atan(x.i/x.r) : atan(x.i/x.r) + 2*PI ;
  }
  friend FLOAT warg(const complex& x) {  return arg(x)/PI*180; }
  // Wurzel:
  friend complex sqrt(const complex& x) {
    FLOAT a = arg(x)/2, b = sqrt(abs(x));
    return complex( REAL(cos(a)*b), REAL(sin(a)*b)); }
  // Potenzieren: Achtung: Nach der Prioritaetenliste fuer
  // Operatoren hat der Operator ^ nicht die Prioritaet, die man fuer das
  // Potenzieren erwartet, man muss also klammern: (u^v) * w
  friend complex log(complex z) {
    return complex(log(abs(z)),arg(z));
  }
  friend complex operator^(complex u, complex v) {
    complex z = v*log(u);  REAL  a = exp(z.r);
    return complex(a*cos(z.i), a*sin(z.i));
  }
  friend complex operator^(complex z, REAL d) { return z^complex(d);  }
  friend complex operator^(REAL u, complex z) { return complex(u)^z;  }
  // Boolsche Operatoren:
  friend bool operator==(const complex& x, const complex& y) {
    return (x.r==y.r && x.i==y.i); }
  friend bool operator==(const complex& x, REAL y) {
    return (x.i == 0 && x.r==y); }
  friend bool operator==(REAL x, const complex& y) {
    return (y.i == 0 && x==y.r); }
  friend bool operator!=(const complex& x, const complex& y) {
    return (x.r!=y.r || x.i!=y.i); }
  friend bool operator!=(const complex& x, REAL y) {
    return (x.r!=y || x.i!=0); }
  friend bool operator!=(REAL x, const complex& y) {
    return (x!=y.r || y.i!=0); }
  // Ausgabeoperator:
  friend ostream& operator<<(ostream& os, const complex z) {
    REAL rr, ii;
    rr = z.r; ii = z.i;
    if ( absval(rr) < 1.0e-15) rr=0;
    if ( absval(ii) < 1.0e-15) ii=0;
    if (ii == 0.0) {
      os<<setw(6)<<setprecision(2)<< rr;
      os<<"    ";  os<<setw(6)<<setprecision(2)<< ' ';
      return os;
    }
    if (rr == 0.0) {
      os<<setw(6)<<' '<<"   ";
      os<<setw(6)<<setprecision(2)<< ii << 'i';
      return os;
    }
    os<<setw(6)<<setprecision(2)<< rr;
    os << ( ii>=0 ? " + ":" - ");
    os<<setw(6)<<setprecision(2)<< ( ii>=0 ? ii : -ii ) << 'i';
    return os; }
}; // Jetzt hat man z. B. folgende Moeglichkeiten:
typedef complex<float,double>  complex_f;
typedef complex<double,double> complex_d;
typedef complex<long double,long double> complex_ld;
typedef complex<int,double>    complex_i;
typedef complex_ld comp;
int main() {
  comp u,v;
  v = u = sqrt(sqrt(comp(0, 1))); // 4-te Wurzel aus i
  cout << u << '\n';
  cout << "\n";
  for (int i=0; i<16; i++) {
    cout << u << " hoch " << setw(2) << i+1 << " : " << v << "   warg = ";
    cout << setw(6) << setprecision(4) << warg(v) << "°\n";
    v *= u;
  }
  cout << "\n";
  v = u = 1/u;
  for (int i=0; i<16; i++) {
    cout << u << " hoch " << setw(2) << i+1 << " : " << v << "   warg = ";
    cout << setw(6) << setprecision(4) << warg(v) << "°\n";
    v *= u;
  }
  cout << "\n";
  cout << " e ^ (2 pi i) = " << ( comp(E,0)^comp(0, 2*PI) ) << '\n';
  cout << " e ^   (pi i) = " << ( E^comp(0, PI) )           << '\n';
  cout << " e ^ (pi/2 i) = " << ( E^comp(0, PI)/comp(2) )   << '\n';
  cout << " e ^ (pi/4 i) = " << ( E^comp(0, PI)/comp(4) )   << '\n';
  cout << " e ^ (pi/8 i) = " << ( E^comp(0, PI)/comp(8) )   << '\n';
  cout << " i ^ i        = " << ( comp(0,1)^comp(0,1) )     << '\n';
  cout << " e ^ (-pi/2)  = " << ( comp(E,0)^comp(-PI/2,0) ) << '\n';
  cout << "-1 ^ 0.5      = " << ( comp(-1,0)^comp(.5,0) )   << '\n';
  cout << "-2i^ 0.5      = " << ( comp(0,-2)^comp(.5,0) )   << '\n';
  cout << " 2i^ 0.5      = " << ( comp(0,2)^comp(.5,0) )    << '\n';
}
