#include <assert.h>
#include <stdio.h>
#include <errno.h>

typedef struct {
	double real;
	double imag;
} Py_complex;

Py_complex
c_quot(Py_complex a, Py_complex b)
{
	/******************************************************************
	This was the original algorithm.  It's grossly prone to spurious
	overflow and underflow errors.  It also merrily divides by 0 despite
	checking for that(!).  The code still serves a doc purpose here, as
	the algorithm following is a simple by-cases transformation of this
	one:

	Py_complex r;
	double d = b.real*b.real + b.imag*b.imag;
	if (d == 0.)
		errno = EDOM;
	r.real = (a.real*b.real + a.imag*b.imag)/d;
	r.imag = (a.imag*b.real - a.real*b.imag)/d;
	return r;
	******************************************************************/

	/* This algorithm is better, and is pretty obvious:  first divide the
	 * numerators and denominator by whichever of {b.real, b.imag} has
	 * larger magnitude.  The earliest reference I found was to CACM
	 * Algorithm 116 (Complex Division, Robert L. Smith, Stanford
	 * University).  As usual, though, we're still ignoring all IEEE
	 * endcases.
	 */
	 Py_complex r;	/* the result */
 	 const double abs_breal = b.real < 0 ? -b.real : b.real;
	 const double abs_bimag = b.imag < 0 ? -b.imag : b.imag;

	 if (abs_breal >= abs_bimag) {
 		/* divide tops and bottom by b.real */
	 	if (abs_breal == 0.0) {
	 		errno = EDOM;
	 		r.real = r.imag = 0.0;
	 	}
	 	else {
	 		const double ratio = b.imag / b.real;
	 		const double denom = b.real + b.imag * ratio;
	 		r.real = (a.real + a.imag * ratio) / denom;
	 		r.imag = (a.imag - a.real * ratio) / denom;
	 	}
	}
	else {
		/* divide tops and bottom by b.imag */
		const double ratio = b.real / b.imag;
		const double denom = b.real * ratio + b.imag;
		assert(b.imag != 0.0);
		r.real = (a.real * ratio + a.imag) / denom;
		r.imag = (a.imag * ratio - a.real) / denom;
	}
	return r;
}

int main(int argc, char* argv[])
{
	Py_complex x = {4.6051701859880918,  0.0}; /* log(complex({100.0, 0.0}) */
	Py_complex y = {0.69314718055994529, 0.0}; /* log(complex({  2.0, 0.0}) */
	Py_complex q = c_quot(x, y);
	/* expect: 6.643856 + j 0.000000 */
	printf("%f + j %f\n", q.real, q.imag);
}

/* compile and run as:

> rm a.out ; cc -xO0 -xlibmieee c_quot.c -lm ; ./a.out
6.643856 + j 0.000000
> rm a.out ; cc -xO5 -xlibmieee c_quot.c -lm ; ./a.out
6.643856 + j 0.000000
> rm a.out ; cc -xO0 c_quot.c -lm ; ./a.out           
6.643856 + j 0.000000
> rm a.out ; cc -xO5 c_quot.c -lm ; ./a.out
6.643856 + j 0.000000

> rm a.out ; cc -xtarget=native64 -xO0 -xlibmieee c_quot.c -lm ; ./a.out
6.643856 + j 0.000000
> rm a.out ; cc -xtarget=native64 -xO5 -xlibmieee c_quot.c -lm ; ./a.out
6.643856 + j 0.000000
> rm a.out ; cc -xtarget=native64 -xO0 c_quot.c -lm ; ./a.out           
6.643856 + j 0.000000
> rm a.out ; cc -xtarget=native64 -xO5 c_quot.c -lm ; ./a.out
6.643856 + j 0.000000

*/