// complex.cc: Complex number class and polynomial evaluation
//
// Written by Michael J. Gourlay
// Fall 1996, University of Colorado at Boulder
// Algorithms
//
// C++ code


/*

Top level description:
----------------------

A class for complex numbers is implemented, including several typical
mathematical operators.

*/

#include <iostream.h>
#include <stdlib.h>
#include <math.h>



// Complex number class:
//
// This could use templates to make the class any available precision.
//
class Complex {
  public:
    // Default constructor: Do nothing extra.
    //
    // I want to do nothing because I do not believe that a new
    // instance should automatically be initialized to some value.
    // Some might say that I should assign the members to zero, but
    // this is inconsistent with the way primative types are
    // handled.  Also, for arrays, I do not want to waste time
    // initializing the values of a bazillion Complex numbers when
    // those values are not going to be used anyway.  The client
    // should have to explicity addign values if he wants them.
    // There is nothing holding the client back from initializing
    // explicit values if they are desired, but this constructor
    // will not do it for them.
    Complex(void) {};

    // Initializing constructor:
    //   Real and imaginary parts may be provided.
    //   Imaginary part is optional.
    //   This constructor also effectively promotes doubles to complex.
    Complex(double r, double i=0.0) : re(r), im(i) {};


    // ---- Mathematical operators ----
    Complex& operator += (const Complex&);
    Complex& operator *= (const Complex&);

    // conj: complex conjugate z* = (z_r - i z_i)
    Complex conj() const { return Complex(re, -im); }

    // arg: "argument" of a complex number, which is the phase
    // This follows from the Euler formula:
    //       exp(z) = |z| (cos(phi) + i sin(phi))
    //   where |z| = (Re^2(z) + Im^2(z))^(1/2) is the Euclidean (2-) norm of z
    //   and where phi = arg(z) =atan(Im/Re) is the argument of z.
    //
    // This is tricky and subtle because arg is not a unique value.
    // In other words, the natural log of a complex number is
    // multi-valued.  This leads to branch cuts in complex analysis.
    double arg() const { return atan2(im, re); }

    // norm: 2-norm of a complex number (the "absolute value")
    // |z| = (z z*)^(1/2) = ((z_r + i z_i) * (z_r - i z_i))^(1/2)
    double norm() const { return sqrt(re*re + im*im); }

    // mag2: square magnitude of a complex number
    // |z|^2 = z z* = (z_r + i z_i) * (z_r - i z_i)
    double mag2() const { return (re*re + im*im); }

    // Accessor functions:
    double real() const { return re; }
    double imag() const { return im; }

  private:
    double re;  // real part
    double im;  // imaginary part
};


// ---- Member methods: ----
// =========================


// Cumulative add operator
//
inline Complex&
Complex::operator += (const Complex& z)
{
  re += z.re;
  im += z.im;
  return *this;
}




// Cumulative multiply operator
//
inline Complex&
Complex::operator *= (const Complex& z)
{
  double r = re * z.re - im * z.im;
  im = re * z.im + im * z.re;
  re = r;
  return *this;
}


// ---- Global methods ----
// ========================


// All of these operators are global (instead of friends) because
// they do not need write access to private members,
// and the private members they need access to have public accessors.


// Complex binary operator + : add two complex numbers.
//
// This routine also elegantly adds a complex number to a float,
// where either addend is complex.  This works because there is a
// Complex constructor which will promote a double to a Complex.
// When is comes time to add a double to a complex, this is the only
// operator which will handle complex numbers, so the machine will
// look for a way to promote or cast a double to a Complex, and the
// Complex constructor does that, so this operator will be called.
// (For some whiney compilers, this will cause a warning:
// value copied to temporary, reference to temporary used)
//
// The Complex constructor is used to create a new Complex instance
// for the return value, since the copy constructor will be used at
// assignment time to copy from the result of this operator.  That
// is, this operator returns an object, not a reference to an
// object.  This is a subtlety because it is somewhat unusual, as
// object operators go.
//
inline Complex
operator + (const Complex& z1, const Complex& z2)
{
  return Complex(z1.real() + z2.real(), z1.imag() + z2.imag()) ;
}



