#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
// remove CMPLX definition in C11
static c CMPLX(d r, d i){c x; __real__ x = r; __imag__ x = i; return x;}

c csqrt(c x){
d A = creal(x), B = cimag(x);
d aA = fabs(A), aB = fabs(B);
if(A < 0) {c r = csqrt(-x); if(B<0) r = -r; return CMPLX(-cimag(r), creal(r));}
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 CMPLX(r, (B/(2*r)));}}
// We assume that underflow for multipliction yields 0.

/*
// unit test:
#include <stdio.h>
void pc(c x, char* m){printf("%25.17e %25.17e %s\n", __real__ x, __imag__ x, m);}
void tr(c x){c r = csqrt(x); printf("\n");
   pc(x, "x"); pc(r, "r"); pc(r*r-x, "r*r-x");}
int main(){pc(2i, "2i"); tr(2+3i); tr(-1-0.001i); tr(-1); 
  tr(0.003 - 4i); tr(-0.003 - 4i); tr(1e122-2e122i);
  tr(-0.7 + .002i); tr(-0.7 - .002i); tr(3e-122-2e-122i); return 0;}
*/