clang is almost useless for complex arithmetic

[Steve Kargl]

It appears that clang developers have chosen the naive
complex division algorithm, and it does not matter whether
one turns CX_LIMITED_RANGE on or off.  This means that 
if one uses clang with complex types, one must be careful
with the range of values allowed in complex division.  In
other words, implementation of complex libm routines cannot
use complex data types and must fallback to a decomposition
into real and imaginary components. 

For the record, I am using

% cc --version
FreeBSD clang version 3.4 (tags/RELEASE_34/final 197956) 20140216
Target: x86_64-unknown-freebsd11.0

With the program that follows my sig, if I do

% cc -o z -O -DOFF cmplx.c -lm && ./z > pragma_off
% cc -o z -O cmplx.c -lm && ./z > pragma_on
% diff -u pragma_off pragma_on | wc -l
0
% cat pragma_on

Expected: 0x1p-1023 -0x1p-1023
Compiler: 0x0p+0 -0x0p+0
Improved: 0x1p-1023 -0x1p-1023

Expected: 0x1p+1023 0x0p+0
Compiler: inf nan
Improved: 0x1p+1023 0x0p+0

Expected: 0x1p+346 -0x1p-1008
Compiler: nan -0x0p+0
Improved: 0x1p+346 -0x1p-1008

Expected: 0x1p+1023 0x0p+0
Compiler: inf 0x0p+0
Improved: 0x1p+1023 0x0p+0

Expected: 0x1p+364 -0x1p-1072
Compiler: nan -0x0p+0
Improved: 0x1p+364 -0x1p-1072

Expected: 0x1p-1072 0x1p+20
Compiler: 0x0p+0 nan
Improved: 0x1p-1072 0x1p+20

Expected: 0x1.ffffffffff8p+961 0x1.ffffffffff8p+982
Compiler: nan nan
Improved: 0x1.ffffffffff8p+961 0x1.ffffffffff8p+982

Expected: 0x1.3333333333333p-1 0x1.999999999999ap-3
Compiler: nan nan
Improved: 0x1.3333333333333p-1 0x1.999999999999ap-3

Expected: 0x1p-9 -0x1p-9
Compiler: nan nan
Improved: 0x1p-9 -0x1p-9

Expected: 0x1p-279 0x1.f8p-729
Compiler: 0x1p-279 0x0p+0
Improved: 0x1p-279 0x1.f8p-729

Here, "Expected" is the expected result.  "Compiler" is what clang
produces. "Improved" implements the algorithm suggested in 
Michael Baudin, Robert L. Smith, "A Robust Complex Division in Scilab,"
http://arxiv.org/abs/1210.4539

-- 
Steve

#include <complex.h>
#include <float.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

#ifdef OFF
#pragma CX_LIMITED_RANGE off
#else
#pragma CX_LIMITED_RANGE on
#endif

typedef union {
        double complex f;
        double a[2];
} double_complex;

#define REALPART(z)     ((z).a[0])
#define IMAGPART(z)     ((z).a[1])

static __inline double complex
cpack(double x, double y)
{
        double_complex z;

        REALPART(z) = x;
        IMAGPART(z) = y;
        return (z.f);
}

/*
 * For details of internal() and divi() see:
 * Michael Baudin, Robert L. Smith, "A Robust Complex Division in Scilab,"
 * http://arxiv.org/abs/1210.4539
 */
double complex
internal(double a, double b, double c, double d, int l)
{

	static const double ti = 0x1.p-1022, rti = 0x1.p1022, hu = 0x1.p1023,
	    rhu = 0x1.p-1023;

	double e, f, r, t, ma, mi;

	r = d / c;

	if (r != 0) {
		if (fabs(c) < rhu) {
			t = c * hu + d * hu * r;
			e = (a * hu + b * hu * r) / t;
			f = (b * hu - a * hu * r) / t;
		} else if (fabs(c) > rti) {
			t = 1 / (c * ti + d * ti * r);
			e = (a * ti + b * ti * r) * t;
			f = (b * ti - a * ti * r) * t;
		} else {
			t = 1 / (c + d * r);
			ma = fmax(fabs(a), fabs(b));
			mi = fmin(fabs(a), fabs(b));
			if (ma < rti && mi * r > rhu) {
				e = (a + b * r) * t;
				f = (b - a * r) * t;
			} else {
				r = r * t;
				e = a * t + b * r;
				f = b * t - a * r;
			}
		}
	} else {
		t = 1 / c;
		e = (a + d * (b / c)) * t;
		f = (b - d * (a / c)) * t;
	}
    
	if (l == 1) f = -f;

   	return (cpack(e, f));

}

double complex
divi(double complex z1, double complex z2)
{
	double complex r;
	double a, b, c, d;

	a = creal(z1);
	b = cimag(z1);
	c = creal(z2);
	d = cimag(z2);

	if (fabs(d) <= fabs(c))
		r = internal(a, b, c, d, 0);
	else
		r = internal(b, a, d, c, 1);

	return (r); 
}

int
main(void) 
{
	/* Numerators. */
	double complex x[10] = {
		cpack(1, 1), cpack(1, 1), cpack(0x1.p1023, 0x1.p-1023),
		cpack(0x1.p1023, 0x1.p1023), cpack(0x1.p1020, 0x1.p-844),
		cpack(0x1.p-71, 0x1.p1021), cpack(0x1.p-347, 0x1.p-54),
		cpack(0x1.p-1074, 0x1.p-1074), cpack(0x1.p1015, 0x1.p-989),
		cpack(0x1.p-622, 0x1.p-1071)
	};

	/* Denominators. */
	double complex y[10] = {
		cpack(1, 0x1.p1023), cpack(0x1.p-1023, 0x1.p-1023),
    		cpack(0x1.p677, 0x1.p-677), cpack(1, 1),
		cpack(0x1.p656, 0x1.p-780), cpack(0x1.p1001, 0x1.p-323),
		cpack(0x1.p-1037, 0x1.p-1058), cpack(0x1.p-1073, 0x1.p-1074),
		cpack(0x1.p1023, 0x1.p1023), cpack(0x1.p-343, 0x1.p-798)
	};

	double complex r[10] = {
		cpack(0x1.p-1023, -0x1.p-1023), cpack(0x1.p1023, 0),
		cpack(0x1.p346, -0x1.p-1008), cpack(0x1.p1023, 0),
		cpack(0x1.p364, -0x1.p-1072), cpack(0x1.p-1072, 0x1.p20),
		cpack(3.898125604559113300E289, 8.174961907852353577E295),
		cpack(0.6, 0.2),
		cpack(0.001953125, -0.001953125),
		cpack(1.02951151789360578E-84, 6.97145987515076231E-220)
	};

	double complex a, z[10];
	int j;

	for (j = 0; j < 10; j++) {
		z[j] = x[j] / y[j];
		a = divi(x[j], y[j]);
		printf("Expected: %la %la\n", creal(r[j]), cimag(r[j]));
		printf("Compiler: %la %la\n", creal(z[j]), cimag(z[j]));
		printf("Improved: %la %la\n", creal(a),    cimag(a));
		printf("\n");
	}

	return 0;
}