#include <math.h>
typedef long double ld;
#include <complex.h>
typedef long double _Complex c;
// The following 4 values are good for the X86 extended precision.
#define bg2 0x1.p8000L
#define bg  0x1.p4000L
#define sm2 0x1.p-8000L
#define sm  0x1.p-4000L
c csqrtl(c x){
ld A = __real__ x, B = __imag__ x;
ld aA = fabsl(A), aB = fabsl(B);
if(A < 0) return (B<0?-1.i:1.i)*csqrtl(-x);
if(aA > bg2 || aB > bg2) return
    (aA==INFINITY || aB==INFINITY) ? x :  bg*csqrtl(sm2*x);
if(aA < sm2 && aB < sm2) return (aA==0 && aB==0) ? x : sm*csqrtl(bg2*x);
{ld r = sqrtl((A + sqrtl(A*A + B*B))/2); return r + 1i*(B/(2*r));}}
// We assume that underflow for multipliction yields 0.