/* fcomplex.h */

/* C++ program (class) for operations with complex numbers
   For a and b complex, and x,y real:
   Absolute value of a: a.Cabs() 
   Conjugate of a: a.Conjg() 
   Square root of a: a.Csqrt()
   Real part of a: a.r()
   Imaginary part of a: a.i()
   Transforming x,y pair into complex: a = fcomplex(x, y) 
   Adding a + b: a.Cadd(b)
   Subtracting a - b: a.Csub(b)
   Multiplying a * b: a.Cmul(b)
   Dividing a / b: a.Cdiv(b)
   multiplying by a real number a * x: a.RCmul(a)
   written by C.A. Bertulani - 03/08/2000     */ 

#include <math.h>

typedef double Number;

class fcomplex{
public:

	fcomplex(Number r, Number i);
	fcomplex();
	Number r();
	Number i();
	Number Cabs();
	fcomplex Conjg();
	fcomplex Csqrt();
	fcomplex Cadd(fcomplex);
	fcomplex Csub(fcomplex);
	fcomplex Cmul(fcomplex);
	fcomplex Cdiv(fcomplex);
	fcomplex RCmul(Number);
private:
	Number the_r;
	Number the_i;
};

fcomplex::fcomplex(Number r, Number i){
	the_r = r;
	the_i = i;
}

fcomplex::fcomplex(){
}

Number fcomplex::r(){return the_r;}

Number fcomplex::i(){return the_i;}

fcomplex fcomplex::Cadd(fcomplex a)
{	Number x; Number y;

	x = the_r + a.the_r;
	y = the_i + a.the_i;
	return fcomplex(x,y);
}

fcomplex fcomplex::Csub(fcomplex a)
{	Number x; Number y;
	x = the_r - a.the_r;
	y = the_i - a.the_i;
	return fcomplex(x,y);
}


fcomplex fcomplex::Cmul(fcomplex a)
{   Number x, y; 	
	x = the_r*a.the_r - the_i*a.the_i;
	y = the_i*a.the_r + the_r*a.the_i;
	return fcomplex(x,y);
}


fcomplex fcomplex::Conjg()
{	 
	return fcomplex(the_r, - the_i);
}

fcomplex fcomplex::Cdiv(fcomplex a)
{	fcomplex c;	
	Number r,den;
	if (fabs(the_r) >= fabs(the_i)) {
		r = the_i / the_r;
		den = the_r + r*the_i;
		c.the_r = (a.the_r + r*a.the_i) / den;
		c.the_i=(a.the_i-r*a.the_r) / den;
	} else {
		r = the_r/the_i;
		den = the_i+r*the_r;
		c.the_r=(a.the_r*r + a.the_i) / den;
		c.the_i=(a.the_i*r - a.the_r) / den;
	}
	return c;
}

Number fcomplex::Cabs()
{	Number x,y,ans,temp;
	x=fabs(the_r);
	y=fabs(the_i);
	if (x == 0.0)
		ans=y;
	else if (y == 0.0)
		ans=x;
	else if (x > y) {
		temp=y/x;
		ans=x*sqrt(1.0+temp*temp);
	} else {
		temp=x/y;
		ans=y*sqrt(1.0+temp*temp);
	}
	return ans;
}

fcomplex fcomplex::Csqrt()
{	fcomplex c;
	Number x,y,w,r;
	if ((the_r == 0.0) && (the_i == 0.0)) {
		c.the_r=0.0;
		c.the_i=0.0;
		return c;
	} else {
		x=fabs(the_r);
		y=fabs(the_i);
		if (x >= y) {
			r=y/x;
			w=sqrt(x)*sqrt(0.5*(1.0+sqrt(1.0+r*r)));
		} else {
			r=x/y;
			w=sqrt(y)*sqrt(0.5*(r+sqrt(1.0+r*r)));
		}
		if (the_r >= 0.0) {
			c.the_r=w;
			c.the_i=the_i/(2.0*w);
		} else {
			c.the_i=(the_i >= 0) ? w : -w;
			c.the_r=the_i/(2.0*c.the_i);
		}
		return c;
	}
}

fcomplex fcomplex::RCmul(Number x)
{	fcomplex c;
	c.the_r=x*the_r;
	c.the_i=x*the_i;
	return c;
}