//---------------------------------------------------------------------------
// ŠNikolai V. Shokhirev, 2004-2008  <nikolai@shokhirev.com> http://www.shokhirev.com/nikolai.html
//---------------------------------------------------------------------------

#include "complex.h"

//---------------------------------------------------------------------------

Complex cmplx(real r, real i)
{
  Complex result;
  result.Re = r;
  result.Im = i;
  return result;
}

Complex Real(real r)
{
  return cmplx(r, 0.0);
}

Complex Img(real i)
{
  return cmplx(0.0, i);
}

Complex Cadd(Complex c1, Complex c2)
{
  return cmplx(c1.Re + c2.Re, c1.Im + c2.Im);
}

Complex Cadd(Complex c1, real r)
{
  return cmplx(c1.Re + r, c1.Im);
}

Complex Cadd(real r, Complex c2)
{
  return cmplx(r + c2.Re, c2.Im);
}


Complex Csub(Complex c1, Complex c2)
{
  return cmplx(c1.Re - c2.Re, c1.Im - c2.Im);
}

Complex Csub(Complex c1, real r)
{
  return cmplx(c1.Re - r, c1.Im);
}

Complex Csub(real r, Complex c2)
{
  return cmplx(r - c2.Re, - c2.Im);
}

Complex Cmul(Complex c1, Complex c2)
{
  return cmplx(c1.Re*c2.Re - c1.Im*c2.Im, c1.Re*c2.Im + c1.Im*c2.Re);
}

Complex Cmul(Complex c1, real r)
{
  return cmplx(c1.Re*r,c1.Im*r);
}

Complex Cmul(real r, Complex c2)
{
  return cmplx(c2.Re*r,c2.Im*r);
}

Complex Cdiv(Complex c1, Complex c2)
{
  Complex result;
  result.Re = (c1.Re*c2.Re + c1.Im*c2.Im)/mod2(c2);
  result.Im = (c1.Im*c2.Re - c1.Re*c2.Im)/mod2(c2);
  return result;
}

Complex Cdiv(Complex c1, real r)
{
  return cmplx(c1.Re/r, c1.Im/r);
}

Complex Cdiv(real r, Complex c2)
{
  Complex result;
  result.Re = r*c2.Re/mod2(c2);
  result.Im = -r*c2.Im/mod2(c2);
  return result;
}

Complex conjug(Complex c)
{
  return cmplx(c.Re,-c.Im);
}

real mod2(Complex c)
{
  return (c.Re*c.Re + c.Im*c.Im);
}

real mod(Complex c)
{
  return sqrt(mod2(c));
}

bool Equal(Complex c1, Complex c2, real eps)
{
  return (fabs(c1.Re-c2.Re) < eps && fabs(c1.Im-c2.Im) < eps);
}

bool Equal(Complex c1, real r, real i, real eps)
{
  return (fabs(c1.Re-r) < eps && fabs(c1.Im-i) < eps);
}

bool IsZero(Complex c)
{
  return (c.Re == 0.0) && (c.Im == 0.0);
}

Complex exp(Complex c)
{
  Complex result;
  result.Re = exp(c.Re)*cos(c.Im);
  result.Im = exp(c.Re)*sin(c.Im);
  return result;
}

// result = sqrt(z), result.Re > 0
Complex Csqrt(Complex z)
{
  Complex result;
  result.Re = sqrt(abs(0.5*(mod(z) + z.Re)));
  result.Im = sqrt(abs(0.5*(mod(z) - z.Re)));
  return result;
};