#pragma STDC CX_LIMITED_RANGE OFF
#pragma STDC FP_CONTRACT OFF
#include <math.h>
#include <stdio.h>
#include <complex.h>

#ifdef BROKEN_C99
#define creal(x) ((double)(x))
#define cimag(x) ((double)(-I*(x)))
#define INFINITY (1.0/0.0)
#define isfinite(x) (! isinf(x) && ! isnan(x))
extern double fmax(double,double);
#endif



double complex cdivd99(double complex z, double complex w)
{
    double a, b, c, d, logbw, denom, x, y;
    int ilogbw = 0;
    a = creal(z); b = cimag(z);
    c = creal(w); d = cimag(w);
    logbw = logb(fmax(fabs(c), fabs(d)));
    if (isfinite(logbw)) {
        ilogbw = (int)logbw;
        c = scalbn(c, -ilogbw); d = scalbn(d, -ilogbw);
    }
    denom = c * c + d * d;
    x = scalbn((a * c + b * d) / denom, -ilogbw);
    y = scalbn((b * c - a * d) / denom, -ilogbw);
    if (isnan(x) && isnan(y)) {
        if ((denom == 0.0) &&
            (!isnan(a) || !isnan(b))) {
            x = copysign(INFINITY, c) * a;
            y = copysign(INFINITY, c) * b;
        }
        else if ((isinf(a) || isinf(b)) &&
            isfinite(c) && isfinite(d)) {
            a = copysign(isinf(a) ? 1.0 : 0.0, a);
            b = copysign(isinf(b) ? 1.0 : 0.0, b);
            x = INFINITY * ( a * c + b * d );
            y = INFINITY * ( b * c - a * d );
        }
        else if (isinf(logbw) &&
            isfinite(a) && isfinite(b)) {
            c = copysign(isinf(c) ? 1.0 : 0.0, c);
            d = copysign(isinf(d) ? 1.0 : 0.0, d);
            x = 0.0 * ( a * c + b * d );
            y = 0.0 * ( b * c - a * d );
        }
    }
    return x + I * y;
}



double complex cdivd11 (double complex z, double complex w) {
    double a, b, c, d, logbw, denom, x, y;
    int ilogbw = 0;

    a = creal(z); b = cimag(z);
    c = creal(w); d = cimag(w);
    logbw = logb(fmax(fabs(c), fabs(d)));
    if (isfinite(logbw)) {
        ilogbw = (int)logbw;
        c = scalbn(c, -ilogbw); d = scalbn(d, -ilogbw);
    }
    denom = c * c + d * d;
    x = scalbn((a * c + b * d) / denom, -ilogbw);
    y = scalbn((b * c - a * d) / denom, -ilogbw);
    /* Recover infinities and zeros that computed as NaN+iNaN;
        the only cases are nonzero/zero, infinite/finite, and
        finite/infinite, ... */
    if (isnan(x) && isnan(y)) {
        if ((denom == 0.0) && (!isnan(a) || !isnan(b))) {
            x = copysign(INFINITY, c) * a;
            y = copysign(INFINITY, c) * b;
        } else if ((isinf(a) || isinf(b)) && isfinite(c) && isfinite(d)) {
            a = copysign(isinf(a) ? 1.0 : 0.0, a);
            b = copysign(isinf(b) ? 1.0 : 0.0, b);
            x = INFINITY * ( a * c + b * d );
            y = INFINITY * ( b * c - a * d );
        } else if ((logbw == INFINITY) && isfinite(a) && isfinite(b)) {
            c = copysign(isinf(c) ? 1.0 : 0.0, c);
           d = copysign(isinf(d) ? 1.0 : 0.0, d);
            x = 0.0 * ( a * c + b * d );
            y = 0.0 * ( b * c - a * d );
        }
    }
    return x + I * y;
}



int main() {
    int i;
    double d, r, z;
    double complex c1, c2, c3;

    for (i = 0; i < 20; ++i) {
        d = i*0.1e308;
        c1 = (1.0e308+I*1.0e308)/(1.0e308+I*d);
        c2 = cdivd99(1.0e308+I*1.0e308,1.0e308+I*d);
        c3 = cdivd11(1.0e308+I*1.0e308,1.0e308+I*d);
        if (1.0e308 > d) {
            r = d/1.0e308;
            z = 1.0e308+d*r;
            printf("%.2e (%.3f,%.3f) (%.3f,%.3f) (%.3f,%.3f) (%.3f,%.3f)\n",
                d,creal(c1),cimag(c1),creal(c2),cimag(c2),
                creal(c3),cimag(c3),
                (1.0e308+1.0e308*r)/z,(1.0e308-1.0e308*r)/z);
        } else {
            r = 1.0e308/d;
            z = d+1.0e308*r;
            printf("%.2e (%.3f,%.3f) (%.3f,%.3f) (%.3f,%.3f) (%.3f,%.3f)\n",
                d,creal(c1),cimag(c1),creal(c2),cimag(c2),
                creal(c3),cimag(c3),
                (1.0e308*r+1.0e308)/z,(1.0e308*r-1.0e308)/z);
        }
    }
    return 0;
}