#include <math.h>
typedef double d;
#include <complex.h>
typedef double _Complex c;
#define bg2 0x1.p300
#define bg  0x1.p150
#define sm2 0x1.p-300
#define sm  0x1.p-150
c csqrt(c x){
d A = creal(x), B = cimag(x);
d aA = fabs(A), aB = fabs(B);
if(A < 0) return (B<0?-1.i:1.i)*csqrt(-x);
if(aA > bg2 || aB > bg2) return
    (aA==INFINITY || aB==INFINITY) ? x :  bg*csqrt(sm2*x);
if(aA < sm2 && aB < sm2) return (aA==0 && aB==0) ? x : sm*csqrt(bg2*x);
{d r = sqrt((A + sqrt(A*A + B*B))/2); return r + 1i*(B/(2*r));}}
// We assume that underflow for multipliction yields 0.