#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.

// 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); return 0;}

/*
  0.00000000000000000e+00   2.00000000000000000e+00 2i

  2.00000000000000000e+00   3.00000000000000000e+00 x
  1.67414922803554012e+00   8.95977476129838135e-01 r
  0.00000000000000000e+00   0.00000000000000000e+00 r*r-x

 -1.00000000000000000e+00  -1.00000000000000002e-03 x
  4.99999937500027257e-04  -1.00000012499996105e+00 r
 -2.22044604925031308e-16   0.00000000000000000e+00 r*r-x

 -1.00000000000000000e+00   0.00000000000000000e+00 x
  0.00000000000000000e+00   1.00000000000000000e+00 r
  0.00000000000000000e+00   0.00000000000000000e+00 r*r-x
*/