#include <math.h>
#include "libgfortran.h"
GFC_REAL_8
cabs (GFC_COMPLEX_8 z)
{
return hypot (REALPART (z), IMAGPART (z));
}
GFC_REAL_8
carg (GFC_COMPLEX_8 z)
{
GFC_REAL_8 arg;
arg = atan2 (IMAGPART (z), REALPART (z));
if (arg < 0)
return arg + 2 * M_PI;
else
return arg;
}
GFC_COMPLEX_8
cexp (GFC_COMPLEX_8 z)
{
GFC_REAL_8 a;
GFC_REAL_8 b;
GFC_COMPLEX_8 v;
a = REALPART (z);
b = IMAGPART (z);
COMPLEX_ASSIGN (v, cos (b), sin (b));
return exp (a) * v;
}
GFC_COMPLEX_8
clog (GFC_COMPLEX_8 z)
{
GFC_COMPLEX_8 v;
COMPLEX_ASSIGN (v, log (cabs (z)), carg (z));
return v;
}
GFC_COMPLEX_8
clog10 (GFC_COMPLEX_8 z)
{
GFC_COMPLEX_8 v;
COMPLEX_ASSIGN (v, log10 (cabs (z)), carg (z));
return v;
}
GFC_COMPLEX_8
cpow (GFC_COMPLEX_8 base, GFC_COMPLEX_8 power)
{
return cexp (power * clog (base));
}
GFC_COMPLEX_8
csqrt (GFC_COMPLEX_8 z)
{
GFC_REAL_8 re;
GFC_REAL_8 im;
GFC_COMPLEX_8 v;
re = REALPART (z);
im = IMAGPART (z);
if (im == 0.0)
{
if (re < 0.0)
{
COMPLEX_ASSIGN (v, 0.0, copysign (sqrt (-re), im));
}
else
{
COMPLEX_ASSIGN (v, fabs (sqrt (re)),
copysign (0.0, im));
}
}
else if (re == 0.0)
{
GFC_REAL_8 r;
r = sqrt (0.5 * fabs (im));
COMPLEX_ASSIGN (v, copysign (r, im), r);
}
else
{
GFC_REAL_8 d, r, s;
d = hypot (re, im);
if (re > 0)
{
r = sqrt (0.5 * d + 0.5 * re);
s = (0.5 * im) / r;
}
else
{
s = sqrt (0.5 * d - 0.5 * re);
r = fabs ((0.5 * im) / s);
}
COMPLEX_ASSIGN (v, r, copysign (s, im));
}
return v;
}