// 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
syntax highlighted by Code2HTML, v. 0.9.1