/*********************************************************************
 * This file is part of the cpplibs suite.
 *
 * Copyright (C) 2001 Topi Mäenpää and Jaakko Viertola
 * All rights reserved.
 *
 * This program is free software. You can redistribute and/or modify
 * it under the terms of the free software licence found in the
 * accompanying file "COPYING". The licence terms must always be
 * redistributed with this source file. The above copyright notice
 * must be reproduced in all modified and unmodified copies of this
 * source file.
 *
 * $Revision: 1.12 $
 *********************************************************************/

#ifndef _COMPLEX_H
#define _COMPLEX_H

#include <math.h>
#include "SmartPtr.h"
#include "xml/XMLParser.h"

namespace util
{
	/**
	 * A class for complex calculations. Compared to the STL
	 * correspondent, this class provides more support for mathematical
	 * operations and XML serialization.
	 **/
	template <class T> class Complex : virtual public Object
	{
	public:
		/**
		 * Create a new complex number with the given real and imaginary
		 * parts.
		 **/
		Complex(T r = 0, T i = 0) : real(r), imag(i) {}
		/**
		 * Create a copy of a complex number.
		 **/
		Complex(const Complex& other) : real(other.real), imag(other.imag) {}
		/**
		 * Create a typecast copy of a complex number.
		 **/
		template <class U> Complex(const Complex<U>& other) : real((T)other.real), imag((T)other.imag) {}

		/**
		 * The real part of this complex number.
		 **/
		T real;
		/**
		 * The imaginary part of this complex number.
		 **/
		T imag;

		/**
		 * Convert this number to its complex conjugate.
		 **/
		Complex& conjugate(void) { imag = -imag; return *this; }
		/**
		 * Get a new complex number that is the conjucate of this number.
		 **/
		Complex getConjugate(void) const { return Complex(real,imag); }

		/**
		 * Get the modulus of this complex number (sqrt(real^2 + imag^2)).
		 **/
		T getModulus(void) const { return sqrt(real*real + imag*imag); }
		/**
		 * Get the phase of this complex number (tan^-1(imag/real)).
		 **/
		T getPhase(void) const { return T(atan2(imag,real)); }

		/**
		 * Copy a complex number.
		 **/
		Complex& operator= (const Complex& other) { real=other.real; imag=other.imag; return *this;}
		/**
		 * Compare two complex numbers.
		 **/
		template <class U> inline Complex& operator= (const Complex<U>& other);
		/**
		 * Cast this number to another type.
		 **/
		template <class U> operator Complex<U>() { return Complex<U>(U(real), U(imag)); }
		/**
		 * Insert value T to complex number. Means that the value is given for
		 * the real part of complex number and the imag is set to zero.
		 **/
		Complex& operator= (T value) { real=value;imag=0.0; return *this;}

		/**
		 * Add a complex number to this number.
		 **/
		void operator+= (const Complex& other) { real += other.real; imag += other.imag; }
		/**
		 * Subtract a complex number from this number.
		 **/
		void operator-= (const Complex& other) { real -= other.real; imag -= other.imag; }
		/**
		 * Add a real number to this complex number (affects only real
		 * part):
		 **/
		void operator+= (T value) { real += value; }
		/**
		 * Subtract a real number from this complex number (affects only
		 * real part):
		 **/
		void operator-= (T value) { real -= value; }
		/**
		 * Multiply this complex number by a real value.
		 **/
		void operator*= (T value) { real *= value; imag *= value; }
		/**
		 * Divide this complex number by a real value.
		 **/
		void operator/= (T value) { real /= value; imag /= value; }
		/**
		 * Multiply this complex number by a complex value.
		 **/
		inline void operator*= (const Complex& other);
		/**
		 * Divide this complex number by a complex value.
		 **/
		inline void operator/= (const Complex& other);
	};

	template <class T>
	template <class U> Complex<T>& Complex<T>::operator= (const Complex<U>& other)
	{
		real = T(other.real);
		imag = T(other.imag);
		return *this;
	}

	template <class T> void Complex<T>::operator*= (const Complex<T>& other)
	{
		T tmp = real*other.real-imag*other.imag;
		imag = real*other.imag + imag*other.real;
		real = tmp;
	}
	