// Complex binary operator -: subtract two complex numbers.
inline Complex
operator - (const Complex& z1, const Complex& z2)
{
  return Complex(z1.real() - z2.real(), z1.imag() - z2.imag()) ;
}




// Complex operator * : multiply two complex numbers.
//
// result = z1 * z2
//        = (z1.re + i z1.im) * (z2.re + i z2.im)
//        = (z1.re * z2.re - z1.im * z2.im)
//          + i * (z1.re * z2.im + z1.im * z2.re)
//
// As for the + operator, this will work for doubles multiplying Complex.
//
inline Complex
operator * (const Complex& z1, const Complex& z2)
{
  return Complex(z1.real() * z2.real() - z1.imag() * z2.imag(),
                 z1.real() * z2.imag() + z1.imag() * z2.real() );
}




// Complex operator / : divide a complex number by a double.
//
// This is reduntant with the complex/complex operator,
// but this version is more efficient.
//
inline Complex
operator / (const Complex& znum, const double div)
{
  return Complex(znum.real() / div, znum.imag() / div );
}




// Complex operator / : divide a complex number by a complex number
//
inline Complex
operator / (const Complex& znum, const Complex& zdiv)
{
  return znum * zdiv.conj() / zdiv.mag2();
}




// exp: natural (base e) exponentiation: e^z
//
// Uses exp(z) = exp(z_r + i z_i) = exp(z_r) * exp(i z_i)
//   where exp(i z_i) = cos(z_i) + i sin(z_i) by the Euler formula
//
Complex
exp(const Complex& z)
{
  const double exp_real = exp(z.real());
  return Complex(exp_real * cos(z.imag()), exp_real * sin(z.imag()));
}




// log: natural logarithm
//
// uses ln(z) = ln(|z| * exp(i phi)) = ln|z| + ln(exp(i phi)) = ln|z| + i phi
//   where |z| = 2-norm of z
//   where phi = arg(z)
//
Complex
log(const Complex& z)
{
  return Complex(log(z.norm()), z.arg());
}




// pow: raise a complex number to a real power: z^y
//
// uses z^y = exp(ln(z^y)) = exp(y ln(z))
//
inline Complex
pow(const Complex& z, const double y)
{
  return exp(y * log(z));
}




// input stream operator: input a complex number
//
// Input can either be a real floating point value,
// or a parenthetical comma separated pair: (real,imag)
//
// Acceptable formats:
//   real
//   ( real )
//   ( real , imag )
//
istream&
operator >> (istream& is, Complex& zz)
{
  double real_in;       // real part.
  double imag_in=0.0;   // imaginary part. default zero in case not provided
  char char_in='\0';    // input char for parens and comma.  null by default.

  is >> char_in;

  if(char_in == '(') {
    is >> real_in >> char_in;
    if(char_in == ',') {
      is >> imag_in >> char_in;
    }
    if (char_in != ')') {
      is.clear(ios::badbit);
    }
  } else {
    is.putback(char_in);
    is >> real_in;
  }

  // Compose the complex value
  if(is) {
    zz = Complex(real_in, imag_in);
  }

  return is;
}




// output stream operator: print a complex number
ostream&
operator << (ostream& os, const Complex& zz)
{
  return os << '(' << zz.real() << ',' << zz.imag() << ')';
}




// polyRead: read a list of polynomial coefficients
//
// is: input stream to read coefficients from
// polynomial_order (out): order of the polynomial
// zc (out): polynomial coefficients -- allocated here -- must free later
//
void
polyRead(istream& is, int *polynomial_order, Complex **zc)
{
  // ---- Read the order of the polynomial ----
  is >> *polynomial_order ;

  if(!is) {
    cout << "\npolyRead: Bad input while reading order\n" ;
    exit(1);
  }

  // ---- Read the coefficients ----

  int coeff_idx;                           // Coefficient index
  *zc = new Complex[*polynomial_order+1];  // complex coefficients

  for(coeff_idx=0; coeff_idx <= *polynomial_order; coeff_idx++) {
    // Read the value
    is >> (*zc)[coeff_idx];

    if(!is) {
      cout << "\npolyRead: Bad input while reading coefficient\n" ;
      exit(2);
    }
  }
}




