// Rascal, the Advanced Scientific CALculator
// Copyright (C) 2001, Sebastian Ritterbusch (Rascal@Ritterbusch.de)
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License
// as published by the Free Software Foundation; either version 2
// of the License, or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more detauls.
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
//

#ifndef MCOMPLEX_HPP_INCLUDED
#define MCOMPLEX_HPP_INCLUDED

#include <iostream.h>

template <class T>
class mcomplex
{
   private:
      T r,i;
   public:
      mcomplex<T>(void) {}
      mcomplex<T>(const T & a,const T & b=0) : r(a),i(b) {}
      mcomplex<T>(const mcomplex<T> & x) : r(x.r),i(x.i) {}
      
      const T & real(void) const { return r; }
      const T & imag(void) const { return i; }
      
      friend ostream & operator <<(ostream & o,const mcomplex<T> & a) { o << a.r << ' ' << a.i; return o; }


      friend  inline mcomplex<T> operator -(const mcomplex <T> & a)                     { return mcomplex<T>(-a.r,-a.i); }
      friend  inline mcomplex<T> operator +(const mcomplex <T> & a)                     { return a; }

      friend  inline mcomplex<T> & operator +=(mcomplex<T> & a,const mcomplex<T> & b)    { a.r+=b.r;a.i+=b.i;return a; }
      friend  inline mcomplex<T> & operator +=(mcomplex<T> & a,const T & b)             { a.r+=b;return a; }

      friend  inline mcomplex<T> operator +(const mcomplex<T> & a,const mcomplex<T> & b) { return mcomplex<T>(a.r+b.r,a.i+b.i); }
      friend  inline mcomplex<T> operator +(const T & a,const mcomplex<T> & b)          { return mcomplex<T>(a+b.r,b.i); }
      friend  inline mcomplex<T> operator +(const mcomplex<T> & b,const T & a)          { return mcomplex<T>(b.r+a,b.i); }

      friend  inline mcomplex<T> & operator -=(mcomplex<T> & a,const mcomplex<T> & b)    { a.r-=b.r;a.i-=b.i;return a; }
      friend  inline mcomplex<T> & operator -=(mcomplex<T> & a,const T & b)             { a.r-=b;return a; }

      friend  inline mcomplex<T> operator -(const mcomplex<T> & a,const mcomplex<T> & b) { return mcomplex<T>(a.r-b.r,a.i-b.i); }
      friend  inline mcomplex<T> operator -(const T & a,const mcomplex<T> & b)          { return mcomplex<T>(a-b.r,-b.i); }
      friend  inline mcomplex<T> operator -(const mcomplex<T> & b,const T & a)          { return mcomplex<T>(b.r-a,b.i); }

      friend  inline mcomplex<T> & operator *=(mcomplex<T> & a,const mcomplex<T> & b)    { T r=a.r*b.r-a.i*b.i;a.i=a.r*b.i+a.i*b.r;a.r=r;return a; }
      friend  inline mcomplex<T> & operator *=(mcomplex<T> & a,const T & b)             { a.r*=b;a.i*=b; return a; }

      friend  inline mcomplex<T> operator *(const mcomplex<T> & a,const mcomplex<T> & b) { return mcomplex<T>(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r); }
      friend  inline mcomplex<T> operator *(const T & a,const mcomplex<T> & b)          { return mcomplex<T>(a*b.r,a*b.i); }
      friend  inline mcomplex<T> operator *(const mcomplex<T> & b,const T & a)          { return mcomplex<T>(b.r*a,b.i*a); }

      friend  inline mcomplex<T> & operator /=(mcomplex<T> & a,const mcomplex<T> & b)    { T q=b.r*b.r+b.i*b.i,r=(a.r*b.r+a.i*b.i)/q;a.i=(a.i*b.r-a.r*b.i)/q;a.r=r;return a; }
      friend  inline mcomplex<T> & operator /=(mcomplex<T> & a,const T & b)             { a.r/=b;a.i/=b; return a; }

      friend  inline mcomplex<T> operator /(const mcomplex<T> & a,const mcomplex<T> & b) { T q=b.r*b.r+b.i*b.i; return mcomplex<T>((a.r*b.r+a.i*b.i)/q,(a.i*b.r-a.r*b.i)/q); }
      friend  inline mcomplex<T> operator /(const T & a,const mcomplex<T> & b)          { T q=b.r*b.r+b.i*b.i; return mcomplex<T>((a*b.r)/q,(-a*b.i)/q); }
      friend  inline mcomplex<T> operator /(const mcomplex<T> & b,const T & a)          { return mcomplex<T>(b.r/a,b.i/a); }


      friend  inline int operator ==(const mcomplex<T> & a,const mcomplex<T> & b)   { return (a.r==b.r) && (a.i==b.i); }
      friend  inline int operator !=(const mcomplex<T> & a,const mcomplex<T> & b)   { return (a.r!=b.r) || (a.i!=b.i); }

      friend  inline T abs(const mcomplex<T> & a) { return sqrt((a.r*a.r)+(a.i*a.i)); }
      friend  inline T arg(const mcomplex<T> & a) 
      { 
         if(a.r>0)
            return asin(a.i/abs(a)); 
         else
            if(a.i>0)
               return 4*atan(T(1))-asin(a.i/abs(a));
            else if(a.i<0)
               return -4*atan(T(1))-asin(a.i/abs(a));
            else if(abs(a)==0)
               return T(0);
            else
               return 4*atan(T(1));
      }
      

