#include <complex>
extern "C++" {
template <class FLOAT> complex<FLOAT>
cos (const complex<FLOAT>& x)
{
return complex<FLOAT> (cos (real (x)) * cosh (imag (x)),
- sin (real (x)) * sinh (imag (x)));
}
template <class FLOAT> complex<FLOAT>
cosh (const complex<FLOAT>& x)
{
return complex<FLOAT> (cosh (real (x)) * cos (imag (x)),
sinh (real (x)) * sin (imag (x)));
}
template <class FLOAT> complex<FLOAT>
exp (const complex<FLOAT>& x)
{
return polar (FLOAT (exp (real (x))), imag (x));
}
template <class FLOAT> complex<FLOAT>
log (const complex<FLOAT>& x)
{
return complex<FLOAT> (log (abs (x)), arg (x));
}
template <class FLOAT> complex<FLOAT>
pow (const complex<FLOAT>& x, const complex<FLOAT>& y)
{
FLOAT logr = log (abs (x));
FLOAT t = arg (x);
return polar (FLOAT (exp (logr * real (y) - imag (y) * t)),
FLOAT (imag (y) * logr + real (y) * t));
}
template <class FLOAT> complex<FLOAT>
pow (const complex<FLOAT>& x, FLOAT y)
{
return exp (FLOAT (y) * log (x));
}
template <class FLOAT> complex<FLOAT>
pow (FLOAT x, const complex<FLOAT>& y)
{
return exp (y * FLOAT (log (x)));
}
template <class FLOAT> complex<FLOAT>
sin (const complex<FLOAT>& x)
{
return complex<FLOAT> (sin (real (x)) * cosh (imag (x)),
cos (real (x)) * sinh (imag (x)));
}
template <class FLOAT> complex<FLOAT>
sinh (const complex<FLOAT>& x)
{
return complex<FLOAT> (sinh (real (x)) * cos (imag (x)),
cosh (real (x)) * sin (imag (x)));
}
#include <iostream.h>
template <class FLOAT> istream&
operator >> (istream& is, complex<FLOAT>& x)
{
FLOAT re, im = 0;
char ch = 0;
if (is.ipfx0 ())
{
if (is.peek () == '(')
is >> ch;
is >> re;
if (ch == '(')
{
is >> ch;
if (ch == ',')
is >> im >> ch;
}
}
is.isfx ();
if (ch != 0 && ch != ')')
is.setstate (ios::failbit);
else if (is.good ())
x = complex<FLOAT> (re, im);
return is;
}
template <class FLOAT> ostream&
operator << (ostream& os, const complex<FLOAT>& x)
{
return os << '(' << real (x) << ',' << imag (x) << ')';
}
template <class FLOAT> complex<FLOAT>&
__doadv (complex<FLOAT>* ths, const complex<FLOAT>& y)
{
FLOAT ar = abs (y.re);
FLOAT ai = abs (y.im);
FLOAT nr, ni;
FLOAT t, d;
if (ar <= ai)
{
t = y.re / y.im;
d = y.im * (1 + t*t);
nr = (ths->re * t + ths->im) / d;
ni = (ths->im * t - ths->re) / d;
}
else
{
t = y.im / y.re;
d = y.re * (1 + t*t);
nr = (ths->re + ths->im * t) / d;
ni = (ths->im - ths->re * t) / d;
}
ths->re = nr;
ths->im = ni;
return *ths;
}
template <class FLOAT> complex<FLOAT>
operator / (const complex<FLOAT>& x, const complex<FLOAT>& y)
{
FLOAT ar = abs (real (y));
FLOAT ai = abs (imag (y));
FLOAT nr, ni;
FLOAT t, d;
if (ar <= ai)
{
t = real (y) / imag (y);
d = imag (y) * (1 + t*t);
nr = (real (x) * t + imag (x)) / d;
ni = (imag (x) * t - real (x)) / d;
}
else
{
t = imag (y) / real (y);
d = real (y) * (1 + t*t);
nr = (real (x) + imag (x) * t) / d;
ni = (imag (x) - real (x) * t) / d;
}
return complex<FLOAT> (nr, ni);
}
template <class FLOAT> complex<FLOAT>
operator / (FLOAT x, const complex<FLOAT>& y)
{
FLOAT ar = abs (real (y));
FLOAT ai = abs (imag (y));
FLOAT nr, ni;
FLOAT t, d;
if (ar <= ai)
{
t = real (y) / imag (y);
d = imag (y) * (1 + t*t);
nr = x * t / d;
ni = -x / d;
}
else
{
t = imag (y) / real (y);
d = real (y) * (1 + t*t);
nr = x / d;
ni = -x * t / d;
}
return complex<FLOAT> (nr, ni);
}
template <class FLOAT> complex<FLOAT>
pow (const complex<FLOAT>& xin, int y)
{
if (y == 0)
return complex<FLOAT> (1.0);
complex<FLOAT> r (1.0);
complex<FLOAT> x (xin);
if (y < 0)
{
y = -y;
x = 1/x;
}
for (;;)
{
if (y & 1)
r *= x;
if (y >>= 1)
x *= x;
else
return r;
}
}
template <class FLOAT> complex<FLOAT>
sqrt (const complex<FLOAT>& x)
{
FLOAT r = abs (x);
FLOAT nr, ni;
if (r == 0.0)
nr = ni = r;
else if (real (x) > 0)
{
nr = sqrt (0.5 * (r + real (x)));
ni = imag (x) / nr / 2;
}
else
{
ni = sqrt (0.5 * (r - real (x)));
if (imag (x) < 0)
ni = - ni;
nr = imag (x) / ni / 2;
}
return complex<FLOAT> (nr, ni);
}
}