// polyEval: evaluate a polynomial using the Horner rule
//
// zc (in): list of complex coefficients, where zc[i] is the coefficient of x^i
// order (in): order of the polynomial: zc[i] : 0 <= i <= order
// abcissa (in): evaluation place
//
Complex
polyEval(const Complex zc[], int order, const Complex& abcissa)
{
  int ci;
  Complex result(0.0,0.0);

  // Evaluate polynomial
  for(ci = order; ci >0; ci --) {
    result += zc[ci];
    result *= abcissa;
  }
  // Add constant term
  result += zc[0];

  return result;
}




int
main(int argc, char **argv)
{


  const double pi = 3.14159265358979323844;
  const Complex ii(0.0, 1.0);

#ifdef ARITH_TESTS
  // ---- Test data variables ----
  const double ee = 2.71828182845904523536;

  Complex z1;
  Complex z2(2.0);
  Complex z3(3.0, pi);

  // ---- A bunch of tests of Complex operators ----

  cout << "--------------------\n";
  cout << "Real variable values\n" ;
  cout << "--------------------\n";

  cout << "pi=" << pi << '\n';
  cout << "ee=" << ee << '\n';

  cout << "-----------------------\n";
  cout << "Complex variable values\n" ;
  cout << "-----------------------\n";

  cout << "z1=" << z1 << '\n';  // Not initialized -- should be garbage
  cout << "z2=" << z2 << '\n';
  cout << "z3=" << z3 << '\n';
  cout << "ii=" << ii << '\n';

  cout << "----------\n";
  cout << "   Sums   \n";
  cout << "----------\n";

  // Sum of complex variables
  z1 = z2 + z3;
  cout << "z1 = z2 + z3 =" << z1 << '\n';

  // Sum of complex literals
  z1 = Complex(1.2, 3.4) + Complex(5.6, 7.8);
  cout << "z1 = (1.2,3.4) + (5.6,7.8) =" << z1 << '\n';

  // Sum of real variable + complex variable
  z1 = pi + z3;
  cout << "z1 = pi + z3 =" << z1 << '\n';

  // Sum of complex variable + real variable
  z1 = z2 + ee;
  cout << "z1 = z2 + ee =" << z1 << '\n';

  // Sum of real literal + complex variable
  z1 = 4.0 + z3;
  cout << "z1 = 4.0 + z3 =" << z1 << '\n';

  // Sum of complex variable + real literal
  z1 = z2 + 5.0;
  cout << "z1 = z2 + 5.0 =" << z1 << '\n';

  cout << "---------------\n";
  cout << "   Self-Sums   \n";
  cout << "---------------\n";

  // Sum of complex variable
  z1 += z2 ;
  cout << "z1 += z2 =" << z1 << '\n';

  // Sum of complex literal
  z1 += Complex(1.2, 3.4);
  cout << "z1 += (1.2,3.4) =" << z1 << '\n';

  // Sum of real variable
  z1 += pi ;
  cout << "z1 += pi =" << z1 << '\n';

  // Sum of real literal
  z1 += 4.0;
  cout << "z1 += 4.0 =" << z1 << '\n';

  cout << "------------------\n";
  cout << "   Self-Product   \n";
  cout << "------------------\n";

  // Product of complex variable
  z1 *= z2 ;
  cout << "z1 += z2 =" << z1 << '\n';

  // Product of complex literal
  z1 += Complex(1.2, 3.4);
  cout << "z1 *= (1.2,3.4) =" << z1 << '\n';

  // Product of real variable
  z1 *= pi ;
  cout << "z1 *= pi =" << z1 << '\n';

  // Product of real literal
  z1 *= 4.0;
  cout << "z1 *= 4.0 =" << z1 << '\n';

  cout << "-----------------\n";
  cout << "   Differences   \n";
  cout << "-----------------\n";

  // Difference of complex variables
  z1 = z2 - z3;
  cout << "z1 = z2 - z3 =" << z1 << '\n';

  // Difference of complex literals
  z1 = Complex(1.2, 3.4) - Complex(5.6, 7.8);
  cout << "z1 = (1.2,3.4) - (5.6,7.8) =" << z1 << '\n';

  // Difference of real variable - complex variable
  z1 = pi - z3;
  cout << "z1 = pi - z3 =" << z1 << '\n';

  // Difference of complex variable - real variable
  z1 = z2 - ee;
  cout << "z1 = z2 - ee =" << z1 << '\n';

  // Difference of real literal - complex variable
  z1 = 4.0 - z3;
  cout << "z1 = 4.0 - z3 =" << z1 << '\n';

  // Difference of complex variable - real literal
  z1 = z2 - 5.0;
  cout << "z1 = z2 - 5.0 =" << z1 << '\n';

  cout << "--------------\n";
  cout << "   Products   \n";
  cout << "--------------\n";

  // Product of complex variables
  z1 = z2 * z3;
  cout << "z1 = z2 * z3 =" << z1 << '\n';

  // Product of complex literals
  z1 = Complex(1.2, 3.4) * Complex(5.6, 7.8);
  cout << "z1 = (1.2,3.4) * (5.6,7.8) =" << z1 << '\n';

  // Product of real variable * complex variable
  z1 = pi * z3;
  cout << "z1 = pi * z3 =" << z1 << '\n';

  // Product of complex variable * real variable
  z1 = z2 * ee;
  cout << "z1 = z2 * ee =" << z1 << '\n';

  // Product of real literal * complex variable
  z1 = 4.0 * z3;
  cout << "z1 = 4.0 * z3 =" << z1 << '\n';

  // Product of complex variable * real literal
  z1 = z2 * 5.0;
  cout << "z1 = z2 * 5.0 =" << z1 << '\n';

  cout << "---------------\n";
  cout << "   Quotients   \n";
  cout << "---------------\n";

  // Quotient of complex variables
  z1 = z2 / z3;
  cout << "z1 = z2 / z3 =" << z1 << '\n';

  // Quotient of complex variables
  z1 = z2 / z2;
  cout << "z1 = z2 / z2 =" << z1 << '\n';

  // Quotient of complex variables
  z1 = z3 / z3;
  cout << "z1 = z3 / z3 =" << z1 << '\n';

  // Quotient of complex literals
  z1 = Complex(1.2, 3.4) / Complex(5.6, 7.8);
  cout << "z1 = (1.2,3.4) / (5.6,7.8) =" << z1 << '\n';

  // Quotient of real variable / complex variable
  z1 = pi / z3;
  cout << "z1 = pi / z3 =" << z1 << '\n';

  // Quotient of complex variable / real variable
  z1 = z2 / ee;
  cout << "z1 = z2 / ee =" << z1 << '\n';

  // Quotient of real literal / complex variable
  z1 = 4.0 / z3;
  cout << "z1 = 4.0 / z3 =" << z1 << '\n';

  // Quotient of complex variable / real literal
  z1 = z2 / 5.0;
  cout << "z1 = z2 / 5.0 =" << z1 << '\n';

  cout << "-------------------\n";
  cout << "   Miscellaneous   \n";
  cout << "-------------------\n";

  // Conjugate
  cout << "conj(z1) = " << z1.conj() << '\n' ;
  cout << "conj(z2) = " << z2.conj() << '\n' ;
  cout << "conj(z3) = " << z3.conj() << '\n' ;
  cout << "conj(ii) = " << ii.conj() << '\n' ;

  // Argument
  cout << "arg(z1) = " << z1.arg() << '\n' ;
  cout << "arg(z2) = " << z2.arg() << '\n' ;
  cout << "arg(z3) = " << z3.arg() << '\n' ;
  cout << "arg(ii) = " << ii.arg() << '\n' ;

  // Norm
  cout << "abs(z1) = " << z1.norm() << '\n' ;
  cout << "abs(z2) = " << z2.norm() << '\n' ;
  cout << "abs(z3) = " << z3.norm() << '\n' ;
  cout << "abs(ii) = " << ii.norm() << '\n' ;

  // Square Magnitude
  cout << "mag2(z1) = " << z1.mag2() << '\n' ;
  cout << "mag2(z2) = " << z2.mag2() << '\n' ;
  cout << "mag2(z3) = " << z3.mag2() << '\n' ;
  cout << "mag2(ii) = " << ii.mag2() << '\n' ;

  // Exponent
  cout << "exp(z1) = " << exp(z1) << '\n' ;
  cout << "exp(z2) = " << exp(z2) << '\n' ;
  cout << "exp(z3) = " << exp(z3) << '\n' ;
  cout << "exp(ii) = " << exp(ii) << '\n' ;

  // Log
  cout << "log(z1) = " << log(z1) << '\n' ;
  cout << "log(z2) = " << log(z2) << '\n' ;
  cout << "log(z3) = " << log(z3) << '\n' ;
  cout << "log(ii) = " << log(ii) << '\n' ;

  // Power
  cout << "pow(z1, ee) = " << pow(z1, ee) << '\n' ;
  cout << "pow(z2, pi) = " << pow(z2, pi) << '\n' ;
  cout << "pow(z3, ee) = " << pow(z3, ee) << '\n' ;
  cout << "pow(ii, pi) = " << pow(ii, pi) << '\n' ;


  // ---- Other tests ----
  // Should fail: there is no casting function from Complex to double
  // double d1 = z1;  // Fails to compile, as it should

#endif ARITH_TESTS


  // ---- Read polynomial coefficients ----

  int polynomial_order;
  Complex *zc;

  polyRead(cin, &polynomial_order, &zc);

  // Print out the coefficients to assure user that input was read correctly
  int coeff_idx;
  for(coeff_idx=0; coeff_idx <= polynomial_order; coeff_idx++) {
    cout << "coefficient " << coeff_idx << " == " << zc[coeff_idx] << '\n';
  }

  cout << "\n-------------------------------------\n";


  // ---- Read the evaluation points for the polynomial ----
  Complex abcissa_start;
  Complex abcissa_end;
  int num_evals;

  cout << "\nEnter the (complex) polynomial evaluation start abcissa: ";
  cin >> abcissa_start;
  if(!cin) exit(3);  // handle bad input

  cout << "\nEnter the (complex) polynomial evaluation end abcissa:   ";
  cin >> abcissa_end;
  if(!cin) exit(4);  // handle bad input

  cout << "\nEnter the number of evaluation points:                   ";
  cin >> num_evals;
  if(!cin) exit(5);  // handle bad input

  cout << "\n-------------------------------------\n";


  // ---- Evaluate the polynomial ----

  int eval_iter;
  Complex poly_value;
  Complex abcissa;
  const Complex abcissa_range = abcissa_end - abcissa_start;

  for(eval_iter = 0; eval_iter < num_evals ; eval_iter++) {
    abcissa = abcissa_start + abcissa_range * eval_iter / (num_evals - 1);
    poly_value = polyEval(zc, polynomial_order, abcissa);
    cout << "P(" << abcissa << ")\t= " << poly_value << '\n';
  }

  cout << "\n-------------------------------------\n";


  // ---- Evaluate the polynomial at complex roots of one ----
  // This is actually a Discrete Fourier Transform of the polynomial

  cout << "\nPolynomial evaluated at complex roots of one (DFT of P):\n";
  cout <<   "--------------------------------------------------------\n";

  //int eval_iter;
  //Complex poly_value;
  //Complex abcissa;

  for(eval_iter = 0; eval_iter < 2*polynomial_order ; eval_iter++) {
    abcissa = exp(pi * ii * eval_iter / polynomial_order);
    poly_value = polyEval(zc, polynomial_order, abcissa);
    cout << "P(" << abcissa << ")\t= " << poly_value << '\n';
  }


  // Clean up and exit
  delete [] zc;

  return 0;
}