      // Standard functions

      friend inline mcomplex<T> exp(const mcomplex<T> & a)
      {
         // e^(x+iy)=e^x * (cos(y)+i*sin(y))
         T ex(exp(a.r));
         return mcomplex<T>(ex*cos(a.i),ex*sin(a.i));
      }
      
      friend inline mcomplex<T> log(const mcomplex<T> & a)
      {
         // ln(x+iy)=ln(abs(x+iy))+i*arg(x+iy)
         return mcomplex<T>(log(abs(a)),arg(a));
      }
      
      friend inline mcomplex<T> pow(const mcomplex<T> & a,const mcomplex<T> & b)
      {
         // a^b=e^(b*ln(a));
         return exp(b*log(a));
      }
      
      friend inline mcomplex<T> sqrt(const mcomplex<T> & a)
      {
         // really slow... using polar representation would be better
         return pow(a,mcomplex<T>(.5));
      }

      friend inline mcomplex<T> sqr(const mcomplex<T> & a)
      {
         return mcomplex<T>(sqr(a.r)-sqr(a.i),2*a.r*a.i);
      }

      friend inline mcomplex<T> sin(const mcomplex<T> & a)
      {
         // sin(x+iy)=sin(x)*cosh(y)+i*cos(x)*sinh(y)
         return mcomplex<T>(sin(a.r)*cosh(a.i),cos(a.r)*sinh(a.i));
      }

      friend inline mcomplex<T> cos(const mcomplex<T> & a)
      {
         // cos(x+iy)=cos(x)*cosh(y)-i*sin(x)*sinh(y)
         return mcomplex<T>(cos(a.r)*cosh(a.i),-sin(a.r)*sinh(a.i));
      }
      
      friend inline mcomplex<T> tan(const mcomplex<T> & a)
      {
         // tan(x+iy)=(sin(2*x)+i*sinh(2*y))/(cos(2*x)+cosh(2*y))
         T q(cos(2*a.r)+cosh(2*a.i));
         return mcomplex<T>(sin(a.r+a.r)/q,sinh(a.i+a.i)/q);
      }
      
      friend inline mcomplex<T> cot(const mcomplex<T> & a)
      {
         return T(1)/tan(a);
      }
      
      friend inline mcomplex<T> sinh(const mcomplex<T> & a)
      {
         // sinh(x+iy)=sin(x)*cos(y)+i*cosh(x)*sin(y)
         return mcomplex<T>(sinh(a.r)*cos(a.i),cosh(a.r)*sin(a.i));
      }

      friend inline mcomplex<T> cosh(const mcomplex<T> & a)
      {
         // cosh(x+iy)=cosh(x)*cos(y)+i*sinh(x)*sin(y)
         return mcomplex<T>(cosh(a.r)*cos(a.i),sinh(a.r)*sin(a.i));
      }
      
      friend inline mcomplex<T> tanh(const mcomplex<T> & a)
      {
         // tan(x+iy)=(sinh(2*x)+i*sin(2*y))/(cosh(2*x)+cos(2*y))
         T q(cosh(a.r+a.r)+cos(a.i+a.i));
         return mcomplex<T>(sinh(a.r+a.r)/q,sin(a.i+a.i)/q);
      }
      
      friend inline mcomplex<T> coth(const mcomplex<T> & a)
      {
         return T(1)/tan(a);
      }
      
      friend inline mcomplex<T> asin(const mcomplex<T> & z)
      {
         // arcsin(a)=-i*ln(i*z+sqrt(1-z*z))
         mcomplex<T> i(0,1);
         return -i*log(i*z+sqrt(1-z*z));
      }      
      
      friend inline mcomplex<T> acos(const mcomplex<T> & z)
      {
         // arccos(a)=-i*ln(z+sqrt(z*z-1))
         mcomplex<T> i(0,1);
         return -i*log(z+sqrt(z*z-1));
      }      
      
      friend inline mcomplex<T> atan(const mcomplex<T> & z)
      {
         // arctan(a)=(1/2*i)*ln((1+i*z)/(1-i*z))
         mcomplex<T> i(0,1);
         return (1/(2*i))*log((1+i*z)/(1-i*z));
      }
      
      friend inline mcomplex<T> acot(const mcomplex<T> & z)
      {
         // arccot(a)=-(1/2*i)*ln((i*z+1)/(i*z-1))
         mcomplex<T> i(0,1);
         return -(1/(2*i))*log((i*z+1)/(i*z-1));
      }
      
      friend inline mcomplex<T> asinh(const mcomplex<T> & z)
      {
         // arsinh(z)=ln(z+sqrt(z*z+1))
         return log(z+sqrt(z*z+1));
      }
      
      friend inline mcomplex<T> acosh(const mcomplex<T> & z)
      {
         // arcosh(z)=ln(z+sqrt(z*z-1))
         return log(z+sqrt(z*z+1));
      }
      
      friend inline mcomplex<T> atanh(const mcomplex<T> & z)
      {
         // artanh(z)=ln((1+z)/(1-z))/2
         return log((1+z)/(1-z))/2;
      }
      
      friend inline mcomplex<T> acoth(const mcomplex<T> & z)
      {
         // arcoth(z)=ln((z+1)/(z-1))/2
         return log((z+1)/(z-1))/2;
      }
};


#endif