#include <math.h>
#include <stdio.h> 
typedef long double ld;
#include <complex.h>
#if 0
typedef long double br;
typedef long double _Complex c;
#define sqrtv sqrtl
#define bg2 0x1.p8000L
#define bg  0x1.p4000L
#define sm2 0x1.p-8000L
#define sm  0x1.p-4000L
#define fabsv fabsl
#else
typedef double br;
typedef double _Complex c;
#define sqrtv sqrt
#define bg2 0x1.p300
#define bg  0x1.p150
#define sm2 0x1.p-300
#define sm  0x1.p-150
#define fabsv fabs
#endif
void px(ld r){
  typedef union{ld v; struct{int a; int b; int C; int d;} fs;} ww;
  ww x; x.v = r;
  printf("%Le %04x %08x%08x; ", r, 0xffff&x.fs.C, x.fs.b, x.fs.a);}
void pC(c x, char * s){
   px(__real__ x); px(__imag__ x); printf(" %s\n", s);}
c pv(c cr){
if(1) pC(cr, "out");
 return cr;}
br prv(br cr){if(1) px(cr); return cr;}
c Sqrt(c x){
pC(x, "in");
br A = __real__ x, B = __imag__ x;
br aA = fabsv(A), aB = fabsv(B);
if(A < 0) return pv((B<0?-1.i:1.i)*Sqrt(-x));
if(aA > bg2 || aB > bg2)
   {if(aA==INFINITY || aB==INFINITY) return pv(x);
    else return pv(bg*Sqrt(sm2*x));}
if(aA < sm2 && aB < sm2)
   {if(aA==0 && aB==0) return pv(x); else return pv(sm*Sqrt(bg2*x));}
{br r = sqrtv((A + sqrtv(A*A + B*B))/2); return pv(r + 1i*(B/(2*r)));}}
// We assume that underflow for multipliction yields 0.

void pc(c x){br mx = fmaxl(fabsl(__real__ x), fabsl(__imag__ x));
  c y = Sqrt(x), z = x==0?0:(y*y - x)*(1./mx);
printf ("%9Le %9Le rel err\n", (ld) __real__ z, (ld) __imag__ z);}
int main(){
  if(1) {int j = 14; c r = 2; while(j--) {r = r*r; pc (r);}}
  if(1) {pc(3); pc(1.e-132L + 2.e-132iL);}
  if(1) {pc(2i); pc(1. + 2.i); pc(1.e201 + 2.e201i); pc(-3.e-155 + 4.e-154i);}
  if(1) {pc(0x1.p-16382L + 2.i); pc(0); pc(0x1.p-16382L);}
return 0;}