	template <class T> void Complex<T>::operator/= (const Complex<T>& other)
	{
		T divisor = other.real * other.real + other.imag * other.imag;
		if(divisor != 0)
			{
				T tmp = (real*other.real + imag*other.imag) / divisor;
				imag = (imag*other.real - real*other.imag) / divisor;
				real = tmp;
			}
		else
			{
				imag=0;
				real=0;
			}
	}

	/**
	 * Calculate the sum of two complex numbers.
	 **/
	template <class T> inline Complex<T> operator+ (const Complex<T>& c1, const Complex<T>& c2)
	{
		Complex<T> result(c1.real+c2.real, c1.imag+c2.imag);
		return result;
	}

	/**
	 * Calculate the difference of two complex numbers.
	 **/
	template <class T> inline Complex<T> operator- (const Complex<T>& c1, const Complex<T>& c2)
	{
		Complex<T> result(c1.real-c2.real, c1.imag-c2.imag);
		return result;
	}
	
	/**
	 * Calculate the product of two complex numbers.
	 **/
	template <class T> inline Complex<T> operator* (const Complex<T>& c1, const Complex<T>& c2)
	{
		Complex<T> result(c1.real*c2.real - c1.imag*c2.imag,
											c1.imag*c2.real + c1.real*c2.imag);
		return result;
	}

	/**
	 * Calculate the quotient of two complex numbers.
	 **/
	template <class T> inline Complex<T> operator/ (const Complex<T>& c1, const Complex<T>& c2)
	{
		T divisor = c2.real * c2.real + c2.imag * c2.imag;
		if(divisor != 0)
			{
				Complex<T> result((c1.real*c2.real + c1.imag*c2.imag) / divisor,
													(c1.imag*c2.real - c1.real*c2.imag) / divisor);
				return result;
			}
		else
			{
				Complex<T> result(0,0);
				return result;
			}
	}
	
	/**
	 * Add a real value to a complex number (affects only real part).
	 **/
	template <class T> inline Complex<T> operator+ (const Complex<T>& c1, T value)
	{
		Complex<T> result(c1.real+value, c1.imag);
		return result;
	}
	/**
	 * Subtract a real value from a complex number (affects only real part).
	 **/
	template <class T> inline Complex<T> operator- (const Complex<T>& c1, T value)
	{
		Complex<T> result(c1.real-value, c1.imag);
		return result;
	}
	/**
	 * Multiply a complex number by a real value.
	 **/
	template <class T> inline Complex<T> operator* (const Complex<T>& c1, T value)
	{
		Complex<T> result(c1.real*value, c1.imag*value);
		return result;
	}
	/**
	 * Divide a complex number by a real value.
	 **/
	template <class T> inline Complex<T> operator/ (const Complex<T>& c1, T value)
	{
		Complex<T> result(c1.real/value, c1.imag/value);
		return result;
	}

	/**
	 * Compare two compled numbers.
	 **/
	template <class T> inline bool operator== (const Complex<T>& c1, const Complex<T>& c2)
	{
		return c1.real == c2.real && c1.imag == c2.imag;
	}

	/**
	 * Compare two compled numbers.
	 **/
	template <class T> inline bool operator!= (const Complex<T>& c1, const Complex<T>& c2)
	{
		return c1.real != c2.real || c1.imag != c2.imag;
	}

	template <class T> std::ostream& operator<< (std::ostream& sout, const Complex<T>& c)
	{
		sout << "<complex real='" << c.real << "' imag='" << c.imag << "'/>";
		return sout;
	}

	template <class T> std::istream& operator>> (std::istream& sin, Complex<T>& c)
	{
		using namespace util::xml;
		
		XMLParser parser;
		sin >> std::ws;
		
		SmartPtr<Node> node(parser.readFragment(sin));
		if (node.get() && node->getNodeType() == Node::ELEMENT_NODE)
			{
				Element* element = (Element*)node.get();
				const Node *child = element->getChildNode("complex.real");
				if (child && child->getNodeType() == Node::ATTRIBUTE_NODE)
					c.real = String::parse<T>(((Attr*)child)->value);
				else
					throw io::IOException("operator>>(istream&,Complex<T>&): input stream does not contain a properly formatted complex number.");
				child = element->getChildNode("complex.imag");
				if (child && child->getNodeType() == Node::ATTRIBUTE_NODE)
					c.imag = String::parse<T>(((Attr*)child)->value);
				else
					throw io::IOException("operator>>(istream&,Complex<T>&): input stream does not contain a properly formatted complex number.");
			}
		return sin;
	}
}

#endif