// filename:c2011-G-5-1-ex.c
// original examples and/or notes:
// 		(c) ISO/IEC JTC1 SC22 WG14 N1570, April 12, 2011
// 			C2011 G.5.1 Multiplicative operators
// compile and output mecG.5.1 Multiplicative operatorshanism:
// 		(c) Ogawa Kiyoshi, kaizen@gifu-u.ac.jp, December.29, 2013
// compile errors and/or wornings:
// 		(c) Apple LLVM version 4.2 (clang-425.0.27) (based on LLVM 3.2svn)
// 			Target: x86_64-apple-darwin11.4.2 //Thread model: posix
// 		(c) LLVM 2003-2009 University of Illinois at Urbana-Champaign.
// gcc-4.9 (GCC) 4.9.0 20131229 (experimental)
// Copyright (C) 2013 Free Software Foundation, Inc.
#include <stdio.h>

// Example 1
#include <math.h>
#include <complex.h>
/* Multiplyz * w... */
double complex _Cmultd(double complex z, double complex w)
{
#pragma STDC FP_CONTRACT OFF
double a, b, c, d, ac, bd, ad, bc, x, y;
a = creal(z); b = cimag(z);
c = creal(w); d = cimag(w);
ac = a * c; bd = b * d;
ad = a * d; bc = b * c;
x = ac - bd; y = ad + bc;
if (isnan(x) && isnan(y)) {
/* Recover infinities that computed as NaN+iNaN ... */
int recalc = 0;
if (isinf(a) || isinf(b)) { // z is infinite
/* "Box" the infinity and change NaNs in the other factor to 0 */
a = copysign(isinf(a) ? 1.0 : 0.0, a);
b = copysign(isinf(b) ? 1.0 : 0.0, b);
if (isnan(c)) c = copysign(0.0, c);
		if (isnan(d)) d = copysign(0.0, d);
		recalc = 1;
}
if (isinf(c) || isinf(d)) { // w is infinite
/* "Box" the infinity and change NaNs in the other factor to 0 */
c = copysign(isinf(c) ? 1.0 : 0.0, c);
d = copysign(isinf(d) ? 1.0 : 0.0, d);
if (isnan(a)) a = copysign(0.0, a);
if (isnan(b)) b = copysign(0.0, b);
recalc = 1;
}
if (!recalc && (isinf(ac) || isinf(bd) ||
isinf(ad) || isinf(bc))) {
/* Recover infinities from overflow by changing NaNs to 0 ... */
if (isnan(a)) a = copysign(0.0, a);
if (isnan(b)) b = copysign(0.0, b);
if (isnan(c)) c = copysign(0.0, c);
if (isnan(d)) d = copysign(0.0, d);
recalc = 1;
}
if (recalc) {
x = INFINITY * ( a * c - b * d );
y = INFINITY * ( a * d + b * c );
}
}
return x + I * y;
}		

// Example 2		
/* Dividez / w ... */
double complex _Cdivd(double complex z, double complex w)
{
#pragma STDC FP_CONTRACT OFF
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(void)
{
	double complex z, w;
	double complex d = _Cmultd( z, w);
	double complex e =_Cdivd( z,  w);
return printf("G.5.1 Multiplicative operators %f %f\n", d, e);
}
// output may be 
//G.5.1 Multiplicative operators 0.000000 0.000000