JTC1/SC22/WG14
N1399
Complex Multiply and Complex Divide, taking into account IEEE-754 (IEC
60559) signed zeros, signed infinities, NaN, and C99 _Imaginary_I.
by
Fred J. Tydeman
WG14 N1399
2009/09/25
* Problems with C99 Annex G.5.1 examples
Yuri Akutin of Intel (in 2003) discovered that C99 Annex G.5.1 _Cdivd
of (1.0+1.0*I)/(0.0+NAN*I) produced (0.0+0.0*I) instead of the
expected (NAN+NAN*I). This can be fixed by changing "isinf(logbw)" to
"(INFINITY==logbw)". This was brought to the attention of Jim Thomas
and Fred Tydeman (the two main authors of Annex G IEC 60559-compatible
complex arithmetic in C99). Jim and Fred each then created initial
reference implementations of cmult (complex multiply) and cdiv
(complex divide). Using his reference implementations, Fred then
found other cases where Annex G.5.1 produces the "wrong" answer.
Here are specific operands that show the problems. They use the
following shorthands:
nan = NAN
inf = INFINITY
max = DBL_MAX
min = DBL_MIN
fin = any non-zero finite number
** complex multiply:
Operands as by C99 correct Comments
Need to add special case code for overflow-overflow that exactly cancel.
(max+max*I)*(max+max*I) (nan+inf*I) (0.0+inf*I) overflow-overflow
(max+max*I)*(max-max*I) (inf+nan*I) (inf+0.0*I) overflow-overflow
Need to add special case code for inf-overflow.
(inf+max*I)*(inf+max*I) (nan+inf*I) (inf+inf*I) inf-overflow
(inf+max*I)*(-max+inf*I) (inf+nan*I) (inf-inf*I) inf-overflow
(min+max*I)*(inf+max*I) (nan+inf*I) (inf+inf*I) inf-overflow
(min+max*I)*(max-inf*I) (inf+nan*I) (inf-inf*I) inf-overflow
Need to add special case code for nan with inf.
(nan+inf*I)*(inf+inf*I) (-inf+inf*I) (-inf+nan*I) NaN in => NaN out
(nan+inf*I)*(1.0+max*I) (-inf+inf*I) (-inf+nan*I) NaN in => NaN out
(inf+nan*I)*(1.0+1.0*I) (inf+inf*I) (inf+nan*I) NaN in => NaN out
Need to add special case code for nan with finite.
(nan-max*I)*(max+max*I) (inf-inf*I) (inf+nan*I) NaN in => NaN out
(nan+max*I)*(max+1.0*I) (-inf+inf*I) (nan+inf*I) NaN in => NaN out
(1.0+nan*I)*(inf+inf*I) (inf+inf*I) (inf+nan*I) NaN in => NaN out
Need to add special case code for either being on an axis.
The following are axis * zero. [3 zeros]
(inf+0.0*I)*(0.0+0.0*I) (nan+nan*I) (nan+0.0*I) x-axis*zero
(0.0+nan*I)*(0.0+0.0*I) (nan+nan*I) (0.0+nan*I) y-axis*zero
(0.0+0.0*I)*(nan+0.0*I) (nan+nan*I) (nan+0.0*I) zero*x-axis
(0.0+0.0*I)*(0.0+inf*I) (nan+nan*I) (0.0+nan*I) zero*y-axis
The following are non-axis * zero. [2 zeros]
(inf+1.0*I)*(0.0+0.0*I) (nan+nan*I) (nan+0.0*I) non-axis*zero
(max+nan*I)*(0.0+0.0*I) (nan+nan*I) (0.0+nan*I) non-axis*zero
(0.0+0.0*I)*(nan+max*I) (nan+nan*I) (nan+0.0*I) zero*non-axis
(0.0+0.0*I)*(-1.+inf*I) (nan+nan*I) (-0.+nan*I) zero*non-axis
The following are axis * axis. [2 zeros]
(0.0+max*I)*(0.0+inf*I) (-inf+nan*I) (-inf+0.0*I) y-axis*y-axis
(0.0-1.0*I)*(0.0+inf*I) (inf+nan*I) (inf+0.0*I) y-axis*y-axis
(0.0+1.0*I)*(inf+0.0*I) (nan+inf*I) (0.0+inf*I) y-axis*x-axis
(inf+0.0*I)*(0.0+1.0*I) (nan+inf*I) (0.0+inf*I) x-axis*y-axis
(nan+0.0*I)*(0.0+1.0*I) (nan+nan*I) (0.0+nan*I) x-axis*y-axis
(nan+0.0*I)*(1.0+0.0*I) (nan+nan*I) (nan+0.0*I) x-axis*x-axis
The following are non-axis * axis. [1 zero]
(0.0+1.0*I)*(1.0+inf*I) (-inf+nan*I) (-inf+1.0*I) y-axis*non-axis
(0.0+1.0*I)*(inf+1.0*I) (nan+inf*I) (-1.+inf*I) y-axis*non-axis
(max+0.0*I)*(inf+max*I) (inf+nan*I) (inf+inf*I) x-axis*non-axis
(nan+1.0*I)*(0.0+1.0*I) (nan+nan*I) (-1.+nan*I) non-axis*y-axis
(inf+1.0*I)*(0.0+inf*I) (nan+inf*I) (-inf+inf*I) non-axis*y-axis
(inf+1.0*I)*(inf+0.0*I) (inf+nan*I) (inf+inf*I) non-axis*x-axis
(nan+1.0*I)*(1.0+0.0*I) (nan+nan*I) (nan+1.0*I) non-axis*x-axis
** complex divide:
Operands as by C99 correct Comments
The (nan+0*I) is done as if (inf+0*I) due to isinf(logbw);
should be INFINITY==logbw.
(1.0+0.0*I)/(nan+0.0*I) (0.0+0.0*I) (nan+0.0*I) x-axis / x-axis
(1.0+0.0*I)/(0.0+nan*I) (0.0+0.0*I) (0.0+nan*I) x-axis / y-axis
(0.0+1.0*I)/(nan+0.0*I) (0.0+0.0*I) (0.0+nan*I) y-axis / x-axis
(0.0+1.0*I)/(0.0+nan*I) (0.0+0.0*I) (nan+0.0*I) y-axis / y-axis
(1.0+1.0*I)/(nan+0.0*I) (0.0+0.0*I) (nan+nan*I) any / x-axis
(1.0+1.0*I)/(0.0+nan*I) (0.0+0.0*I) (nan+nan*I) any / y-axis
Need to add special case code for (0+0*I)/(finite+nan*I) and (0+0*I)/(nan+finite*I)
(0.0+0.0*I)/(1.0+nan*I) (nan+nan*I) (0.0+0.0*I) zero / non-zero
(0.0+0.0*I)/(nan+1.0*I) (nan+nan*I) (0.0+0.0*I) zero / non-zero
Need to add special case code for (inf+nan*I)/finite or (nan+inf*I)/finite.
(nan+inf*I)/(1.0+1.0*I) (inf+inf*I) (nan+inf*I) NaN in => NaN out
(nan+inf*I)/(1.0+max*I) (inf+inf*I) (nan+inf*I) NaN in => NaN out
(inf+nan*I)/(1.0+1.0*I) (inf-inf*I) (inf+nan*I) NaN in => NaN out
(inf+nan*I)/(1.0+max*I) (inf-inf*I) (inf+nan*I) NaN in => NaN out
These are due to isnan(x) && isnan(y);
should be isnan(x) || isnan(y).
(inf+inf*I)/(1.0+max*I) (inf+nan*I) (inf-inf*I) pure inf/finite => no NaN
(inf+inf*I)/(-1.+max*I) (nan-inf*I) (inf-inf*I) pure inf/finite => no NaN
(inf+inf*I)/(max+1.0*I) (inf+nan*I) (inf+inf*I) pure inf/finite => no NaN
(inf+inf*I)/(max-1.0*I) (nan+inf*I) (inf+inf*I) pure inf/finite => no NaN
The (min+max*I) is done as if (0+max*I) due to limited exponent range;
make sure have non-zero.
(0.0+inf*I)/(min+max*I) (inf+nan*I) (inf+inf*I) inf/finite => no NaN
(min+inf*I)/(min+max*I) (inf+nan*I) (inf+inf*I) inf/finite => no NaN
(max+inf*I)/(max+min*I) (nan+inf*I) (inf+inf*I) inf/finite => no NaN
(inf+0.0*I)/(max+min*I) (inf+nan*I) (inf-inf*I) inf/finite => no NaN
(inf+1.0*I)/(min+max*I) (nan-inf*I) (inf-inf*I) inf/finite => no NaN
Need to add special case code for (nan with large)/small
(max+nan*I)/(min+min*I) (nan+nan*I) (inf+nan*I) large/small => overflow to inf
(nan+max*I)/(min+min*I) (nan+nan*I) (nan+inf*I) large/small => overflow to inf
Need to add special case code for small/(nan w/ large)
(min+min*I)/(nan+max*I) (nan+nan*I) (0.0-0.0*I) small/large => underflow to 0
(min+min*I)/(max+nan*I) (nan+nan*I) (0.0+0.0*I) small/large => underflow to 0
Need to add special case code for all the following (vector(s) on an axis).
The following is zero / zero. [4 zeros]
(0.0+0.0*I)/(0.0+0.0*I) (nan+nan*I) (nan+0.0*I) zero/zero on x-axis
The following are axis / zero. [3 zeros]
(1.0+0.0*I)/(0.0+0.0*I) (inf+nan*I) (inf+0.0*I) x-axis/zero
(nan+0.0*I)/(0.0+0.0*I) (nan+nan*I) (nan+0.0*I) x-axis/zero
(0.0+1.0*I)/(0.0+0.0*I) (nan+inf*I) (0.0+inf*I) y-axis/zero
(0.0+nan*I)/(0.0+0.0*I) (nan+nan*I) (0.0+nan*I) y-axis/zero
The following are axis / axis. [2 zeros]
(inf+0.0*I)/(inf+0.0*I) (nan+nan*I) (nan+0.0*I) x-axis / x-axis
(inf+0.0*I)/(1.0+0.0*I) (inf+nan*I) (inf+0.0*I) x-axis / x-axis
(nan+0.0*I)/(nan+0.0*I) (nan+nan*I) (nan+0.0*I) x-axis / x-axis
(inf+0.0*I)/(nan+0.0*I) (nan+nan*I) (nan+0.0*I) x-axis / x-axis
(inf+0.0*I)/(0.0+inf*I) (nan+nan*I) (0.0+nan*I) x-axis / y-axis
(inf+0.0*I)/(0.0+1.0*I) (nan-inf*I) (0.0-inf*I) x-axis / y-axis
(nan+0.0*I)/(0.0+nan*I) (nan+nan*I) (0.0+nan*I) x-axis / y-axis
(inf+0.0*I)/(0.0+nan*I) (nan+nan*I) (0.0+nan*I) x-axis / y-axis
(0.0+inf*I)/(inf+0.0*I) (nan+nan*I) (0.0+nan*I) y-axis / x-axis
(0.0+inf*I)/(1.0+0.0*I) (nan+inf*I) (0.0+inf*I) y-axis / x-axis
(0.0+nan*I)/(nan+0.0*I) (nan+nan*I) (0.0+nan*I) y-axis / x-axis
(0.0+nan*I)/(inf+0.0*I) (nan+nan*I) (0.0+nan*I) y-axis / x-axis
(0.0+inf*I)/(0.0+inf*I) (nan+nan*I) (nan+0.0*I) y-axis / y-axis
(0.0+inf*I)/(0.0+1.0*I) (inf+nan*I) (inf+0.0*I) y-axis / y-axis
(0.0+nan*I)/(0.0+nan*I) (nan+nan*I) (nan+0.0*I) y-axis / y-axis
(0.0+nan*I)/(0.0+inf*I) (nan+nan*I) (nan+0.0*I) y-axis / y-axis
The following are non-axis / axis. [1 zero]
(1.0+inf*I)/(inf+0.0*I) (nan+nan*I) (0.0+nan*I) non-axis/x-axis
(inf+1.0*I)/(inf+0.0*I) (nan+nan*I) (nan+0.0*I) non-axis/x-axis
(inf+1.0*I)/(1.0+0.0*I) (inf+nan*I) (inf+1.0*I) non-axis/x-axis
(1.0+inf*I)/(1.0+0.0*I) (nan+inf*I) (1.0+inf*I) non-axis/x-axis
(1.0+nan*I)/(1.0+0.0*I) (nan+nan*I) (1.0+nan*I) non-axis/x-axis
(nan+1.0*I)/(1.0+0.0*I) (nan+nan*I) (nan+1.0*I) non-axis/x-axis
(1.0+nan*I)/(inf+0.0*I) (nan+nan*I) (0.0+nan*I) non-axis/x-axis
(nan+1.0*I)/(inf+0.0*I) (nan+nan*I) (nan+0.0*I) non-axis/x-axis
(inf+1.0*I)/(0.0+inf*I) (nan+nan*I) (0.0+nan*I) non-axis/y-axis
(1.0+inf*I)/(0.0+inf*I) (nan+nan*I) (nan-0.0*I) non-axis/y-axis
(1.0+inf*I)/(0.0+1.0*I) (inf+nan*I) (inf-1.0*I) non-axis/y-axis
(inf+1.0*I)/(0.0+1.0*I) (nan-inf*I) (1.0-inf*I) non-axis/y-axis
(nan+1.0*I)/(0.0+1.0*I) (nan+nan*I) (1.0+nan*I) non-axis/y-axis
(1.0+nan*I)/(0.0+1.0*I) (nan+nan*I) (nan-1.0*I) non-axis/y-axis
(nan+1.0*I)/(0.0+inf*I) (nan+nan*I) (0.0+nan*I) non-axis/y-axis
(1.0+nan*I)/(0.0+inf*I) (nan+nan*I) (nan-0.0*I) non-axis/y-axis
Proposal number 1: Remove G.5.1, Example 1 and 2.
* Different ways of viewing complex numbers
The rectangular (x,y) is equivalent to x+y*I, just a different
notation for this writeup. There are two other ways to look at a
complex number: polar and hybrid. Polar format consists of a length
and a direction angle. Hybrid format consists of a length and a
direction vector on the unit square. That is,
Rectangular (x,y), eg,
3+4*I => (3, 4)
Polar (hypot(x,y),atan2(y,x)), aka (cabs,carg), eg,
3+4*I => (5, 53.13deg) or (5, 0.927rad)
Hybrid len*(x/len,y/len), where len=max(fabs(x),fabs(y)), eg,
3+4*I => [4*(.75, 1.)]
It helps to consider both the rectangular and the polar viewpoints of
complex numbers when determining the best way to compute a result.
Using both, results in the hybrid form that uses the polar form for
the length and the rectangluar form for the direction.
z * w = (a + b*I)*(c + d*I) = (a*c-b*d) + (b*c+a*d)*I
polar: (cabsz,cargz)*(cabsw,cargw) => (cabsz*cabsw,cargz+cargw)
hybrid: [lenz*(a,b)]*[lenw*(c,d)] => lenz*lenw*((ac-bd) + (bc+ad)*I)
z a + b*I (a+b*I)(c-d*I) (a*c+b*d) + (b*c-a*d)*I
- = ------- = -------------- = -----------------------
w c + d*I (c+d*I)(c-d*I) c*c + d*d
(cabsz,cargz) cabsz
polar: ------------- => (-----,cargz-cargw)
(cabsw,cargw) cabsw
lenz*(a,b) lenz (a*c+b*d) + (b*c-a*d)*I
hybrid: ---------- = ---- * -----------------------
lenw*(c,d) lenw c*c + d*d
* Assumptions on classification of complex numbers
A complex number, a+I*b, has two components.
A real or imaginary component is classified into four ways:
NaN, infinity, zero, or finite,
(that is, finite in the following excludes zero).
If both components are (signed) zero, then it is a complex zero.
If either component is an INFINITY, then the complex number is
considered a complex infinity (even if the other component is a NaN).
If both components are a NaN, then the number is a complex NaN.
If one component is a NaN and the other component is zero, then the
number is considered a polar vector with length that could be zero.
So, this kind of NaN * complex infinity is NaN (as 0*INF is NaN).
If one component is a NaN and the other component is finite, then the
number is considered a polar vector with length as least as long as
the finite. So, this kind of NaN * complex infinity is complex
infinity (as finite*INF is INF; also (max,NaN)*(4,3) is (INF,NaN) as a
length of max times a length of 5 is overflow to INF).
In most cases, if any of the four input components is a NaN, then the
output shall contain at least one NaN component. One exception to the
NaN-in implies NaN-out rule is (0,0)/(nan,fin) which is (0,0), not
(0,nan) nor (nan,nan), because zero/finite and zero/infinite are zero,
and (nan,fin) is a vector at least as long as fin. Another exception
is (fin,fin)/(nan,fin), if replacing the nan with zero results in
(zero,zero) due to underflows.
Use the 'real' cases to guide the 'complex' cases:
0.0 Fin INF NaN * 0.0 Fin INF NaN /
---------------+ ---------------+
0.0 0.0 NaN NaN| 0.0 NaN INF INF NaN| 0.0
0.0 Fin INF NaN| Fin 0.0 Fin INF NaN| Fin
NaN INF INF NaN| INF 0.0 0.0 NaN NaN| INF
NaN NaN NaN NaN| NaN NaN NaN NaN NaN| NaN
Zero components in the result have their sign determined by the usual
rectangular formula with any non-zero input component replaced with
copysign(0.0,non-zero). This allows directed rounding and/or (-0)+(-0)
to do the right thing.
A NaN component can be replaced with 0, +inf, -inf or any finite
number in calculations to see if (part of the) result is independent
of the NaN, or, to see if the length of the result is independent of
the NaN.
If the intermediate calculations (with NaNs treated as zero) results
in (inf,inf), but there is a NaN in one of the input components, then
there is a need to have one of the result components set to NaN.
Originaly, Jim and Fred set the imaginary one to NaN (this is
arbitrary, but kind of matches cproj). But, after some more thinking
about specific cases, the reference implementations have been upgraded
to consider the directions of both the input vectors. One example is:
(inf,nan)*(fin,inf). If one assumes that inf is stronger than nan and
fin, this is as if (inf,0)*(0,inf), or x-axis*y-axis which is y-axis.
Therefore, the output should be of the form (nan,inf), that is, the
inf is on the y-axis. Taking the next step, consider
(inf,nan)*(min,max). This is "close" to (inf,0)*(0,max), or
x-axis*y-axis which is y-axis. Therefore, the output should be of the
form (nan,inf).
Since there are four components in (a+I*b)*(c+I*d), and each component
can have four classifications, there are 4**4 or 256 cases to
consider. The same is true for divide. That is why the reference
implementation is so large.
* Assumptions about numbers being on an axis.
A complex number of the form a+0*I is a vector on the x-axis. A
complex number of the form 0+b*I is a vector on the y-axis.
Thinking in polar terms, the following (ideas about vectors on an
axis) are easy to see.
Vector on x-axis times vector on x-axis is on the x-axis.
(a,0)*(c,0) => (a*c,0)
Vector on y-axis times vector on y-axis is on the x-axis.
(0,b)*(0,d) => (-b*d,0)
Vector on x-axis times vector on y-axis is on the y-axis.
(a,0)*(0,d) => (0,a*d)
(0,b)*(c,0) => (0,b*c)
Vector on x-axis times arbitrary vector is scaled length and either
same direction, or rotated by 180 degrees.
(a,0)*(c,d) => (a*c,a*d)
(a,b)*(c,0) => (a*c,b*c)
Vector on y-axis times arbitrary vector is scaled length and rotated
by +/- 90 degrees.
(0,b)*(c,d) => (-b*d,b*c)
(a,b)*(0,d) => (-b*d,a*d)
Vector on x-axis divided by vector on x-axis is on the x-axis.
(a,0)/(c,0) => (a/c,0)
Vector on y-axis divided by vector on y-axis is on the x-axis.
(0,b)/(0,d) => (b/d,0)
Vector on x-axis divided by vector on y-axis is on the y-axis.
(a,0)/(0,d) => (0,-a/d)
Vector on y-axis divided by vector on x-axis is on the y-axis.
(0,b)/(c,0) => (0,b/c)
Arbitrary vector divided by vector on x-axis is scaled length and
either same direction, or rotated by 180 degrees.
(a,b)/(c,0) => (a/c,b/c)
Arbitrary vector divided by vector on y-axis is scaled length and
rotated by +/- 90 degrees.
(a,b)/(0,d) => (b/d,-a/d)
The above are true even if one or more of a,b,c,d are INFINITY or NAN.
If one is thinking in rectangular terms, then one must assume that
zero*inf and zero*nan are zero for the above ideas to be true.
* NODIRECTION
Assuming that (+/-0, +/-0) are consistent with atan2() and carg() leads
to signed zeros being on the x-axis:
(+0,+0) => +0 deg => zero*(+1,+0)
(+0,-0) => -0 deg => zero*(+1,-0)
(-0,+0) => +180 deg => zero*(-1,+0)
(-0,-0) => -180 deg => zero*(-1,-0)
It is as if they were a vector of length zero and direction (+/-1,+/-0).
A possible advantage to using x-axis directions for zeros would be
creal(z op w) = creal(z) op creal(w) for z, w on the x-axis
But, another viewpoint is (+/-0, +/-0) represents a vector (with zero
length) in one of the quadrants:
(+0,+0) => 1st quad, +0 deg to +90 deg, with average of +45 deg
(+0,-0) => 4th quad, -0 deg to -90 deg, with average of -45 deg
(-0,+0) => 2nd quad, +90 deg to +180 deg, with average of +135 deg
(-0,-0) => 3rd quad, -90 deg to -180 deg, with average of -135 deg
This is meaning taken for "directionless". This makes them act more
like atan2(+/-INF,+/-INF). So, perhaps, atan2(+/-0, +/-0) needs to be
revisited.
If using the usual rectangular formula (with unlimited precision and
range) results in no NaNs, then use that result. This happens if all
of a,b,c,d are zero or finite or INF (as long as there is no 0*INF or
INF-INF). But, what about the case of (0,0)*(0,-1)? The usual formula
results in (0,0), while the x-axis*y-axis special case formula results
in (0,-0). The same sort of difference appears with (1,1)*(-0,0),
(1,1)*(-0,-0), and (0,-0)*(-1,1). This problem only happens when one
of the inputs is ( +/-0, +/-0 ) and the other input is ( finite,
finite ). This is one of two cases where NODIRECTION result differs
from !NODIRECTION. This also is why the reference code has
EXCLUDE_NORMAL.
The other case where NODIRECTION result differs from !NODIRECTION is
when result is (nan,0) or (0,nan) versus (nan,nan). This happens when
one input is (nan,0) or (0,nan) and other is (0,non-finite) or
(non-finite,0); or when one input is (0,0) and other is
(non-finite,finite) or (finite,non-finite).
* CLOSE
Assuming that (+/-INF, +/-finite), (+/-finite, +/-INF), and (+/-INF,
+/-INF) are consistent with atan2() and carg() leads to them being
"exactly" N*45 degrees (where N is +/-0, +/-1, +/-2, +/-3, +/-4).
This assumption was based upon the existing Annex G code that "boxes"
the infinities. That is,
(+INF,+INF) => +45 deg => INF*(+1,+1)
(+INF,-INF) => -45 deg => INF*(+1,-1)
(-INF,+INF) => +135 deg => INF*(-1,+1)
(-INF,-INF) => -135 deg => INF*(-1,-1)
(+INF,+fin) => +0 deg => INF*(+1,+0)
(+INF,-fin) => -0 deg => INF*(+1,-0)
(-INF,+fin) => +180 deg => INF*(-1,+0)
(-INF,-fin) => -180 deg => INF*(-1,-0)
(+fin,+INF) => +90 deg => INF*(+0,+1)
(+fin,-INF) => -90 deg => INF*(+0,-1)
(-fin,+INF) => +90 deg => INF*(-0,+1)
(-fin,-INF) => -90 deg => INF*(-0,-1)
This, along with the assumption that N*90 degrees for a direction
component result means a zero component result, appears to be a bad
pair of assumptions. Consider (inf+inf*I)*(inf+inf*I) and
(inf+inf*I)/(inf+inf*I). With these assumptions, an inf-inf is done
and it results in zero, instead of NaN.
Consider (inf+I)*(inf+I). Using the normal rectangular formula gets
one (inf+inf*I). But, if one treats (inf+I) as "exactly" N*45
degrees, then this becomes, in hybrid terms, inf*(1+0*I)*inf*(1+0*I),
which is inf+0*I; which is inconsistent. Therefore, (+/-INF, +/-INF),
(+/-INF, +/-finite), (+/-finite, +/-INF) are "close" to N*45 degrees
is a better assumption then they are "exactly" N*45 degrees.
Both Fred and Jim have settled on CLOSE being 1 as the most useful.
* NAN_SIGN
In complex multiply, the expression (a*c-b*d) appears. This can be
computed in three ways:(a*c - b*d), (a*c + (-b)*d) and (a*c + b*(-d)).
They appear to produce the same result for all non-NaN results; but
they can produce different NaN results (depends upon how the sign bits
of NaNs are done in the hardware, which is something not specified in
IEEE-754 nor C99).
If a result component is non-NAN, but involves a NaN input, there are
three ways to treat that case:
Use the sign of the NaN (NAN_SIGN is 0 in sample code)
Ignore the term that has the NaN (NAN_SIGN is 1 in sample code)
Ignore the expression that has the NaN, use +0 instead (NAN_SIGN is 2).
Both Fred and Jim have settled on NAN_SIGN being 1 as the most useful.
Some IEEE-754 implementations set the sign of the generated NaN from
0*INF, 0/0, INF/INF as the XOR of the input operand signs. This
retains a bit of information in the output result. Other IEEE-754
implementations use the same bit pattern, including the sign bit, for
the generated NaNs from those invalid operations. Therefore, we cannot
count on the sign bit of a NaN to carry any useful information in all
cases. Hence, it is best to ignore the NaN sign bit in computing
non-NaN results (that is, do not use NAN_SIGN being 0). One such case
is (finite,finite)/(INF,NaN). Using the terms without the NaN sign do
use the signs of the other operands, so produces zero with a sign that
has some information. Ignoring all terms that containt the NaN sign
and just using +0 appears to throw away information. So, NAN_SIGN
being 1 appears to be the portable method that maintains the most
information.
Another problem with the sign bit of NaNs is:
b*c - a*d,
b*c +(-a)*d,
0 - a*d,
0 +(-a)*d,
-(a*d),
(-a)*d,
are not all the same if one of 'a' or 'd' is a NaN and b*c is not a
NaN. A negate usually flips the sign bit of a NaN, while a subtract
usually does not alter the NaN's sign bit. Using b*c-a*d is the most
consistent.
* Identities checked
What effect do any of the above assumptions have on some identities?
Consider:
w*z == z*w
This one passes our reference implementations independent of
assumptions.
w*z == (-w)*(-z)
This fails with (nan,nan)*(nan,nan) if one considers the sign of
NaN. That is w*z is (nan,nan), while (-w)*(-z) is (-nan,-nan) [at
least on machines where (-nan)*(-nan) is -nan, rather than +nan].
If you ignore the sign of NaN, then it passes, independent of
assumptions.
-(w*z) == (-w)*z
This fails for (1,1)*(1,1) => (-0,-2) != (0,-2).
This fails for (nan,nan)*(nan,nan) => (-nan,-nan) != (nan,nan).
It passes if you ignore signs of NaNs and signs of zeros.
conj(w*z) == conj(w) * conj(z)
This one has problems with w=(a+b*I) and z=(a-b*I), even if a is 0.
That is, conj(w*z) is (a*a+b*b-0*I), while conj(w)*conj(z) is
(a*a+b*b+0*I). Here, conj(w) is z and conj(z) is w, so,
conj(w)*conj(z) is z*w, which is same as w*z, so cannot equal
conj(w*z). If you ignore the signs of NaNs and zeros, then it
passes, independent of assumptions.
cimag(z*conj(z)) == +0
If one assumes !CLOSE, then (1,inf)*(1,-inf), (inf,1)*(inf,-1),
(inf,inf)*(inf,-inf) pass; otherwise they fail with nan != +0.
If one assumes !NODIRECTION, then (0,nan)*(0,-nan), and
(nan,0)*(nan,-0) pass; otherwise they fail with nan != +0.
Independent of CLOSE and NODIRECTION, (1,nan)*(1,-nan),
(nan,1)*(nan,-1), (inf,nan)*(inf,-nan), and (nan,inf)*(nan,-inf)
all fail. If you ignore the cases where the cimag() is nan, then
it passes.
carg(z*w) == carg(z) + carg(w); with both sides mod 2*pi
If one assumes !CLOSE, then (1,1)*(inf,inf) and (1,inf)*(1,-inf)
pass; otherwise they fail with nan != finite number.
If one assumes !NODIRECTION, then (0,0)*(-0,-0) and (0,0)*(-1,-0)
pass; otherwise they fail with -0 != -pi.
Independent of CLOSE and NODIRECTION, (0,0)*(0,-1), (0,0)*(1,-1),
(0,0)*(0,1), (0,0)*(0,inf), (0,0)*(1,1), (0,0)*(-1,-1),
(0,inf)*(1,inf), and (1,inf)*(inf,0) all fail.
If one assumes carg((zero,zero)), carg((inf,inf)), carg((inf,fin)),
and carg((fin,inf)) are all NaN, then this identity is true.
z/z == 1+0i
This fails for (nan,any), (any,nan), (0,0), (any,inf), (inf,any)
z*(1/z) == 1+0i
(nan,any), (any,nan), (0,0), (any,inf), (inf,any), (0,den), (den,0),
(den,den) all fail [some because of limited range].
z/(1+0i) == z
(-0,fin), (-0,inf), (-0,-0), (-fin,-0), (-inf,-0) with sign of zero
problems.
1/(1/z) == z
Some common (independent of assumptions) "failures"
(nan,fin) => (nan,nan)
(fin,nan) => (nan,nan)
(+0.,den) => (+0,+0)
(+0.,max) => (+0,inf)
(den,den) => (+0,+0)
(den,max) => (+0,inf)
(max,1.0) => (inf,+0)
(max,max) => (inf,inf)
(inf,max) => (inf,+0)
(1.0,-0.) => (1.0,+0)
(max,-0.) => (inf,+0)
(-1.,max) => (-0.,inf)
On the other hand, do any of the specific failures matter for real
world applications? We do not know.
* Identities still to be checked
z*(1+0i) == z
z/w == (-z)/(-w)
-(z/w) == (-z)/w
(z/w)*w = z = (z*w)/w
conj(z/w) == conj(z)/conj(w)
cabs(z/w) == cabs(z) / cabs(w)
carg(z/w) == carg(z) - carg(w); with both sides mod 2*pi
* Proposal number 2: Replace G.5.1, Example 1 and 2 with the following.
** Proposed replacement for Annex G.5.1, Example 1, _Cmultf
#include <float.h> /* FLT_MAX_EXP */
#include <math.h> /* copysign(), fabs(), fmax(), logb(), isinf(),
* nextafter(), scalbnl(), isnan(), isfinite() */
#include <complex.h> /* creal(), cimag() */
#ifndef _Imaginary_I
/*
* Create a complex number from a pair of real numbers.
* Needed since x+I*y fails if I is complex and y is NaN or Inf.
* This shows that C really needs a better way to create a complex number.
*/
float complex cmplx_f( float re, float im ){
union {
float complex c;
float r[2];
} pair;
pair.r[0] = re;
pair.r[1] = im;
return pair.c;
}
#endif
#define SZ(x) copysignf(0.f,x) /* needed for Signed Zeros */
/*
* Multiply z * w...; Improved code for C99 Annex G
*/
float complex _Cmultf(float complex z, float complex w){
#pragma STDC FP_CONTRACT OFF
/*
* This value is chosen so that (scaled max)*(scaled max) is max/2.
* Needed for (min,max)*(1,max) and (max,max)*(max,max).
*/
static const float half = (FLT_MAX_EXP-3)/2;
float a, b, c, d, ac, bd, ad, bc, x, y;
float logbw, logbz, savednan;
long int ilogbz = 0L; /* z scale factor */
long int ilogbw = 0L; /* w scale factor */
int infnan = 0; /* (inf,nan) or (nan,inf) indicator */
int nanfin = 0; /* (fin,nan) or (nan,fin) indicator */
a = creal(z); b = cimag(z);
c = creal(w); d = cimag(w);
if( 0.f == d ){ /* w on x-axis; includes (+/-0,+/-0) */
if( 0.f == b ){ /* z on x-axis; includes (+/-0,+/-0) => (ac,0) */
x = a*c; y = b*SZ(c)+SZ(a)*d;
}else if( 0.f == a ){ /* z on y-axis => (0,bc) */
x = a*SZ(c)-SZ(b)*d; y = b*c;
}else{ /* z off axis => (ac,bc) */
x = a*c; y = b*c;
}
}else if( 0.f == c ){ /* w on y-axis */
if( 0.f == b ){ /* z on x-axis; includes (+/-0,+/-0) => (0,ad) */
x = SZ(a)*c-b*SZ(d); y = a*d;
}else if( 0.f == a ){ /* z on y-axis => (-bd,0) */
x = -b*d; y = SZ(b)*c+a*SZ(d);
}else{ /* z off axis => (-bd,ad) */
x = -b*d; y = a*d;
}
}else if( 0.f == b ){ /* z on x-axis; includes (+/-0,+/-0) => (ac,ad) */
x = a*c; y = a*d;
}else if( 0.f == a ){ /* z on y-axis => (-bd,bc) */
x = -b*d; y = b*c;
}else{
/* At this point, there are no zero components in either z or w.
* NaN could be replaced with +/-INF, +/-finite, +/-0. Use 0 as
* "middle" of direction choices.
*/
if ( isinf(a) || isinf(b) ) { /* z is infinite (INF,INF), (INF,NaN), (INF,finite) */
if (isnan(a) || isnan(b)){
savednan = a*b;
infnan = 1; /* inf with nan */
if (isnan(a)) a = copysignf(0.f, a); /* needed for: (nan,inf)*(1,1) */
if (isnan(b)) b = copysignf(0.f, b); /* needed for: (inf,nan)*(1,1) */
}else{
if (isfinite(a)) a = copysignf(1.f, a); /* needed for: (max,inf)*(1,-max) */
if (isfinite(b)) b = copysignf(1.f, b); /* needed for: (inf,max)*(-max,max) */
}
}else{ /* z is (finite,finite) or (NaN,finite) */
if (isnan(a) || isnan(b)) { /* z is (NaN,finite) */
savednan = a*b;
nanfin = 1; /* nan with fin */
if (isnan(a)) a = copysignf(0.f, a); /* needed for: (nan,1)*(1,inf) */
if (isnan(b)) b = copysignf(0.f, b); /* needed for: (max,nan)*(max,1) */
}
logbz = logb(fmax(fabs(a), fabs(b)));
ilogbz = (long)(logbz-half);
a = scalbnl(a, -ilogbz); b = scalbnl(b, -ilogbz); /* scale finites */
if(0.f==a && !isnan(creal(z))) a=nextafterf(a,creal(z)); /* (small,large) => (min,sqrt(max)); not (0,sqrt(max)) */
if(0.f==b && !isnan(cimag(z))) b=nextafterf(b,cimag(z)); /* needed for: (max,min)*(1,inf) */
}
/* Do same transformations to w as to z */
if ( isinf(c) || isinf(d) ) { /* w is infinite (INF,INF), (INF,NaN), (INF,finite) */
if (isnan(c) || isnan(d)){
savednan = c*d;
infnan = 1; /* inf with nan */
if (isnan(c)) c = copysignf(0.f, c); /* needed for: (1,1)*(nan,inf) */
if (isnan(d)) d = copysignf(0.f, d); /* needed for: (1,inf)*(inf,nan) */
}else{
if (isfinite(c)) c = copysignf(1.f, c); /* needed for: (1,max)*(max,-inf) */
if (isfinite(d)) d = copysignf(1.f, d); /* needed for: (nan,max)*(inf,max) */
}
}else{ /* w is (finite,finite) or (NaN,finite) */
if (isnan(c) || isnan(d)) { /* w is (NaN,finite) */
savednan = c*d;
nanfin = 1; /* nan with fin */
if (isnan(c)) c = copysignf(0.f, c); /* needed for: (1,inf)*(nan,1) */
if (isnan(d)) d = copysignf(0.f, d); /* needed for: (max,1)*(max,nan) */
}
logbw = logb(fmax(fabs(c), fabs(d)));
ilogbw = (long)(logbw-half);
c = scalbnl(c, -ilogbw); d = scalbnl(d, -ilogbw); /* scale finites */
if(0.f==c && !isnan(creal(w))) c=nextafterf(c,creal(w)); /* (small,large) => (min,sqrt(max)); not (0,sqrt(max)) */
if(0.f==d && !isnan(cimag(w))) d=nextafterf(d,cimag(w)); /* needed for: (inf,1)*(max,min) */
}
/* Compute direction; avoid 0*inf; 0 comes from nan */
ac = (a && c) ? a*c : 0.f;
ad = (a && d) ? a*d : 0.f;
bc = (b && c) ? b*c : 0.f;
bd = (b && d) ? b*d : 0.f;
x = ac - bd; /* may do INF-INF */
y = bc + ad;
if(ilogbw+ilogbz){ /* rescale the finites */
x = scalbnl(x, ilogbw+ilogbz);
y = scalbnl(y, ilogbw+ilogbz);
}
if( infnan /* add back in NaN from inf+nan*I or nan+inf*I */
|| nanfin ){ /* add back in NaN from fin+nan*I or nan+fin*I */
if( !isinf(x) )
x = savednan; /* needed for: (NaN,1)*(1,NaN); (NaN,1)*(1,1);
* (NaN,max)*(max,NaN); (NaN,max)*(max,1) */
if( !isinf(y) || isinf(x) ) /* (inf,inf) => (inf,nan) */
y = savednan; /* needed for: (NaN,1)*(NaN,1); (NaN,1)*(1,1);
* (NaN,max)*(NaN,max); (NaN,max)*(1,max) */
}
}
#ifdef _Imaginary_I
return x + _Imaginary_I * y;
#else
return cmplx_f(x,y); /* was x + I*y; but has problems with non-imaginary I */
#endif
}
** Proposed replacement for Annex G.5.1, Example 2, _Cdivf
#include <float.h> /* FLT_MAX_EXP */
#include <math.h> /* copysign(), fabs(), fmax(), logb(), isinf(),
* nextafter(), scalbnl(), isnan(), isfinite() */
#include <complex.h> /* creal(), cimag() */
#ifndef _Imaginary_I
/*
* Create a complex number from a pair of real numbers.
* Needed since x+I*y fails if I is complex and y is NaN or Inf.
* This shows that C really needs a better way to create a complex number.
*/
float complex cmplx_f( float re, float im ){
union {
float complex c;
float r[2];
} pair;
pair.r[0] = re;
pair.r[1] = im;
return pair.c;
}
#endif
#define SZ(x) copysignf(0.f,x) /* needed for Signed Zeros */
/*
* Divide z / w...; Improved code for C99 Annex G
*/
float complex _Cdivf(float complex z, float complex w){
#pragma STDC FP_CONTRACT OFF
/*
* This value is chosen so that (scaled max)*(scaled max) is max/2.
* So that c*c + d*d is as large as possible, but does not overflow.
*/
static const float half = (FLT_MAX_EXP-3)/2;
float a, b, c, d, logbw, denom, x, y;
int ilogbw = 0;
a = creal(z); b = cimag(z);
c = creal(w); d = cimag(w);
if( 0.f == d ){ /* w on x-axis; includes (+/-0,+/-0) */
if( 0.f == b ){ /* z on x-axis; includes (+/-0,+/-0) => (a/c,0) */
x = a/c; y = b*SZ(c)-SZ(a)*d;
}else if( 0.f == a ){ /* z on y-axis => (0,b/c) */
x = a*SZ(c)+SZ(b)*d; y = b/c;
}else{ /* z off axis => (a/c,b/c) */
x = a/c; y = b/c;
}
}else if( 0.f == c ){ /* w on y-axis */
if( 0.f == b ){ /* z on x-axis; includes (+/-0,+/-0) => (0,-a/d) */
x = SZ(a)*c+b*SZ(d); y = -a/d;
}else if( 0.f == a ){ /* z on y-axis => (b/d,0) */
x = b/d; y = SZ(b)*c-a*SZ(d);
}else{ /* z off axis => (b/d,-a/d) */
x = b/d; y = -a/d;
}
}else{
logbw = logb(fmax(fabs(c), fabs(d))); /* fmax: NAN < 0.f < denorm < normal < INFINITY */
if (isfinite(logbw)) {
ilogbw = (int)(logbw-half);
c = scalbnl(c, -ilogbw); d = scalbnl(d, -ilogbw);
if(0.f==c) c=nextafterf(c,creal(w)); /* (small,large) => (min,sqrt(max)); not (0,sqrt(max)) */
if(0.f==d) d=nextafterf(d,cimag(w)); /* needed for (0,inf)/(min,max) */
}
denom = c * c + d * d;
x = scalbnl((a * c + b * d) / denom, -ilogbw);
y = scalbnl((b * c - a * d) / denom, -ilogbw);
/* Recover infinities and zeros that computed as NaN and/or iNaN; */
/* the only cases are zero/nan, nan/finite, infinite/finite, and finite/infinite, ... */
if (isnan(x) || isnan(y)) { /* || needed for (inf,inf)/(1,max) */
if( (0.f==a) && (0.f==b) && isnan(d) ){ /* (0,0)/(non-nan,nan) */
x = a/c; y = b/c; /* div by "x-axis" */
}else if( (0.f==a) && (0.f==b) && isnan(c) ){ /* (0,0)/(nan,non-nan) */
x = b/d; y = -a/d; /* div by "y-axis" => rotate by 90 deg */
}else if( isnan(a) && isfinite(b) && isfinite(c) && isfinite(d) ){/* (nan,fin)/(fin,fin) */
x = scalbnl((b*d) / denom, -ilogbw); /* treat 'a' as zero for direction */
y = scalbnl((b*c) / denom, -ilogbw); /* may overflow to infinity */
if( isinf(x) && isinf(y) ){ /* (inf,inf) => (nan,inf) or (inf,nan) */
if( fabs(d) < fabs(c) ){
x = a*c; /* div by "x-axis" => (nan,inf) */
}else{
y = -a*d; /* div by "y-axis" => rotate by 90 deg (inf,nan) */
}
}else{
if( !isinf(x) ) x = a;
if( !isinf(y) ) y = a;
}
}else if( isfinite(a) && isnan(b) && isfinite(c) && isfinite(d) ){/* (fin,nan)/(fin,fin) */
x = scalbnl(( a*c) / denom, -ilogbw); /* treat 'b' as zero for direction */
y = scalbnl((-a*d) / denom, -ilogbw); /* may overflow to infinity */
if( isinf(x) && isinf(y) ){ /* (inf,inf) => (nan,inf) or (inf,nan) */
if( fabs(d) <= fabs(c) ){
y = b*c; /* div by "x-axis" => (inf,nan) */
}else{
x = b*d; /* div by "y-axis" => rotate by 90 deg (nan,inf) */
}
}else{
if( !isinf(x) ) x = b;
if( !isinf(y) ) y = b;
}
}else if( (isinf(a) || isinf(b)) &&
isfinite(c) && isfinite(d) ){
if( isnan(a) || isnan(b) ){ /* (inf,nan)/(fin,fin) or (nan,inf)/(fin,fin) */
if( (fabs(d) < fabs(c)) || ((fabs(d) == fabs(c))&&(isinf(a))) ){
x = a/c; y = b/c; /* div by "x-axis" */
}else{
x = b/d; y = -a/d; /* div by "y-axis" => rotate by 90 deg */
}
}else{ /* inf/fin */
a = copysignf(isinf(a) ? 1.f : 0.f, a);
b = copysignf(isinf(b) ? 1.f : 0.f, b);
x = INFINITY * ( a * c + b * d );
y = INFINITY * ( b * c - a * d );
}
}else if( (INFINITY==logbw) && /* (inf,non-zero) or (non-zero,inf) */
isfinite(a) && isfinite(b) ){ /* (fin,fin), (fin,0), (0,fin), (0,0) */
if( (0.f==a) && (0.f==b) ) a = copysignf(1.f,a); /* treat (0,0) as on x-axis */
c = copysignf(isinf(c) ? 1.f : 0.f, c);
d = copysignf(isinf(d) ? 1.f : 0.f, d);
x = 0.f * ( a * c + b * d );
y = 0.f * ( b * c - a * d );
}
}
}
#ifdef _Imaginary_I
return x + _Imaginary_I * y;
#else
return cmplx_f(x,y); /* was x + I*y; but has problems with non-imaginary I */
#endif
}
* Proposal number 3: Replace G.5.1, Example 1 and 2 with the following code.
** Reference implementation of cmult and cdiv
In comparing the two initial reference implementations done by Jim and
Fred, Fred discovered that they produced different results for the
same operands in a few cases, and after much analysis, he determined
that the differences were due to three sets of initial assumptions:
1) (+/-0, +/-0), (NaN, +/-0) are on the x-axis and (+/-0, NaN) is on
the y-axis [Fred] versus they are "directionless" [Jim]. We both
agree that (+/-finite,0) and (+/-INF,0) are on the x-axis, and that
(0,+/-finite) and (0,+/-INF) are on the y-axis.
2) (+/-INF, +/-INF), (+/-INF, +/-finite), (+/-finite, +/-INF) are
"close" to N*45 degrees [Jim] versus "exactly" N*45 degrees [Fred].
We have since concluded that "close" makes the most sense for most
cases (even though there are a few cases where "exact" gives a better
answer).
3) The sign bit of NaNs is [Fred] or is not [Jim] used. We have
since concluded that ignoring the NaN's sign bit makes the most sense
for most implementations.
Fred then modified his reference implementations to work with any
combination of those assumptions. That code is found next.
In the following code, intermediate calculations are done using more
precision and more exponent range so that there are no spurious
overflows, no spurious underflows, and (almost) no spurious inexacts.
For this specific reference code, use float for input and output, and
double for intermediate.
#if 1 /* Fred's current assumptions */
#define NODIRECTION 0 /* (0,0), (NaN,0) are x-axis, (0,NaN) is y-axis */
#define CLOSE 1 /* (INF,INF), (INF,fin), (fin,INF) are "close" to N*45 degrees */
#define NAN_SIGN 1 /* Ignore term with signed NaN */
#else /* Jim's current(?) assumptions */
#define NODIRECTION 1 /* (0,0), (NAN,0), (0,NAN) are "directionless" */
#define CLOSE 1 /* (INF,INF), (INF,fin), (fin,INF) are "close" to N*45 degrees */
#define NAN_SIGN 1 /* Ignore term with signed NaN */
#endif
#include <float.h> /* FLT_MAX */
#include <math.h> /* fabs, isinf, isnan, isfinite */
#include <complex.h> /* complex, _Imaginary_I, creal, cimag */
#ifndef _Imaginary_I
/*
* Create a complex number from a pair of real numbers.
* Needed since x+I*y fails if I is complex and y is NaN or Inf.
* This shows that C really needs a better way to create a complex number.
*/
float complex cmplx_f( float re, float im ){
union {
float complex c;
float r[2];
} pair;
pair.r[0] = re;
pair.r[1] = im;
return pair.c;
}
#endif
/* Classification as Zero, Finite, Infinite, or Nan */
#define Z_ 1u
#define F_ 2u
#define I_ 4u
#define N_ 8u
#define CLASSIFY(x) (isfinite(x)?(((x)==0.0)?Z_:F_):(isinf(x)?I_:N_))
/* For classification of four components */
#define P(i,j,k,m) ((i<<12)|(j<<8)|(k<<4)|m)
/*****************************************************************
* Multiply z * w...; Fred's reference implementation #2
* (a*c - b*d) + (b*c + a*d)*I
*****************************************************************/
float complex cmult2(float complex z, float complex w){
#pragma STDC FP_CONTRACT OFF
float a, b, c, d, x, y;
double ac, bd, ad, bc; /* extra precision and exponent range */
double aa, bb, cc, dd;
unsigned int cp;
a = creal(z);
b = cimag(z);
c = creal(w);
d = cimag(w);
#ifdef EXCLUDE_NORMAL
aa = a; /* increase precision and exponent range */
bb = b; /* to avoid spurious overflows, underflows */
cc = c;
dd = d;
ac = aa * cc;
bd = bb * dd;
ad = aa * dd;
bc = bb * cc;
x = (ac - bd); /* may overflow or underflow */
y = (bc + ad); /* may overflow or underflow */
if( isnan(x) || isnan(y) ) /* fixup any undeserved NaNs */
#endif
{
/* E.g. P(Z_,F_,I_,N_) indicates a is zero, b finite, c infinite, and d NaN */
cp = P(CLASSIFY(a),CLASSIFY(b),CLASSIFY(c),CLASSIFY(d));
switch (cp) {
case P(F_,F_ , F_,F_): /* fin * fin => no NaN */
{
aa = a; /* increase precision and exponent range */
bb = b; /* to avoid spurious overflows, underflows */
cc = c;
dd = d;
ac = aa * cc;
bd = bb * dd;
ad = aa * dd;
bc = bb * cc;
x = (ac - bd); /* may overflow or underflow */
y = (bc + ad); /* may overflow or underflow */
break;
}
/* (Z ,Z , Z ,Z ) */
case P(Z_,Z_ , Z_,Z_): /* zero * zero => no NaN */
/* (* ,Z , Z ,Z ) */
case P(F_,Z_ , Z_,Z_): /* fin * zero => no NaN */
case P(I_,Z_ , Z_,Z_): /* inf * zero => nan */
/* (Z ,Z , * ,Z ) */
case P(Z_,Z_ , F_,Z_): /* zero * fin => no NaN */
case P(Z_,Z_ , I_,Z_): /* zero * inf => nan */
#if NODIRECTION
x = a*c - b*d;
y = b*c + a*d;
break;
#endif
/* (* ,Z , * ,Z ) */
case P(F_,Z_ , F_,Z_): /* fin * fin => no NaN */
case P(F_,Z_ , I_,Z_):
case P(I_,Z_ , F_,Z_):
case P(I_,Z_ , I_,Z_):
/* x-axis * x-axis => x-axis; (a*c, y is signed zero) */
x = a*c;
y = b*O(c) + O(a)*d;
break;
/* (* ,Z , * ,Z ) */
case P(N_,Z_ , Z_,Z_): /* "inf" * zero => nan */
#if NODIRECTION
x = a*c - b*d;
y = b*c + a*d;
#else /* x-axis * x-axis => x-axis; (a*c, y is signed zero) */
x = a*c;
#if 2 == NAN_SIGN
y = 0.f; /* +zero */
#elif 1 == NAN_SIGN
y = b*O(c); /* signed zero w/o NaN term */
#else
y = b*O(c) + O(a)*d; /* signed zero */
#endif
#endif
break;
case P(N_,Z_ , F_,Z_): /* "zero" * fin */
case P(N_,Z_ , I_,Z_): /* "zero" * inf => nan */
x = a*c;
#if NODIRECTION
y = b*c; /* as if any * x-axis */
#else
#if 2 == NAN_SIGN
y = 0.f; /* +zero */
#elif 1 == NAN_SIGN
y = b*O(c); /* signed zero w/o NaN term */
#else
y = b*O(c) + O(a)*d; /* signed zero */
#endif
#endif
break;
case P(Z_,Z_ , N_,Z_): /* zero * "inf" => nan */
#if NODIRECTION
x = a*c - b*d;
y = b*c + a*d;
#else /* x-axis * x-axis => x-axis; (a*c, y is signed zero) */
x = a*c;
#if 2 == NAN_SIGN
y = 0.f; /* +zero */
#elif 1 == NAN_SIGN
y = O(a)*d; /* signed zero w/o NaN term */
#else
y = b*O(c) + O(a)*d; /* signed zero */
#endif
#endif
break;
case P(F_,Z_ , N_,Z_): /* fin * "zero" */
case P(I_,Z_ , N_,Z_): /* inf * "zero" => nan */
x = a*c;
#if NODIRECTION
y = a*d; /* as if x-axis * any */
#else
#if 2 == NAN_SIGN
y = 0.f; /* +zero */
#elif 1 == NAN_SIGN
y = O(a)*d; /* signed zero w/o NaN term */
#else
y = b*O(c) + O(a)*d; /* signed zero */
#endif
#endif
break;
case P(N_,Z_ , N_,Z_): /* "zero" * "inf" => nan */
#if NODIRECTION
x = a*c - b*d;
y = b*c + a*d;
#else /* x-axis * x-axis => x-axis; (a*c, y is signed zero) */
x = a*c;
#if (2 == NAN_SIGN) || (1 == NAN_SIGN)
y = 0.f; /* +zero */
#else
y = b*O(c) + O(a)*d; /* signed zero */
#endif
#endif
break;
/* (Z ,Z , Z ,* ) */
case P(Z_,Z_ , Z_,F_): /* zero * fin => no NaN */
case P(Z_,Z_ , Z_,I_): /* zero * inf => nan */
/* x-axis * y-axis => y-axis; (x is signed zero, a*d) */
#if NODIRECTION
x = a*c - b*d;
y = b*c + a*d;
break;
#endif
/* (* ,Z , Z ,* ) */
case P(F_,Z_ , Z_,F_): /* fin * fin => no NaN */
case P(F_,Z_ , Z_,I_):
case P(I_,Z_ , Z_,F_):
case P(I_,Z_ , Z_,I_):
/* x-axis * y-axis => y-axis; (x is signed zero, a*d) */
x = O(a)*c - b*O(d);
y = a*d;
break;
case P(Z_,Z_ , Z_,N_): /* zero * "inf" => nan */
#if NODIRECTION
x = a*c - b*d;
y = b*c + a*d;
#else /* x-axis * y-axis => y-axis; (x is signed zero, a*d) */
#if 2 == NAN_SIGN
x = 0.f; /* +zero */
#elif 1 == NAN_SIGN
x = O(a)*c; /* signed zero w/o NaN term */
#else
x = O(a)*c - b*O(d); /* signed zero */
#endif
y = a*d;
#endif
break;
case P(F_,Z_ , Z_,N_): /* fin * "zero" */
case P(I_,Z_ , Z_,N_): /* inf * "zero" => nan */
/* x-axis * y-axis => y-axis; (x is signed zero, a*d) */
#if NODIRECTION
x = a*c; /* as if x-axis * any */
#else
#if 2 == NAN_SIGN
x = 0.f; /* +zero */
#elif 1 == NAN_SIGN
x = O(a)*c; /* signed zero w/o NaN term */
#else
x = O(a)*c - b*O(d); /* signed zero */
#endif
#endif
y = a*d;
break;
/* (N ,Z , Z ,* ) */
case P(N_,Z_ , Z_,I_): /* "zero" * inf => nan */
case P(N_,Z_ , Z_,F_): /* "zero" * fin */
#if NODIRECTION
x = - (b*d); /* as if any * y-axis */
#else /* x-axis * y-axis => y-axis; (x is signed zero, a*d) */
#if 2 == NAN_SIGN
x = 0.f; /* +zero */
#elif 1 == NAN_SIGN
x = - b*O(d); /* signed zero w/o NaN term */
#else
x = O(a)*c - b*O(d); /* signed zero */
#endif
#endif
y = a*d;
break;
/* (* ,Z , Z ,* ) */
case P(N_,Z_ , Z_,N_): /* "zero" * "inf" => nan */
#if NODIRECTION
x = a*c - b*d;
y = b*c + a*d;
#else /* x-axis * y-axis => y-axis; (x is signed zero, a*d) */
#if (2 == NAN_SIGN) || (1 == NAN_SIGN)
x = 0.f; /* +zero */
#else
x = O(a)*c - b*O(d); /* signed zero */
#endif
y = a*d;
#endif
break;
/* (Z ,Z , * ,* ) */
case P(Z_,Z_ , F_,F_): /* zero * fin => no NaN */
case P(Z_,Z_ , F_,I_): /* zero * inf */
case P(Z_,Z_ , F_,N_): /* zero * "fin" */
case P(Z_,Z_ , I_,F_): /* zero * inf */
case P(Z_,Z_ , N_,F_): /* zero * "fin" */
/* x-axis * any => any; (a*c, a*d) */
#if NODIRECTION
x = a*c - b*d;
y = b*c + a*d;
break;
#endif
case P(I_,I_ , F_,N_): /* inf * "fin" => inf w/ nan */
/* (I ,F , * ,* ) */
case P(I_,F_ , F_,F_): /* inf * fin => no NaN */
case P(I_,F_ , F_,N_): /* inf * "fin" => inf w/ nan */
case P(I_,F_ , N_,F_): /* inf * "fin" => inf w/ nan */
case P(I_,F_ , I_,N_): /* inf * "inf" => inf w/ nan */
/* (* ,Z , * ,* ) */
case P(F_,Z_ , F_,F_): /* fin * fin => no NaN */
case P(F_,Z_ , F_,I_):
case P(F_,Z_ , F_,N_):
case P(F_,Z_ , I_,F_):
case P(F_,Z_ , I_,I_):
case P(F_,Z_ , I_,N_):
case P(F_,Z_ , N_,F_):
case P(F_,Z_ , N_,I_):
case P(F_,Z_ , N_,N_): /* (N,N) */
case P(I_,Z_ , F_,F_): /* inf * fin => no NaN */
case P(I_,Z_ , F_,I_):
case P(I_,Z_ , F_,N_):
case P(I_,Z_ , I_,F_):
case P(I_,Z_ , I_,I_):
case P(I_,Z_ , I_,N_):
case P(I_,Z_ , N_,F_):
case P(I_,Z_ , N_,I_):
case P(I_,Z_ , N_,N_): /* (N,N) */
/* x-axis * any => any; (a*c, a*d) */
x = a*c;
y = a*d;
break;
case P(N_,Z_ , F_,F_): /* "zero" */
#if NODIRECTION
x = a*c - b*d;
y = b*c + a*d;
#else
x = a*c;
y = a*d;
#endif
break;
/* (* ,* , Z ,Z ) */
case P(F_,F_ , Z_,Z_): /* fin * zero => no NaN */
case P(F_,I_ , Z_,Z_): /* inf * zero => nan */
case P(I_,F_ , Z_,Z_): /* inf * zero => nan */
case P(F_,N_ , Z_,Z_): /* "inf" * zero => nan */
case P(N_,F_ , Z_,Z_): /* "inf" * zero => nan */
/* any * x-axis => any; (a*c, b*c) */
#if NODIRECTION
x = a*c - b*d;
y = b*c + a*d;
break;
#endif
case P(F_,N_ , I_,I_): /* "fin" * inf => inf w/ nan */
/* (* ,* , I ,F ) */
case P(F_,F_ , I_,F_): /* fin * fin => no NaN */
case P(F_,N_ , I_,F_): /* "fin" * inf => inf w/ nan */
case P(I_,N_ , I_,F_): /* "inf" * inf => inf w/ nan */
case P(N_,F_ , I_,F_): /* "fin" * inf => inf w/ nan */
/* (* ,* , * ,Z ) */
case P(F_,F_ , F_,Z_): /* fin * fin => no NaN */
case P(F_,F_ , I_,Z_):
case P(F_,I_ , F_,Z_):
case P(F_,I_ , I_,Z_):
case P(F_,N_ , F_,Z_):
case P(F_,N_ , I_,Z_):
case P(I_,F_ , F_,Z_):
case P(I_,F_ , I_,Z_):
case P(I_,I_ , F_,Z_):
case P(I_,I_ , I_,Z_):
case P(I_,N_ , F_,Z_):
case P(I_,N_ , I_,Z_):
case P(N_,F_ , F_,Z_):
case P(N_,F_ , I_,Z_):
case P(N_,I_ , F_,Z_):
case P(N_,I_ , I_,Z_):
case P(N_,N_ , F_,Z_): /* (N,N) */
case P(N_,N_ , I_,Z_): /* (N,N) */
/* any * x-axis => any; (a*c, b*c) */
x = a*c;
y = b*c;
break;
case P(F_,F_ , N_,Z_): /* "zero" */
#if NODIRECTION
x = a*c - b*d;
y = b*c + a*d;
#else
x = a*c;
y = b*c;
#endif
break;
/* (Z ,* , Z ,Z ) */
case P(Z_,F_ , Z_,Z_): /* fin * zero => no NaN */
case P(Z_,I_ , Z_,Z_): /* inf * zero => nan */
/* y-axis * x-axis => y-axis; (x is signed zero, b*c) */
#if NODIRECTION
x = a*c - b*d;
y = b*c + a*d;
break;
#endif
/* (Z ,* , * ,Z ) */
case P(Z_,F_ , F_,Z_): /* fin * fin => no NaN */
case P(Z_,F_ , I_,Z_):
case P(Z_,I_ , F_,Z_):
case P(Z_,I_ , I_,Z_):
/* y-axis * x-axis => y-axis; (x is signed zero, b*c) */
x = a*O(c) - O(b)*d;
y = b*c;
break;
/* (Z ,N , * ,Z ) */
case P(Z_,N_ , Z_,Z_): /* "inf" * zero => nan */
#if NODIRECTION
x = a*c - b*d;
y = b*c + a*d;
#else
#if 2 == NAN_SIGN /* +zero */
x = 0.f;
#elif 1 == NAN_SIGN /* signed zero w/o NaN term */
x = a*O(c);
#else /* signed zero */
x = a*O(c) - O(b)*d;
#endif
y = b*c;
#endif
break;
case P(Z_,N_ , I_,Z_): /* "zero" * inf => nan */
case P(Z_,N_ , F_,Z_): /* "zero" * fin */
#if NODIRECTION
x = a*c; /* as if any * x-axis */
#else
#if 2 == NAN_SIGN /* +zero */
x = 0.f;
#elif 1 == NAN_SIGN /* signed zero w/o NaN term */
x = a*O(c);
#else /* signed zero */
x = a*O(c) - O(b)*d;
#endif
#endif
y = b*c;
break;
/* (Z ,* , N ,Z ) */
case P(Z_,I_ , N_,Z_): /* "inf" * "zero" => nan */
case P(Z_,F_ , N_,Z_): /* fin * "zero" */
#if NODIRECTION
x = - (b*d); /* as if y-axis * any */
#else
#if 2 == NAN_SIGN /* +zero */
x = 0.f;
#elif 1 == NAN_SIGN /* signed zero w/o NaN term */
x = - O(b)*d;
#else /* signed zero */
x = a*O(c) - O(b)*d;
#endif
#endif
y = b*c;
break;
case P(Z_,N_ , N_,Z_): /* "zero" * "inf" => nan */
#if NODIRECTION
x = a*c - b*d;
y = b*c + a*d;
#else
#if (2 == NAN_SIGN) || (1 == NAN_SIGN) /* +zero */
x = 0.f;
#else /* signed zero */
x = a*O(c) - O(b)*d;
#endif
y = b*c;
#endif
break;
/* (Z ,* , Z ,* ) */
case P(Z_,F_ , Z_,F_): /* fin * fin => no NaN */
case P(Z_,F_ , Z_,I_):
case P(Z_,I_ , Z_,F_):
case P(Z_,I_ , Z_,I_):
/* y-axis * y-axis => x-axis; (-b*d, y is signed zero) */
x = -(b*d); /* -b*d != -(b*d) if b is a NaN! */
y = O(b)*c + a*O(d);
break;
case P(Z_,F_ , Z_,N_): /* fin * "zero" */
case P(Z_,I_ , Z_,N_): /* inf * "zero" => nan */
x = -(b*d);
#if NODIRECTION
y = b*c; /* as if y-axis * any */
#else
#if 2 == NAN_SIGN /* +zero */
y = 0.f;
#elif 1 == NAN_SIGN /* signed zero w/o NaN term */
y = O(b)*c;
#else /* signed zero */
y = O(b)*c + a*O(d);
#endif
#endif
break;
case P(Z_,N_ , Z_,F_): /* "zero" * fin */
case P(Z_,N_ , Z_,I_): /* "zero" * inf => nan */
x = -(b*d);
#if NODIRECTION
y = a*d; /* as if any * y-axis */
#else
#if 2 == NAN_SIGN /* +zero */
y = 0.f;
#elif 1 == NAN_SIGN /* signed zero w/o NaN term */
y = a*O(d);
#else /* signed zero */
y = O(b)*c + a*O(d);
#endif
#endif
break;
case P(Z_,N_ , Z_,N_): /* "zero" * "inf" => nan */
#if NODIRECTION
x = a*c - b*d;
y = b*c + a*d;
#else
x = -(b*d);
#if (2 == NAN_SIGN) || (1 == NAN_SIGN) /* +zero */
y = 0.f;
#else /* signed zero */
y = O(b)*c + a*O(d);
#endif
#endif
break;
case P(I_,I_ , N_,F_): /* inf * "fin" => inf w/ nan */
/* (F ,I , * ,* ) */
case P(F_,I_ , F_,F_): /* inf * fin => no NaN */
case P(F_,I_ , N_,F_): /* inf * "fin" => inf w/ nan */
case P(F_,I_ , N_,I_): /* inf * "inf" => inf w/ nan */
/* (Z ,* , * ,* ) */
case P(Z_,F_ , F_,F_): /* fin * fin => no NaN */
case P(Z_,F_ , F_,I_):
case P(Z_,F_ , F_,N_):
case P(Z_,F_ , I_,F_):
case P(Z_,F_ , I_,I_):
case P(Z_,F_ , I_,N_):
case P(Z_,F_ , N_,F_):
case P(Z_,F_ , N_,I_):
case P(Z_,F_ , N_,N_): /* (N,N) */
case P(Z_,I_ , F_,F_): /* inf * fin => no NaN */
case P(Z_,I_ , F_,I_):
case P(Z_,I_ , F_,N_):
case P(Z_,I_ , I_,F_):
case P(Z_,I_ , I_,I_):
case P(Z_,I_ , I_,N_):
case P(Z_,I_ , N_,F_):
case P(Z_,I_ , N_,I_):
case P(Z_,I_ , N_,N_): /* (N,N) */
/* y-axis * any => any; (-b*d, b*c) */
x = -(b*d);
y = b*c;
break;
case P(F_,I_ , F_,N_): /* inf * "fin" => inf w/ nan */
x = O(a)*c - (b*d); /* nan */
y = b*c; /* inf */
break;
case P(Z_,N_ , F_,F_): /* "zero" */
#if NODIRECTION
x = a*c - b*d;
y = b*c + a*d;
#else
x = -(b*d);
y = b*c;
#endif
break;
case P(N_,F_ , I_,I_): /* "fin" * inf => inf */
/* (* ,* , F ,I ) */
case P(F_,F_ , F_,I_): /* fin * inf => no NaN */
case P(N_,F_ , F_,I_): /* "fin" * inf => inf */
case P(N_,I_ , F_,I_): /* "inf" * inf => inf */
/* (* ,* , Z ,* ) */
case P(F_,F_ , Z_,F_): /* fin * fin => no NaN */
case P(F_,F_ , Z_,I_): /* fin * inf => no NaN */
case P(F_,I_ , Z_,F_):
case P(F_,I_ , Z_,I_):
case P(F_,N_ , Z_,F_):
case P(F_,N_ , Z_,I_):
case P(I_,F_ , Z_,F_):
case P(I_,F_ , Z_,I_):
case P(I_,I_ , Z_,F_):
case P(I_,I_ , Z_,I_):
case P(I_,N_ , Z_,F_):
case P(I_,N_ , Z_,I_):
case P(N_,F_ , Z_,F_):
case P(N_,F_ , Z_,I_):
case P(N_,I_ , Z_,F_):
case P(N_,I_ , Z_,I_):
case P(N_,N_ , Z_,F_): /* (N,N) */
case P(N_,N_ , Z_,I_): /* (N,N) */
/* any * y-axis => any; (-b*d,a*d) */
x = -(b*d);
y = a*d;
break;
case P(F_,N_ , F_,I_): /* "fin" * inf => inf w/ NaN */
x = O(a)*c - (b*d); /* nan; ac needed for signed NaN */
y = a*d; /* inf */
break;
case P(F_,F_ , Z_,N_): /* "zero" */
#if NODIRECTION
x = a*c - (b*d);
y = b*c + a*d;
#else
x = -(b*d);
y = a*d;
#endif
break;
/* (I ,N , * ,* ) */
case P(I_,N_ , F_,I_): /* "inf" * inf => inf w/ nan; inf on y-axis */
x = a*c - (b*d); /* nan; ac needed for signed NaN */
y = a*d; /* inf */
break;
case P(I_,N_ , F_,F_): /* "inf" * fin => inf w/ nan */
case P(I_,N_ , I_,I_): /* "inf" * inf => inf w/ nan */
if( fabs(c) < fabs(d) ){ /* "y-axis" */
x = a*c - (b*d); /* nan; ac needed for signed NaN */
y = a*d; /* inf */
}else{ /* "x-axis" */
x = a*c; /* inf */
y = b*c; /* nan */
}
break;
/* (N ,I , * ,* ) */
case P(N_,I_ , I_,F_): /* "inf" * "inf" => inf w/ nan; inf on y-axis */
x = a*c;
y = b*c;
break;
case P(N_,I_ , F_,F_): /* "inf" * fin => inf w/ nan */
case P(N_,I_ , I_,I_): /* "inf" * inf => inf w/ nan */
if( fabs(c) <= fabs(d) ){ /* "y-axis" */
x = -(b*d); /* inf */
y = a*d; /* nan */
}else{ /* "x-axis" */
x = a*c; /* nan */
y = b*c; /* inf */
}
break;
/* (* ,* , I ,N ) */
case P(F_,I_ , I_,N_): /* inf * "inf" => inf w/ nan; inf on y-axis */
x = a*c - (b*d); /* nan; ac needed for signed NaN */
y = b*c; /* inf */
break;
case P(F_,F_ , I_,N_): /* fin * "inf" => inf w/ nan */
case P(I_,I_ , I_,N_): /* inf * "inf" => inf w/ nan */
if( fabs(a) < fabs(b) ){ /* "y-axis" */
x = a*c - (b*d); /* nan; ac needed for signed NaN */
y = b*c; /* inf */
}else{ /* "x-axis" */
x = a*c; /* inf */
y = a*d; /* nan */
}
break;
/* (* ,* , N, I ) */
case P(I_,F_ , N_,I_): /* inf * "inf" => inf w/ nan; inf on y-axis */
x = a*c;
y = a*d;
break;
case P(F_,F_ , N_,I_): /* fin * "inf" => inf w/ nan */
case P(I_,I_ , N_,I_): /* inf * "inf" => inf w/ nan */
if( fabs(a) <= fabs(b) ){ /* "y-axis" */
x = -(b*d); /* inf */
y = b*c; /* nan */
}else{ /* "x-axis" */
x = a*c; /* nan */
y = a*d; /* inf */
}
break;
/* (i ,N , N ,i ) ad is inf */
case P(F_,N_ , N_,I_): /* inf w/ nan; x-axis * y-axis => inf on y-axis */
case P(I_,N_ , N_,F_): /* "inf" * "fin" => inf */
case P(I_,N_ , N_,I_): /* "inf" * "inf" => inf */
x = a*c - b*d; /* nan */
y = a*d; /* inf */
break;
/* (* ,F , Z ,N ) */
case P(I_,F_ , Z_,N_): /* inf * "zero" */
case P(N_,F_ , Z_,N_): /* "zero" */
/* (N ,Z , F ,* ) */
case P(N_,Z_ , F_,I_): /* "zero" * inf */
case P(N_,Z_ , F_,N_): /* "zero" */
/* bc is 0 */
x = a*c - b*d;
#if NODIRECTION
y = b*c + a*d;
#else
y = a*d;
#endif
break;
/* (N ,i , i ,N ) bc is inf */
case P(N_,F_ , I_,N_): /* "fin" * "inf" => inf w/ nan */
case P(N_,I_ , F_,N_): /* "inf" * "fin" => inf w/ nan */
case P(N_,I_ , I_,N_): /* "inf" * "inf" => inf w/ nan */
x = a*c - b*d; /* nan */
y = b*c; /* inf */
break;
/* (F ,* , N ,Z ) */
case P(F_,I_ , N_,Z_): /* inf * "zero" */
case P(F_,N_ , N_,Z_): /* "zero" */
/* (Z ,N , * ,F ) */
case P(Z_,N_ , I_,F_): /* "zero" * inf */
case P(Z_,N_ , N_,F_): /* "zero" */
/* ad is 0 */
x = a*c - b*d;
#if NODIRECTION
y = b*c + a*d;
#else
y = b*c;
#endif
break;
/* (i ,N , i ,N ) ac is inf */
case P(F_,N_ , I_,N_): /* "fin" * "inf" => inf */
case P(I_,N_ , F_,N_): /* "inf" * "fin" => inf */
case P(I_,N_ , I_,N_): /* "inf" * "inf" => inf */
x = a*c; /* inf */
y = b*c + a*d; /* nan */
break;
/* (* ,F , N ,Z ) */
case P(I_,F_ , N_,Z_): /* inf * "zero" */
case P(N_,F_ , N_,Z_): /* "zero" */
/* (N ,Z , * ,F ) */
case P(N_,Z_ , I_,F_): /* "zero" * inf */
case P(N_,Z_ , N_,F_): /* "zero" */
/* bd is 0 */
#if NODIRECTION
x = a*c - b*d;
#else
x = a*c;
#endif
y = b*c + a*d;
break;
/* (N ,i , N ,i ) bd is inf */
case P(N_,F_ , N_,I_): /* "fin" * "inf" => inf */
case P(N_,I_ , N_,F_): /* "inf" * "fin" => inf */
case P(N_,I_ , N_,I_): /* "inf" * "inf" => inf */
x = -(b*d); /* inf */
y = b*c + a*d; /* nan */
break;
/* (F ,* , Z ,N ) */
case P(F_,I_ , Z_,N_): /* inf * "zero" */
case P(F_,N_ , Z_,N_): /* "zero" */
/* (Z ,N , F ,* ) */
case P(Z_,N_ , F_,I_): /* "zero" * inf */
case P(Z_,N_ , F_,N_): /* "zero" */
/* ac is 0 */
#if NODIRECTION
x = a*c - (b*d);
#else
x = -(b*d);
#endif
y = b*c + a*d;
break;
case P(F_,N_ , F_,N_): /* "fin" * "fin" => may overflow w/ nan */
x = 0.f;
ac = a; ac *= c;
if( FLT_MAX < fabs(ac) ) x = ac; /* may overflow */
if( !isinf(x) ) x = ac - b*d; /* inf stronger than nan; ac needed for signed NaN */
y = b*c + a*d;
break;
case P(F_,N_ , N_,F_): /* "fin" * "fin" => may overflow w/ nan */
x = a*c - b*d;
y = 0.f;
ad = a; ad *= d;
if( FLT_MAX < fabs(ad) ) y = ad; /* may overflow */
if( !isinf(y) ) y = b*c; /* inf stronger than nan */
break;
case P(N_,F_ , F_,N_): /* "fin" * "fin" => may overflow w/ nan */
x = a*c - b*d;
y = 0.f;
bc = b; bc *= c;
if( FLT_MAX < fabs(bc) ) y = bc; /* may overflow */
if( !isinf(y) ) y = a*d; /* inf stronger than nan */
break;
case P(N_,F_ , N_,F_): /* "fin" * "fin" => may overflow w/ nan */
x = 0.f;
bd = -b; bd *= d;
if( FLT_MAX < fabs(bd) ) x = bd; /* may overflow */
if( !isinf(x) ) x = a*c; /* inf stronger than nan */
y = b*c + a*d;
break;
case P(F_,F_ , F_,N_): /* fin * "fin" => may overflow w/ nan */
x = y = 0.f;
cc = c;
ac = a*cc;
if( FLT_MAX < fabs(ac) ) x = ac; /* may overflow */
if( !isinf(x) ) x = ac - (b*d); /* inf stronger than nan; ac needed for signed NaN */
bc = b*cc;
if( FLT_MAX < fabs(bc) ) y = bc; /* may overflow */
if( !isinf(y) ) y = a*d; /* inf stronger than nan */
if( isinf(x) && isinf(y) ){
if( fabs(a) < fabs(b) ){ /* "y-axis" */
x = ac - (b*d); /* (nan,inf); ac needed for signed NaN */
}else{ /* "x-axis" */
y = a*d; /* (inf,nan) */
}
}
break;
case P(F_,F_ , N_,F_): /* fin * "fin" => may overflow w/ nan */
x = y = 0.f;
dd = d;
bd = -(b*dd);
if( FLT_MAX < fabs(bd) ) x = bd; /* may overflow */
if( !isinf(x) ) x = a*c; /* inf stronger than nan */
ad = a*dd;
if( FLT_MAX < fabs(ad) ) y = ad; /* may overflow */
if( !isinf(y) ) y = b*c; /* inf stronger than nan */
if( isinf(x) && isinf(y) ){
if( fabs(a) <= fabs(b) ){ /* "y-axis" */
y = b*c; /* (inf,nan) */
}else{ /* "x-axis" */
x = a*c; /* (nan,inf) */
}
}
break;
case P(F_,N_ , F_,F_): /* "fin" * fin => may overflow w/ nan */
x = y = 0.f;
aa = a;
ac = aa*c;
if( FLT_MAX < fabs(ac) ) x = ac; /* may overflow */
if( !isinf(x) ) x = ac - (b*d); /* inf stronger than nan; ac needed for signed NaN */
ad = aa*d;
if( FLT_MAX < fabs(ad) ) y = ad; /* may overflow */
if( !isinf(y) ) y = b*c; /* inf stronger than nan */
if( isinf(x) && isinf(y) ){
if( fabs(c) < fabs(d) ){ /* "y-axis" */
x = ac - (b*d); /* (nan,inf); ac needed for signed NaN */
}else{ /* "x-axis" */
y = b*c; /* (inf,nan) */
}
}
break;
case P(N_,F_ , F_,F_): /* "fin" * fin => may overflow w/ nan */
x = y = 0.f;
bb = b;
bd = -(bb*d);
if( FLT_MAX < fabs(bd) ) x = bd; /* may overflow */
if( !isinf(x) ) x = a*c; /* inf stronger than nan */
bc = bb*c;
if( FLT_MAX < fabs(bc) ) y = bc; /* may overflow */
if( !isinf(y) ) y = a*d; /* inf stronger than nan */
if( isinf(x) && isinf(y) ){
if( fabs(c) <= fabs(d) ){ /* "y-axis" */
y = a*d; /* (inf,nan) */
}else{ /* "x-axis" */
x = a*c; /* (nan,inf) */
}
}
break;
case P(F_,I_ , F_,I_): /* inf * inf => inf; nan may come from inf-inf */
#if !CLOSE
x = -(b*d);
if(b*c == -(a*d)){ /* (1,inf)*(-1,inf) => (-inf,0) */
y = 0.f;
}else{ /* (1,inf)*(1,inf) => (-inf,inf) */
y = b*c + a*d;
}
break;
#endif
case P(F_,I_ , I_,F_): /* inf * inf => inf; nan may come from inf-inf */
#if !CLOSE
if(a*c == b*d){ /* (1,inf)*(inf,1) => (0,inf) */
x = 0.f;
}else{ /* (1,inf)*(inf,-1) => (inf,inf) */
x = a*c - b*d;
}
y = b*c;
break;
#endif
case P(I_,F_ , F_,I_): /* inf * inf => inf; nan may come from inf-inf */
#if !CLOSE
if(a*c == b*d){ /* (inf,1)*(1,inf) => (0,inf) */
x = 0.f;
}else{ /* (inf,-1)*(1,inf) => (inf,inf) */
x = a*c - b*d;
}
y = a*d;
break;
#endif
case P(I_,F_ , I_,F_): /* inf * inf => inf; nan may come from inf-inf */
#if !CLOSE
x = a*c;
if(b*c == -(a*d)){ /* (inf,1)*(inf,-1) => (inf,0) */
y = 0.f;
}else{ /* (inf,1)*(inf,1) => (inf,inf) */
y = b*c + a*d;
}
#else
/* 3 of the 4 products are an inf; 4th may overflow */
aa = a; /* increase precision and exponent range */
bb = b; /* to avoid spurious overflows, underflows */
cc = c;
dd = d;
ac = aa * cc; /* prevent max*max being overflow to inf */
bd = bb * dd;
ad = aa * dd;
bc = bb * cc;
x = (ac - bd); /* inf or nan (inf-inf) */
y = (bc + ad); /* inf or nan (inf-inf) */
#endif
break;
case P(F_,F_ , I_,I_): /* fin * inf => inf; nan may come from inf-inf */
case P(I_,I_ , F_,F_): /* inf * fin => inf; nan may come from inf-inf */
/* each of the 4 products is an inf */
#if !CLOSE
if( isinf(a) ) a = copysignf(1.f, a); /* "box" the infinities */
if( isinf(b) ) b = copysignf(1.f, b);
if( isinf(c) ) c = copysignf(1.f, c);
if( isinf(d) ) d = copysignf(1.f, d);
if(a*c == b*d){
x = 0.f;
}else{
x = INFINITY * (a*c - b*d);
}
if(b*c == -(a*d)){
y = 0.f;
}else{
y = INFINITY * (b*c + a*d);
}
break;
#endif
case P(F_,I_ , I_,I_): /* inf * inf => inf; nan may come from inf-inf */
case P(I_,F_ , I_,I_): /* inf * inf => inf; nan may come from inf-inf */
case P(I_,I_ , F_,I_): /* inf * inf => inf; nan may come from inf-inf */
case P(I_,I_ , I_,F_): /* inf * inf => inf; nan may come from inf-inf */
case P(I_,I_ , I_,I_): /* inf * inf => inf; nan may come from inf-inf */
/* each of the 4 products is an inf */
#if !CLOSE
a = copysignf(isinf(a) ? 1.f : 0.f, a); /* "box" the infinities */
b = copysignf(isinf(b) ? 1.f : 0.f, b);
c = copysignf(isinf(c) ? 1.f : 0.f, c);
d = copysignf(isinf(d) ? 1.f : 0.f, d);
x = a*c - b*d; /* direction of result */
y = b*c + a*d;
if(x) x *= INFINITY; /* treat 0*inf as 0 */
if(y) y *= INFINITY;
break;
#endif
case P(I_,I_ , N_,Z_): /* inf * "zero" */
case P(I_,I_ , Z_,N_): /* inf * "zero" */
case P(I_,I_ , Z_,Z_): /* inf * zero */
case P(I_,N_ , N_,Z_): /* "inf" * "zero" */
case P(I_,N_ , Z_,N_): /* "inf" * "zero" */
case P(I_,N_ , Z_,Z_): /* "inf" * zero */
case P(N_,I_ , N_,Z_): /* "inf" * "zero" */
case P(N_,I_ , Z_,N_): /* "inf" * "zero" */
case P(N_,I_ , Z_,Z_): /* "inf" * zero */
case P(N_,Z_ , I_,I_): /* "zero" * inf */
case P(N_,Z_ , I_,N_): /* "zero" * "inf" */
case P(N_,Z_ , N_,I_): /* "zero" * "inf" */
case P(Z_,N_ , I_,I_): /* "zero" * inf */
case P(Z_,N_ , I_,N_): /* "zero" * "inf" */
case P(Z_,N_ , N_,I_): /* "zero" * "inf" */
case P(Z_,Z_ , I_,I_): /* zero * inf */
case P(Z_,Z_ , I_,N_): /* zero * "inf" */
case P(Z_,Z_ , N_,I_): /* zero * "inf" */
/* each of the 4 products is a nan */
x = a*c - b*d;
y = b*c + a*d;
break;
/* (* ,* , N ,N ) */
case P(F_,F_ , N_,N_): /* (N,N) */
case P(F_,I_ , N_,N_): /* (N,N) */
case P(F_,N_ , N_,N_): /* (N,N) */
case P(I_,F_ , N_,N_): /* (N,N) */
case P(I_,I_ , N_,N_): /* (N,N) */
case P(I_,N_ , N_,N_): /* (N,N) */
case P(N_,F_ , N_,N_): /* (N,N) */
case P(N_,I_ , N_,N_): /* (N,N) */
case P(N_,Z_ , N_,N_): /* (N,N) */
case P(Z_,N_ , N_,N_): /* (N,N) */
case P(Z_,Z_ , N_,N_): /* (N,N) */
/* (N ,N , * ,* ) */
case P(N_,N_ , F_,F_): /* (N,N) */
case P(N_,N_ , F_,I_): /* (N,N) */
case P(N_,N_ , F_,N_): /* (N,N) */
case P(N_,N_ , I_,F_): /* (N,N) */
case P(N_,N_ , I_,I_): /* (N,N) */
case P(N_,N_ , I_,N_): /* (N,N) */
case P(N_,N_ , N_,F_): /* (N,N) */
case P(N_,N_ , N_,I_): /* (N,N) */
case P(N_,N_ , N_,N_): /* (N,N) */
case P(N_,N_ , N_,Z_): /* (N,N) */
case P(N_,N_ , Z_,N_): /* (N,N) */
case P(N_,N_ , Z_,Z_): /* (N,N) */
/* each of the 4 products is a nan */
x = a*c - b*d;
y = b*c + a*d;
break;
default: /* should not happen */
(void)printf("To be done: %#x: a=%g, b=%g, c=%g, d=%g\n",
cp, a, b, c, d);
break;
}/*switch*/
}
#ifdef _Imaginary_I
return x + _Imaginary_I * y;
#else
return cmplx_f(x,y); /* was x + I*y; but has problems with non-imaginary I */
#endif
}
/*****************************************************************
* Divide z / w...; Fred's reference implementation #2
* [(a*c + b*d) + (b*c - a*d)*I] / (c*c + d*d)
*****************************************************************/
float complex cdiv2(float complex z, float complex w){
#pragma STDC FP_CONTRACT OFF
float a, b, c, d;
float x, y;
double aa, bb, cc, dd; /* extra precision and exponent range */
double ac, bd, ad, bc, cc2, dd2, denom, xx, yy;
unsigned int cp;
a = creal(z);
b = cimag(z);
c = creal(w);
d = cimag(w);
#ifdef EXCLUDE_NORMAL
aa = a; /* increase precision and exponent range */
bb = b; /* to avoid spurious overflows, underflows */
cc = c;
dd = d;
ac = aa * cc;
bd = bb * dd;
ad = aa * dd;
bc = bb * cc;
cc2 = cc * cc;
dd2 = dd * dd;
denom = cc2 + dd2;
x = ((ac + bd) / denom); /* may overflow or underflow */
y = ((bc - ad) / denom); /* may overflow or underflow */
if( isnan(x) || isnan(y) ) /* fixup any undeserved NaNs */
#endif
{
/* E.g. P(Z_,F_,I_,N_) indicates a is zero, b finite, c infinite, and d NaN */
cp = P(CLASSIFY(a),CLASSIFY(b),CLASSIFY(c),CLASSIFY(d));
switch (cp) {
/* (* ,* , F ,F ) Following do NOT produce a NaN */
case P(F_,F_ , F_,F_):
case P(F_,I_ , F_,F_):
case P(F_,Z_ , F_,F_):
case P(I_,F_ , F_,F_):
case P(I_,Z_ , F_,F_):
case P(Z_,F_ , F_,F_):
case P(Z_,I_ , F_,F_):
case P(Z_,Z_ , F_,F_):
aa = a; /* increase precision and exponent range */
bb = b; /* to avoid spurious overflows, underflows */
cc = c;
dd = d;
ac = aa * cc;
bd = bb * dd;
ad = aa * dd;
bc = bb * cc;
cc2 = cc * cc;
dd2 = dd * dd;
denom = cc2 + dd2;
x = ((ac + bd) / denom); /* may overflow or underflow */
y = ((bc - ad) / denom); /* may overflow or underflow */
break;
/* (* ,Z , Z ,Z ) */
case P(F_,Z_ , Z_,Z_): /* fin/zero; x-axis/x-axis; (I,0) or (I,N) w/ INV? */
case P(I_,Z_ , Z_,Z_): /* inf/zero; x-axis/x-axis; (I,0) or (I,N) w/ INV? */
#if NODIRECTION /* treat (c,d) as (2n+1)*45 degress => result near 45 deg */
if( O(a)*O(c) == O(b)*O(d) ){ /* |real| > |imag| */
x = a/c; /* inf */
y = (b*c - a*d) / (c*c); /* nan */
}else{ /* |imag| > |real| */
x = (a*c + b*d) / (d*d); /* nan */
y = -(a/d); /* inf */
}
#else /* x-axis / x-axis => x-axis; (a/c, y is signed zero) */
x = a/c;
a = O(a); c = O(c); /* convert to direction vector */
y = (b*c-a*d); /* signed zero */
#endif
break;
case P(N_,Z_ , Z_,Z_): /* "zero" */
#if NODIRECTION /* treat (c,d) as (2n+1)*45 degress => result near 45 deg */
#if (2 == NAN_SIGN) || (1 == NAN_SIGN)
x = a/c; /* "inf" becomes nan */
y = b*c - (a*d); /* nan */
#else
if( O(a)*O(c) == O(b)*O(d) ){ /* |real| > |imag| */
x = a/c; /* "inf" becomes nan */
y = b*c - (a*d); /* nan */
}else{ /* |imag| > |real| */
x = (a*c) + b*d; /* nan */
y = b*c - (a/d); /* "inf" becomes nan */
}
#endif
#else /* x-axis / x-axis => x-axis; (a/c, y is signed zero) */
x = a/c;
#if 2 == NAN_SIGN /* +zero */
y = 0.f;
#elif 1 == NAN_SIGN /* signed zero w/o NaN term */
c = O(c); /* convert to direction vector */
y = (b*c);
#else /* signed zero */
a = O(a); c = O(c); /* convert to direction vector */
y = (b*c-a*d);
#endif
#endif
break;
/* (* ,Z , N ,Z ) */
case P(F_,Z_ , N_,Z_): /* fin / "zero" */
case P(I_,Z_ , N_,Z_): /* inf / "zero" */
case P(Z_,Z_ , N_,Z_): /* zero / "zero"; x-axis/x-axis; (N,0) or (N,N) */
#if NODIRECTION
x = (a*c + b*d) / (c*c);
y = (b*c - a*d) / (c*c); /* raise all invalids */
#else /* x-axis / x-axis => x-axis; (a/c, y is signed zero) */
x = a/c;
#if 2 == NAN_SIGN /* +zero */
y = 0.f;
#elif 1 == NAN_SIGN /* signed zero w/o NaN term */
a = O(a); /* convert to direction vector */
y = -(a*d);
#else /* signed zero */
a = O(a); c = O(c); /* convert to direction vector */
y = (b*c-a*d);
#endif
#endif
break;
case P(N_,Z_ , N_,Z_): /* "zero" / "zero"; x-axis/x-axis; (N,0) or (N,N) */
#if NODIRECTION
x = (a*c + b*d) / (c*c);
y = (b*c - a*d) / (c*c); /* raise all invalids */
#else /* x-axis / x-axis => x-axis; (a/c, y is signed zero) */
x = a/c;
#if (2 == NAN_SIGN) || (1 == NAN_SIGN)/* +zero */
y = 0.f;
#else /* signed zero */
a = O(a); c = O(c); /* convert to direction vector */
y = (b*c-a*d);
#endif
#endif
break;
case P(I_,F_ , I_,F_):
/* (i ,Z , I ,F ) */
case P(I_,Z_ , I_,F_):
x = a/c;
#if CLOSE
y = (b*c - a*d) / (c*c); /* raise all invalids */
#else /* x-axis / x-axis => x-axis; (a/c, y is signed zero) */
b = Z(b); d = Z(d); /* box the inf's */
a = O(a); c = O(c); /* convert to direction vector */
y = (b*c-a*d); /* signed zero */
#endif
break;
case P(N_,Z_ , I_,F_): /* "inf" / inf */
x = a/c;
#if CLOSE
y = (b*c - a*d) / (c*c); /* raise all invalids */
#else /* x-axis / x-axis => x-axis; (a/c, y is signed zero) */
#if 2 == NAN_SIGN /* +zero */
y = 0.f;
#elif 1 == NAN_SIGN /* signed zero w/o NaN term */
c = O(c); /* convert to direction vector */
y = (b*c);
#else /* signed zero */
b = Z(b); d = Z(d); /* box the inf's */
a = O(a); c = O(c); /* convert to direction vector */
y = (b*c-a*d);
#endif
#endif
break;
/* (* ,Z , I ,F ) */
case P(F_,Z_ , I_,F_): /* fin / inf */
/* x-axis / x-axis => x-axis; (a/c, y is signed zero) */
x = a/c;
y = -(0.f*a)*d; /* (-a*d)/(c*c) */
break;
case P(Z_,Z_ , I_,F_): /* zero / inf */
#if NODIRECTION /* treat as any/x-axis */
c = O(c); /* avoid 0*inf */
d = Z(d); /* avoid a*d overflow or underflow */
x = 0.f*(a*c+b*d); /* signed zero */
#else /* x-axis / x-axis => x-axis; (a/c, y is signed zero) */
x = a/c;
c = O(c); /* avoid 0*inf */
d = Z(d); /* avoid a*d overflow or underflow */
#endif
y = 0.f*(b*c-a*d); /* signed zero */
break;
/* (Z ,Z , * ,Z ) */
case P(Z_,Z_ , F_,Z_): /* no NaN */
case P(Z_,Z_ , I_,Z_):
c = O(c); /* avoid 0*inf */
#if NODIRECTION
x = (a*c+b*d); /* signed zero */
#else /* x-axis / x-axis => x-axis; (a/c, y is signed zero) */
x = a/c;
#endif
y = (b*c-a*d); /* signed zero */
break;
/* (* ,Z , * ,Z ) */
case P(F_,Z_ , F_,Z_): /* no NaN */
case P(F_,Z_ , I_,Z_):
case P(I_,Z_ , F_,Z_):
case P(I_,Z_ , I_,Z_):
/* x-axis / x-axis => x-axis; (a/c, y is signed zero) */
x = a/c;
a = O(a); c = O(c); /* convert to direction vector */
y = (b*c-a*d); /* signed zero */
break;
case P(N_,Z_ , F_,Z_): /* "zero" */
case P(N_,Z_ , I_,Z_): /* "zero" */
x = a/c;
#if 2 == NAN_SIGN /* +zero */
y = 0.f;
#elif 1 == NAN_SIGN /* signed zero w/o NaN term */
c = O(c); /* convert to direction vector */
y = (b*c);
#else /* signed zero */
a = O(a); c = O(c); /* convert to direction vector */
y = (b*c-a*d);
#endif
break;
/* (* ,Z , Z ,N ) */
case P(F_,Z_ , Z_,N_): /* "zero" */
case P(I_,Z_ , Z_,N_): /* "zero" */
case P(Z_,Z_ , Z_,N_): /* zero / "zero"; x-axis/y-axis */
#if NODIRECTION
x = (a*c + b*d) / (d*d);
y = (b*c - a*d) / (d*d);
#else /* x-axis / y-axis => y-axis; (x is signed zero, -a/d) */
y = -(a/d);
#if 2 == NAN_SIGN /* +zero */
x = 0.f;
#elif 1 == NAN_SIGN /* signed zero w/o NaN term */
a = O(a); /* convert to direction vector */
x = (a*c);
#else /* signed zero */
a = O(a); d = O(d); /* convert to direction vector */
x = (a*c+b*d);
#endif
#endif
break;
case P(N_,Z_ , Z_,N_): /* "zero" */
#if NODIRECTION
x = (a*c + b*d) / (d*d);
y = (b*c - a*d) / (d*d);
#else /* x-axis / y-axis => y-axis; (x is signed zero, -a/d) */
y = -(a/d);
#if (2 == NAN_SIGN) || (1 == NAN_SIGN)/* +zero */
x = 0.f;
#else /* signed zero */
a = O(a); d = O(d); /* convert to direction vector */
x = (a*c+b*d);
#endif
#endif
break;
case P(I_,F_ , F_,I_):
/* (i ,Z , F ,I ) */
case P(I_,Z_ , F_,I_):
case P(N_,Z_ , F_,I_): /* "inf" / inf */
#if CLOSE
x = (a*c + b*d)/(d*d); /* raise all invalids */
y = (b*c - a*d)/(d*d);
#else /* x-axis / y-axis => y-axis; (x is signed zero, -a/d) */
y = -(a/d);
#if 2 == NAN_SIGN /* +zero */
x = 0.f;
#elif 1 == NAN_SIGN /* signed zero w/o NaN term */
b = Z(b); /* box the inf's */
d = O(d); /* convert to direction vector */
x = (b*d);
#else /* signed zero */
b = Z(b); c = Z(c); /* box the inf's */
a = O(a); d = O(d); /* convert to direction vector */
x = (a*c+b*d);
#endif
#endif
break;
/* (* ,Z , F ,I ) */
case P(F_,Z_ , F_,I_):
/* x-axis / y-axis => y-axis; (x is signed zero, -a/d) */
x = (0.f*a)*c; /* as if (a*c)/(d*d) */
y = -(a/d);
break;
case P(Z_,Z_ , F_,I_):
#if NODIRECTION
d = O(d); /* avoid 0*inf */
y = 0.f*(b*c-a*d); /* signed zero */
#else /* x-axis / y-axis => y-axis; (x is signed zero, -a/d) */
y = -(a/d); /* signed zero */
d = O(d); /* avoid 0*inf */
#endif
x = 0.f*(a*c+b*d); /* signed zero */
break;
/* (Z ,Z , Z ,* ) */
case P(Z_,Z_ , Z_,F_): /* no NaN */
case P(Z_,Z_ , Z_,I_):
#if NODIRECTION
d = O(d); /* avoid 0*inf */
y = (b*c-a*d); /* signed zero */
#else /* x-axis / y-axis => y-axis; (x is signed zero, -a/d) */
y = -(a/d); /* signed zero */
d = O(d); /* avoid 0*inf */
#endif
x = (a*c+b*d); /* signed zero */
break;
/* (* ,Z , Z ,* ) */
case P(F_,Z_ , Z_,F_): /* no NaN */
case P(F_,Z_ , Z_,I_):
case P(I_,Z_ , Z_,F_):
case P(I_,Z_ , Z_,I_):
/* x-axis / y-axis => y-axis; (x is signed zero, -a/d) */
y = -(a/d);
a = O(a); d = O(d); /* convert to direction vector */
x = 0.f*(a*c+b*d); /* signed zero */
break;
case P(N_,Z_ , Z_,F_): /* "zero" */
case P(N_,Z_ , Z_,I_): /* "zero" */
/* x-axis / y-axis => y-axis; (x is signed zero, -a/d) */
y = -(a/d);
#if 2 == NAN_SIGN /* +zero */
x = 0.f;
#elif 1 == NAN_SIGN /* signed zero w/o NaN term */
d = O(d); /* convert to direction vector */
x = (b*d);
#else /* signed zero */
a = O(a); d = O(d); /* convert to direction vector */
x = (a*c+b*d);
#endif
break;
/* (Z ,* , Z ,Z ) */
case P(Z_,F_ , Z_,Z_): /* fin/zero; y-axis/x-axis; (0,I) or (N,I) w/ INV? */
case P(Z_,I_ , Z_,Z_): /* inf/zero; y-axis/x-axis; (0,I) or (N,I) w/ INV? */
#if NODIRECTION /* treat (c,d) as (2n+1)*45 degress => result near 45 deg */
if( O(a)*O(c) == O(b)*O(d) ){ /* |real| > |imag| */
x = b/d; /* inf */
y = (b*c - a*d) / (d*d); /* nan */
}else{ /* |imag| > |real| */
x = (a*c + b*d) / (c*c); /* nan */
y = b/c; /* inf */
}
#else /* y-axis / x-axis => y-axis; (x is signed zero, b/c) */
y = b/c;
b = O(b); c = O(c); /* convert to direction vector */
x = (a*c+b*d); /* signed zero */
#endif
break;
case P(Z_,N_ , Z_,Z_): /* "zero" */
#if NODIRECTION /* treat (c,d) as (2n+1)*45 degress => result near 45 deg */
#if (2 == NAN_SIGN) || (1 == NAN_SIGN)
x = b/d; /* inf becomes nan */
y = (b*c) - a*d; /* nan */
#else
if( O(a)*O(c) == O(b)*O(d) ){ /* |real| > |imag| */
x = b/d; /* inf becomes nan */
y = (b*c) - a*d; /* nan */
}else{ /* |imag| > |real| */
x = a*c + (b*d); /* nan */
y = b/c; /* inf becomes nan */
}
#endif
#else /* y-axis / x-axis => y-axis; (x is signed zero, b/c) */
y = b/c;
#if 2 == NAN_SIGN /* +zero */
x = 0.f;
#elif 1 == NAN_SIGN /* signed zero w/o NaN term */
c = O(c); /* convert to direction vector */
x = (a*c);
#else /* signed zero */
b = O(b); c = O(c); /* convert to direction vector */
x = (a*c+b*d);
#endif
#endif
break;
case P(Z_,I_ , I_,F_):
/* (Z ,i , I ,F ) */
case P(F_,I_ , I_,F_):
y = b/c;
#if CLOSE
x = (a*c + b*d) / (c*c); /* raise all invalids */
#else /* y-axis / x-axis => y-axis; (x is signed zero, b/c) */
a = Z(a); d = Z(d); /* box the inf's; fall thru */
b = O(b); c = O(c); /* convert to direction vector */
x = 0.f*(a*c+b*d); /* signed zero */
#endif
break;
/* (Z ,i , I ,F ) */
case P(Z_,N_ , I_,F_): /* "inf" / inf */
y = b/c;
#if CLOSE
x = (a*c + b*d) / (c*c); /* raise all invalids */
#else /* y-axis / x-axis => y-axis; (x is signed zero, b/c) */
#if 2 == NAN_SIGN /* +zero */
x = 0.f;
#elif 1 == NAN_SIGN /* signed zero w/o NaN term*/
c = O(c); /* avoid 0*inf */
x = (a*c);
#else /* signed zero */
c = O(c); b = O(b); /* avoid 0*inf */
x = 0.f*(a*c+b*d);
#endif
#endif
break;
/* (Z ,* , I ,F ) */
case P(Z_,F_ , I_,F_):
/* y-axis / x-axis => y-axis; (x is signed zero, b/c) */
y = b/c;
x = (0.f*b)*d; /* signed zero; as if (b*d)/(c*c) */
break;
/* (Z ,* , N ,Z ) */
case P(Z_,F_ , N_,Z_): /* "zero" */
case P(Z_,I_ , N_,Z_): /* "zero" */
#if NODIRECTION
x = (a*c + b*d) / (c*c); /* raise all invalids */
y = (b*c - a*d) / (c*c);
#else /* y-axis / x-axis => y-axis; (x is signed zero, b/c) */
y = b/c;
#if 2 == NAN_SIGN /* +zero */
x = 0.f;
#elif 1 == NAN_SIGN /* signed zero w/o NaN term */
b = O(b); /* convert to direction vector */
x = (b*d);
#else /* signed zero */
b = O(b); c = O(c); /* convert to direction vector */
x = (a*c+b*d);
#endif
#endif
break;
case P(Z_,N_ , N_,Z_): /* "zero" */
#if NODIRECTION
x = (a*c + b*d) / (c*c); /* raise all invalids */
y = (b*c - a*d) / (c*c);
#else /* y-axis / x-axis => y-axis; (x is signed zero, b/c) */
y = b/c;
#if (2 == NAN_SIGN) || (1 == NAN_SIGN)/* +zero */
x = 0.f;
#else /* signed zero */
b = O(b); c = O(c); /* convert to direction vector */
x = (a*c+b*d);
#endif
#endif
break;
/* (Z ,* , * ,Z ) */
case P(Z_,F_ , F_,Z_): /* no NaN */
case P(Z_,F_ , I_,Z_):
case P(Z_,I_ , F_,Z_):
case P(Z_,I_ , I_,Z_):
/* y-axis / x-axis => y-axis; (x is signed zero, b/c) */
y = b/c;
b = O(b); c = O(c); /* convert to direction vector */
x = (a*c+b*d); /* signed zero */
break;
case P(Z_,N_ , F_,Z_): /* "zero" */
case P(Z_,N_ , I_,Z_): /* "zero" */
/* y-axis / x-axis => y-axis; (x is signed zero, b/c) */
y = b/c;
#if 2 == NAN_SIGN /* +zero */
x = 0.f;
#elif 1 == NAN_SIGN /* signed zero w/o NaN term */
c = O(c); /* convert to direction vector */
x = (a*c);
#else /* signed zero */
b = O(b); c = O(c); /* convert to direction vector */
x = (a*c+b*d);
#endif
break;
case P(F_,I_ , F_,I_):
x = b/d;
#if CLOSE
y = (b*c-a*d)/(d*d); /* raise all invalids */
#else /* y-axis / y-axis => x-axis; (b/d, y is signed zero) */
a = Z(a); c = Z(c); /* box the inf's */
b = O(b); d = O(d); /* convert to direction vector */
y = 0.f*(b*c-a*d); /* signed zero */
#endif
break;
/* (Z ,i , F ,I ) */
case P(Z_,N_ , F_,I_): /* "inf" / inf */
x = b/d;
#if CLOSE
y = (b*c-a*d)/(d*d); /* raise all invalids */
#else /* y-axis / y-axis => x-axis; (b/d, y is signed zero) */
#if 2 == NAN_SIGN /* +zero */
y = 0.f;
#elif 1 == NAN_SIGN /* signed zero w/o NaN term */
d = O(d);
y = -(a*d);
#else /* signed zero */
b = O(b); d = O(d);
y = 0.f*(b*c-a*d);
#endif
#endif
break;
case P(Z_,I_ , F_,I_):
x = b/d;
#if CLOSE
y = (b*c-a*d)/(d*d); /* raise all invalids */
#else /* y-axis / y-axis => x-axis; (b/d, y is signed zero) */
b = O(b); d = O(d);
y = 0.f*(b*c-a*d); /* signed zero */
#endif
break;
/* (Z ,* , F ,I ) */
case P(Z_,F_ , F_,I_):
/* y-axis / y-axis => x-axis; (b/d, y is signed zero) */
x = b/d;
y = (0.f*b)*c; /* signed zero; as if (b*c)/(d*d) */
break;
/* (Z ,* , Z ,N ) */
case P(Z_,F_ , Z_,N_): /* fin / "zero" */
case P(Z_,I_ , Z_,N_): /* inf / "inf" */
#if NODIRECTION
x = (a*c + b*d) / (d*d);
y = (b*c - a*d) / (d*d);
#else /* y-axis / y-axis => x-axis; (b/d, y is signed zero) */
x = b/d;
#if 2 == NAN_SIGN /* +zero */
y = 0.f;
#elif 1 == NAN_SIGN /* signed zero w/o NaN term */
b = O(b); /* convert to direction vector */
y = (b*c);
#else /* signed zero */
b = O(b); d = O(d); /* convert to direction vector */
y = (b*c-a*d);
#endif
#endif
break;
case P(Z_,N_ , Z_,N_): /* "zero" */
#if NODIRECTION
x = (a*c + b*d) / (d*d);
y = (b*c - a*d) / (d*d);
#else /* y-axis / y-axis => x-axis; (b/d, y is signed zero) */
x = b/d;
#if (2 == NAN_SIGN) || (1 == NAN_SIGN)/* +zero */
y = 0.f;
#else /* signed zero */
b = O(b); d = O(d); /* convert to direction vector */
y = (b*c-a*d);
#endif
#endif
break;
/* (Z ,* , Z ,* ) */
case P(Z_,F_ , Z_,F_): /* no NaN */
case P(Z_,F_ , Z_,I_):
case P(Z_,I_ , Z_,F_):
case P(Z_,I_ , Z_,I_):
/* y-axis / y-axis => x-axis; (b/d, y is signed zero) */
x = b/d;
b = O(b); d = O(d); /* convert to direction vector */
y = (b*c-a*d); /* signed zero */
break;
case P(Z_,N_ , Z_,F_): /* "zero" */
case P(Z_,N_ , Z_,I_): /* "zero" */
/* y-axis / y-axis => x-axis; (b/d, y is signed zero) */
x = b/d;
#if 2 == NAN_SIGN /* +zero */
y = 0.f;
#elif 1 == NAN_SIGN /* signed zero w/o NaN term */
d = O(d); /* convert to direction vector */
y = -(a*d);
#else /* signed zero */
b = O(b); d = O(d); /* convert to direction vector */
y = (b*c-a*d);
#endif
break;
/* (f ,f , I ,* ) */
case P(F_,F_ , I_,F_): /* fin / inf => zero */
x = a/c;
y = b/c;
break;
case P(F_,F_ , I_,N_): /* fin / inf => 0 w/o NaN; any/any */
case P(F_,Z_ , I_,N_):
case P(Z_,F_ , I_,N_):
#if 2 == NAN_SIGN /* +zero */
x = 0.f;
y = 0.f;
#elif 1 == NAN_SIGN /* signed zero w/o NaN term */
x = (a/c);
y = (b/c);
#else /* signed zero */
c = O(c); d = Z(d);
x = 0.f*(a*c+b*d);
y = 0.f*(b*c-a*d);
#endif
break;
/* (f ,f , * ,I ) */
case P(F_,F_ , F_,I_): /* fin / inf => zero */
x = b/d;
y = -(a/d);
break;
case P(F_,F_ , N_,I_): /* fin / inf => 0 w/o NaN; any/any */
case P(F_,Z_ , N_,I_):
case P(Z_,F_ , N_,I_):
#if 2 == NAN_SIGN /* +zero */
x = 0.f;
y = 0.f;
#elif 1 == NAN_SIGN /* signed zero w/o NaN term */
x = (b/d);
y = -(a/d);
#else /* signed zero */
c = Z(c); d = O(d);
x = 0.f*(a*c+b*d);
y = 0.f*(b*c-a*d);
#endif
break;
/* (f ,f , I ,I ) */
case P(F_,F_ , I_,I_): /* fin / inf => zero */
case P(F_,Z_ , I_,I_):
case P(Z_,F_ , I_,I_):
c = O(c); d = O(d);
ac = a*c; bd = b*d; bc = b*c; ad = a*d;
x = 0.f*(ac+bd); /* signed zero */
y = 0.f*(bc-ad); /* signed zero */
break;
case P(Z_,Z_ , I_,I_):
#if !NODIRECTION
a = O(a); /* convert to direction vector */
#endif
c = O(c); d = O(d);
x = 0.f*(a*c+b*d); /* signed zero */
y = 0.f*(b*c-a*d); /* signed zero */
break;
case P(I_,I_ , F_,F_): /* inf / fin => inf (may get w/ nan for CLOSE) */
#if CLOSE /* 45deg / any_deg => 45-any; but 45-45 => NaN */
x = a*c + b*d;
y = b*c - a*d;
#else
a = O(a); b = O(b); /* box the inf */
ac = a*c; bd = b*d; bc = b*c; ad = a*d;
xx = ac+bd; yy = bc-ad;
x = (0.0==xx)? xx : INFINITY*xx; /* signed INF */
y = (0.0==yy)? yy : INFINITY*yy; /* signed INF or zero */
#endif
break;
/* (z ,Z , Z ,Z ) */
case P(Z_,Z_ , Z_,Z_):
#if NODIRECTION
x = (a*c+b*d) / (c*c);
y = (b*c-a*d) / (c*c);
#else /* x-axis / x-axis => x-axis; (a/c, y is signed zero) */
x = a/c;
a = O(a); c = O(c); /* convert to direction vector */
y = (b*c-a*d); /* signed zero */
#endif
break;
/* (I ,* , Z ,Z ) */
case P(I_,F_ , Z_,Z_): /* inf/0 is inf; but direction is unknown=>NaN */
#if NODIRECTION /* treat (c,d) as (2n+1)*45 degress => result is near 45 deg */
if( O(a)*O(c) == O(b)*O(d) ){ /* |real| > |imag| */
x = a/c; /* inf */
y = b*c - a*d; /* nan */
}else{ /* |imag| > |real| */
x = a*c + b*d; /* nan */
y = -(a/d); /* inf */
}
#else /* anydeg / x-axis => anydeg */
x = a/c; /* signed INF */
y = copysignf(INFINITY,b*c); /* INF(b/c) but not div-by-0 */
#endif
break;
/* (* ,I , Z ,Z ) */
case P(F_,I_ , Z_,Z_): /* inf/0 is inf; but direction is unknown=>NaN */
#if NODIRECTION /* treat (c,d) as (2n+1)*45 degress => result is near 45 deg */
if( O(a)*O(c) == O(b)*O(d) ){ /* |real| > |imag| */
x = b/d; /* inf */
y = b*c - a*d; /* nan */
}else{ /* |imag| > |real| */
x = a*c + b*d; /* nan */
y = b/c; /* inf */
}
#else /* anydeg / x-axis => anydeg */
x = copysignf(INFINITY,a*c); /* INF(a/c) but not div-by-0 */
y = b/c; /* signed INF */
#endif
break;
/* (* ,* , Z ,Z ) */
case P(F_,F_ , Z_,Z_): /* fin/0 is inf; but direction is unknown=>NaN */
/* if ac has same sign as bd, then real component of result is longer than imag component */
case P(I_,I_ , Z_,Z_): /* treat (a,b) as (2n+1)*45 degress => result on an axis */
#if NODIRECTION /* treat (c,d) as (2n+1)*45 degress */
if( O(a)*O(c) == O(b)*O(d) ){ /* |real| > |imag| */
x = a/c; /* inf */
y = (b*c - a*d) / (c*c); /* nan */
}else{ /* |imag| > |real| */
x = (a*c + b*d) / (c*c); /* nan */
y = b/c; /* inf */
}
#else /* anydeg / x-axis => anydeg */
x = a/c; /* signed INF */
y = b/c; /* signed INF */
#endif
break;
case P(I_,I_ , I_,I_):
x = a/c; /* (NaN,NaN) */
y = b/d;
#if CLOSE
#else
a = O(a); b = O(b); c = O(c); d = O(d);
xx = (a*c+b*d); /* direction of result */
yy = (b*c-a*d);
if(0.0 == xx) x = xx; /* NaN on an axis */
if(0.0 == yy) y = yy;
#endif
break;
/* (Z ,Z , * ,N ) => zero, w/o NaN, w/o exceptions */
case P(Z_,Z_ , F_,N_): /* 0/fin => 0 w/o NaN */
case P(Z_,Z_ , I_,N_): /* 0/inf => 0 w/o NaN */
#if 2 == NAN_SIGN /* +zero */
x = 0.f;
y = 0.f;
#elif 1 == NAN_SIGN /* signed zero w/o NaN term */
x = a/c;
y = b/c;
#else /* signed zero */
d = Z(d);
#if NODIRECTION
c = O(c); /* avoid 0*inf */
x = 0.f*(a*c+b*d); /* signed zero */
#else /* x-axis / x-axis => x-axis; (a/c, y is signed zero) */
x = a/c;
c = O(c); /* avoid 0*inf */
#endif
y = 0.f*(b*c-a*d);
#endif
break;
/* (Z ,Z , N ,* ) => zero, w/o NaN, w/o exceptions */
case P(Z_,Z_ , N_,F_): /* 0/fin => 0 w/o NaN */
case P(Z_,Z_ , N_,I_): /* 0/inf => 0 w/o NaN */
#if 2 == NAN_SIGN /* +zero */
x = 0.f;
y = 0.f;
#elif 1 == NAN_SIGN /* signed zero w/o NaN term */
x = b/d;
y = -(a/d);
#else /* signed zero */
c = Z(c);
#if NODIRECTION
d = O(d); /* avoid 0*inf */
y = 0.f*(b*c-a*d); /* signed zero */
#else /* x-axis / y-axis => y-axis; (x is signed zero, -a/d) */
y = -(a/d); /* as if x-axis / y-axis */
d = O(d); /* avoid 0*inf */
#endif
x = 0.f*(a*c+b*d);
#endif
break;
/* (* ,* , F ,I ) */
case P(F_,N_ , F_,I_): /* "inf" / inf */
#if CLOSE
x = b*d;
y = b*c;
break;
#endif
case P(N_,F_ , F_,I_): /* "inf" / inf */
#if CLOSE
x = a*c ;
y = 0.f - (a*d); /* 0.f needed for sign of NaN */
break;
#endif
/* (* ,* , Z ,* ) */
case P(F_,F_ , Z_,F_): /* no NaN */
case P(F_,F_ , Z_,I_):
case P(F_,I_ , Z_,F_):
case P(F_,I_ , Z_,I_): /* inf / inf */
case P(F_,N_ , Z_,F_):
case P(F_,N_ , Z_,I_): /* "inf" / inf */
case P(I_,F_ , Z_,F_):
case P(I_,F_ , Z_,I_): /* inf / inf */
case P(I_,I_ , Z_,F_):
case P(I_,I_ , Z_,I_): /* inf / inf */
case P(I_,N_ , Z_,F_):
case P(I_,N_ , Z_,I_): /* inf / inf */
case P(N_,F_ , Z_,F_): /* fin+inf / fin; any/y-axis */
case P(N_,F_ , Z_,I_): /* "inf" / inf */
case P(N_,I_ , Z_,F_):
case P(N_,I_ , Z_,I_): /* inf / inf */
case P(N_,N_ , Z_,F_): /* (N,N) */
case P(N_,N_ , Z_,I_): /* (N,N) */
/* any / y-axis => rotate by 90 deg */
x = b/d;
y = -(a/d);
break;
case P(F_,F_ , Z_,N_): /* fin / "zero"; any/y-axis */
case P(F_,I_ , Z_,N_): /* inf / "inf"; "y-axis"/y-axis */
case P(F_,N_ , Z_,N_): /* "inf" / "inf" */
case P(I_,F_ , Z_,N_): /* inf / "inf"; "x-axis"/y-axis */
case P(I_,I_ , Z_,N_): /* inf / "inf" */
case P(I_,N_ , Z_,N_): /* "zero" */
case P(N_,F_ , Z_,N_): /* "inf" / "inf" */
case P(N_,I_ , Z_,N_): /* inf / "inf" */
#if NODIRECTION
x = (a*c + b*d) / (d*d);
y = (b*c - a*d) / (d*d);
#else
x = b/d;
y = -(a/d);
#endif
break;
/* (* ,* , Z ,Z ) */
case P(N_,I_ , Z_,Z_): /* inf/0=>inf w/ nan */
#if NODIRECTION /* w/ nan on y-axis */
x = b/d;
y = b*c - a*d; /* raise invalid */
break;
#endif
case P(N_,F_ , Z_,Z_): /* fin/0=>inf w/ nan */
#if NODIRECTION /* treat (c,d) as (2n+1)*45 degress => result near 45 deg */
#if (2 == NAN_SIGN) || (1 == NAN_SIGN)
x = b/d; /* inf */
y = b*c - a*d; /* nan */
#else
if( O(a)*O(c) == O(b)*O(d) ){ /* |real| > |imag| */
x = b/d; /* inf */
y = b*c - a*d; /* nan */
}else{ /* |imag| > |real| */
x = a*c + b*d; /* nan */
y = b/c; /* inf */
}
#endif
break;
#endif
case P(I_,N_ , Z_,Z_): /* inf/0=>inf w/ nan */
#if NODIRECTION /* w/ nan on y-axis */
x = a/c;
y = b*c - a*d; /* raise invalid */
break;
#endif
case P(F_,N_ , Z_,Z_): /* fin/0=>inf w/ nan */
#if NODIRECTION /* treat (c,d) as (2n+1)*45 degress => result near 45 deg */
#if (2 == NAN_SIGN) || (1 == NAN_SIGN)
x = a/c; /* inf */
y = b*c - a*d; /* nan */
#else
if( O(a)*O(c) == O(b)*O(d) ){ /* |real| > |imag| */
x = a/c; /* inf */
y = b*c - a*d; /* nan */
}else{ /* |imag| > |real| */
x = a*c + b*d; /* nan */
y = -(a/d); /* inf */
}
#endif
#else /* any / x-axis => any */
x = a/c;
y = b/c;
#endif
break;
/* (* ,* , I ,F ) */
case P(F_,N_ , I_,F_): /* "inf" / inf */
#if CLOSE
x = b*d;
y = b*c;
break;
#endif
case P(N_,F_ , I_,F_): /* "inf" / inf */
#if CLOSE
x = a*c ;
y = b*c - (a*d); /* b*c needed for sign of NaN */
break;
#endif
/* (* ,* , * ,Z ) */
case P(F_,F_ , F_,Z_): /* no NaN */
case P(F_,F_ , I_,Z_):
case P(F_,I_ , F_,Z_):
case P(F_,I_ , I_,Z_): /* inf / inf */
case P(F_,N_ , F_,Z_):
case P(F_,N_ , I_,Z_): /* "inf" / inf */
case P(I_,F_ , F_,Z_): /* inf / fin */
case P(I_,F_ , I_,Z_): /* inf / inf */
case P(I_,I_ , F_,Z_):
case P(I_,I_ , I_,Z_): /* inf / inf */
case P(I_,N_ , F_,Z_):
case P(I_,N_ , I_,Z_): /* inf / inf */
case P(N_,F_ , F_,Z_):
case P(N_,F_ , I_,Z_):
case P(N_,I_ , F_,Z_):
case P(N_,I_ , I_,Z_): /* inf / inf */
case P(N_,N_ , F_,Z_): /* (N,N) */
case P(N_,N_ , I_,Z_): /* (N,N) */
/* any / x-axis => any */
x = a/c;
y = b/c;
break;
case P(F_,F_ , N_,Z_): /* fin / "zero"; any/x-axis */
case P(F_,I_ , N_,Z_): /* inf / "inf"; "y-axis"/x-axis */
case P(F_,N_ , N_,Z_): /* "inf" / "inf" */
case P(I_,F_ , N_,Z_): /* inf / "inf"' "x-axis"/x-axis */
case P(I_,I_ , N_,Z_): /* inf / "inf" */
case P(I_,N_ , N_,Z_): /* inf / "inf" */
case P(N_,F_ , N_,Z_): /* "inf" / "inf" */
case P(N_,I_ , N_,Z_): /* inf / "inf" */
#if NODIRECTION
x = (a*c + b*d) / (c*c);
y = (b*c - a*d) / (c*c);
#else
x = a/c;
y = b/c;
#endif
break;
case P(N_,I_ , F_,F_): /* inf / fin => inf w/ nan */
if( fabs(d) < fabs(c) ){
x = a/c; y = b/c; /* div by "x-axis" => (nan,inf) */
}else{
x = b/d; y = b*c-(a/d); /* div by "y-axis" => rotate by 90 deg (inf,nan) */
}
break;
case P(I_,N_ , F_,F_): /* inf / fin => inf w/ nan */
if( fabs(d) <= fabs(c) ){
x = a/c; y = b/c; /* div by "x-axis" => (inf,nan) */
}else{
x = b/d; y = -(a/d); /* div by "y-axis" => rotate by 90 deg (nan,inf) */
}
break;
/* (fz,fz , F ,N ) => zero or NaN */
case P(F_,F_ , F_,N_): /* f,f/f,z=>(a/c,b/c); f,f/f,i=>(b/d,-a/d) */
case P(F_,Z_ , F_,N_): /* f,z/f,z=>(a/c,0); f,z/f,i=>(0,-a/d); fz/ff=>(ac,-ad)/denom */
case P(Z_,F_ , F_,N_): /* z,f/f,z=>(0,b/c); z,f/f,i=>(b/d,0); zf/ff=>(bd,bc)/denom */
/* Assume d is zero and see if rounded result is zero. */
dd = c*d;
x = y = 1.f;
ac = a; ac /= c;
bc = b; bc /= c;
#if NODIRECTION
#if (2 == NAN_SIGN) || (1 == NAN_SIGN)
#else /* honor sign of NaN */
d = Z(d);
ac += b*d;
bc -= a*d;
#endif
#endif
if( (fabs(ac) < (FLT_MIN*FLT_EPSILON))
&& (fabs(bc) < (FLT_MIN*FLT_EPSILON)) ){ /* both may underflow */
x = ac; /* depends upon rounding direction: either 0 or min den */
y = bc; /* may raise spurious underflow */
}
if(x || y) x = y = dd; /* non-zero => need NaN */
break;
/* (fz,fz , N ,F ) => zero or NaN */
case P(F_,F_ , N_,F_): /* f,f/z,f=>(b/d,-a/d); f,f/i,f=>(a/c,b/c) */
case P(F_,Z_ , N_,F_): /* f,z/z,f=>(0,-a/d); f,z/i,f=>(a/c,0); fz/ff=>(ac,-ad)/denom */
case P(Z_,F_ , N_,F_): /* z,f/z,f=>(b/d,0); z,f/i,f=>(0,b/c); zf/ff=>(bd,bc)/denom */
/* Assume c is zero and see if rounded result is zero. */
cc = c*d;
x = y = 1.f;
bd = b; bd /= d;
ad = -a; ad /= d;
#if NODIRECTION
#if (2 == NAN_SIGN) || (1 == NAN_SIGN)
#else /* honor sign of NaN */
c = Z(c);
bd += a*c;
ad += b*c;
#endif
#endif
if( (fabs(bd) < (FLT_MIN*FLT_EPSILON))
&& (fabs(ad) < (FLT_MIN*FLT_EPSILON)) ){ /* both may underflow */
x = bd; /* depends upon rounding direction: either 0 or min den */
y = ad; /* may raise spurious underflow */
}
if(x || y) x = y = cc; /* non-zero => need NaN */
break;
case P(F_,N_ , F_,F_): /* "fin" / fin; x-axis/any */
x = y = 0.f;
cc = c; /* increase precision and exponent range */
dd = d; /* to avoid spurious underflows */
aa = a; /* treat 'b' as zero to see if overflow to infinity */
ac = aa * cc;
ad = aa * dd;
cc2 = cc * cc;
dd2 = dd * dd;
denom = cc2 + dd2;
xx = ac / denom; /* may raise spurious inexact */
yy = -ad / denom;
if( FLT_MAX < fabs(xx) ) x = xx; /* may overflow */
if( !isinf(x) ) x = b*d; /* inf stronger than nan */
if( FLT_MAX < fabs(yy) ) y = yy; /* may overflow */
if( !isinf(y) ) y = b*c; /* inf stronger than nan */
if( isinf(x) && isinf(y) ){
if( fabs(c) < fabs(d) ){ /* "y-axis" */
x = b*d; /* (nan,inf) */
}else{ /* "x-axis" */
y = b*c; /* (inf,nan) */
}
}
break;
case P(N_,F_ , F_,F_): /* "fin" / fin; y-axis/any */
x = y = 0.f;
cc = c; /* increase precision and exponent range */
dd = d; /* to avoid spurious underflows */
bb = b; /* treat 'a' as zero to see if overflow to infinity */
bd = bb * dd;
bc = bb * cc;
cc2 = cc * cc;
dd2 = dd * dd;
denom = cc2 + dd2;
xx = bd / denom; /* may raise spurious inexact */
yy = bc / denom;
if( FLT_MAX < fabs(xx) ) x = xx; /* may overflow */
if( !isinf(x) ) x = a*c; /* inf stronger than nan */
if( FLT_MAX < fabs(yy) ) y = yy; /* may overflow */
if( !isinf(y) ) y = bc - (a*d); /* inf stronger than nan */
if( isinf(x) && isinf(y) ){
if( fabs(c) <= fabs(d) ){ /* "y-axis" */
y = bc - (a*d); /* (inf,nan) */
}else{ /* "x-axis" */
x = a*c; /* (nan,inf) */
}
}
break;
/* (z ,n , F ,F ) */
case P(N_,Z_ , F_,F_): /* "zero" / fin; x-axis/any */
x = a*c;
y = b*c - (a*d); /* b*c affects sign of NaN */
break;
case P(Z_,N_ , F_,F_): /* "zero" / fin; y-axis/any */
x = b*d;
y = b*c;
break;
/* ("inf" , "inf") */
/* (* ,N , * ,N ) */
case P(F_,N_ , F_,N_): /* "inf" / "inf" */
case P(Z_,N_ , F_,N_): /* "inf" / "inf" */
case P(F_,N_ , I_,N_): /* "inf" / inf */
case P(I_,N_ , F_,N_): /* inf / "inf" */
case P(I_,N_ , I_,N_): /* inf / inf */
denom = d*d;
x = ( b*d) / denom;
y = (b*c - a*d) / denom;
break;
/* (N ,* , * ,N ) */
case P(N_,F_ , F_,N_): /* "inf" / "inf" */
case P(N_,Z_ , F_,N_): /* "inf" / "inf" */
case P(N_,F_ , I_,N_): /* "inf" / inf */
case P(N_,I_ , F_,N_): /* inf / "inf" */
case P(N_,I_ , I_,N_): /* inf / inf */
denom = d*d;
x = (a*c + b*d ) / denom;
y = (0.f -(a*d)) / denom; /* 0.f needed for sign of NaN */
break;
/* (* ,N , N ,* ) */
case P(F_,N_ , N_,F_): /* "inf" / "inf" */
case P(Z_,N_ , N_,F_): /* "inf" / "inf" */
case P(F_,N_ , N_,I_): /* "inf" / inf */
case P(I_,N_ , N_,F_): /* inf / "inf" */
case P(I_,N_ , N_,I_): /* inf / inf */
denom = c*c;
x = (a*c + b*d) / denom;
y = (b*c ) / denom;
break;
/* (N ,* , N ,* ) */
case P(N_,F_ , N_,F_): /* "inf" / "inf" */
case P(N_,Z_ , N_,F_): /* "inf" / "inf" */
case P(N_,F_ , N_,I_): /* "inf" / inf */
case P(N_,I_ , N_,F_): /* inf / "inf" */
case P(N_,I_ , N_,I_): /* inf / inf */
denom = c*c;
x = (a*c ) / denom;
y = (b*c - a*d) / denom;
break;
/* (* ,* , * ,N ) */
case P(F_,I_ , F_,N_): /* inf / "inf" */
case P(I_,F_ , F_,N_): /* inf / "inf" */
case P(I_,I_ , F_,N_): /* inf / "inf" */
case P(I_,Z_ , F_,N_): /* inf / "inf" */
case P(Z_,I_ , F_,N_): /* inf / "inf" */
case P(F_,I_ , I_,N_): /* inf / inf */
case P(I_,F_ , I_,N_): /* inf / inf */
case P(I_,I_ , I_,N_): /* inf / inf */
denom = d*d;
x = ( b*d ) / denom;
y = (0.f-(a*d)) / denom; /* 0.f needed for sign of NaN */
break;
/* (* ,* , N ,* ) */
case P(F_,I_ , N_,F_): /* inf / "inf" */
case P(I_,F_ , N_,F_): /* inf / "inf" */
case P(I_,I_ , N_,F_): /* inf / "inf" */
case P(I_,Z_ , N_,F_): /* inf / "inf" */
case P(Z_,I_ , N_,F_): /* inf / "inf" */
case P(F_,I_ , N_,I_): /* inf / inf */
case P(I_,F_ , N_,I_): /* inf / inf */
case P(I_,I_ , N_,I_): /* inf / inf */
denom = c*c;
x = (a*c) / denom;
y = (b*c) / denom;
break;
/* (* ,n , I ,I ) */
case P(F_,N_ , I_,I_): /* "inf" / inf */
case P(I_,N_ , I_,I_): /* inf / inf */
case P(I_,Z_ , I_,I_): /* inf / inf */
x = b*d;
y = b*c;
break;
/* (n ,* , I ,I ) */
case P(N_,F_ , I_,I_): /* "inf" / inf */
case P(N_,I_ , I_,I_): /* inf / inf */
case P(Z_,I_ , I_,I_): /* inf / inf */
x = a*c ;
y = b*c - (a*d); /* b*c needed for sign of NaN */
break;
/* (* ,* , * ,in) */
case P(I_,I_ , F_,I_): /* inf / inf */
case P(I_,N_ , F_,I_): /* inf / inf */
case P(N_,I_ , F_,I_): /* inf / inf */
case P(I_,Z_ , I_,N_): /* inf / inf */
case P(Z_,I_ , I_,N_): /* inf / inf */
denom = d*d;
x = (a*c + b*d) / denom;
y = (b*c - a*d) / denom;
break;
/* (* ,* , in,* ) */
case P(I_,I_ , I_,F_): /* inf / inf */
case P(I_,N_ , I_,F_): /* inf / inf */
case P(N_,I_ , I_,F_): /* inf / inf */
case P(I_,Z_ , N_,I_): /* inf / inf */
case P(Z_,I_ , N_,I_): /* inf / inf */
case P(F_,I_ , I_,I_): /* inf / inf */
case P(I_,F_ , I_,I_): /* inf / inf */
denom = c*c;
x = (a*c + b*d) / denom;
y = (b*c - a*d) / denom;
break;
case P(N_,Z_ , I_,I_): /* "inf" / inf */
case P(N_,Z_ , I_,N_): /* "inf" / inf */
case P(N_,Z_ , N_,I_): /* "inf" / inf */
case P(Z_,N_ , I_,I_): /* "inf" / inf */
case P(Z_,N_ , I_,N_): /* "inf" / inf */
case P(Z_,N_ , N_,I_): /* "inf" / inf */
denom = c*c + d*d;
x = (a*c + b*d) / denom;
y = (b*c - a*d) / denom;
break;
/* (* ,* , N ,N ) */
case P(F_,F_ , N_,N_): /* (N,N) */
case P(F_,I_ , N_,N_): /* (N,N) */
case P(F_,N_ , N_,N_): /* (N,N) */
case P(F_,Z_ , N_,N_): /* (N,N) */
case P(I_,F_ , N_,N_): /* (N,N) */
case P(I_,I_ , N_,N_): /* (N,N) */
case P(I_,N_ , N_,N_): /* (N,N) */
case P(I_,Z_ , N_,N_): /* (N,N) */
case P(N_,F_ , N_,N_): /* (N,N) */
case P(N_,I_ , N_,N_): /* (N,N) */
case P(N_,Z_ , N_,N_): /* (N,N) */
case P(Z_,F_ , N_,N_): /* (N,N) */
case P(Z_,I_ , N_,N_): /* (N,N) */
case P(Z_,N_ , N_,N_): /* (N,N) */
case P(Z_,Z_ , N_,N_): /* (N,N) */
/* (N ,N , * ,* ) */
case P(N_,N_ , F_,F_): /* (N,N) */
case P(N_,N_ , F_,I_): /* (N,N) */
case P(N_,N_ , F_,N_): /* (N,N) */
case P(N_,N_ , I_,F_): /* (N,N) */
case P(N_,N_ , I_,I_): /* (N,N) */
case P(N_,N_ , I_,N_): /* (N,N) */
case P(N_,N_ , N_,F_): /* (N,N) */
case P(N_,N_ , N_,I_): /* (N,N) */
case P(N_,N_ , N_,N_): /* (N,N) */
case P(N_,N_ , N_,Z_): /* (N,N) */
case P(N_,N_ , Z_,N_): /* (N,N) */
case P(N_,N_ , Z_,Z_): /* (N,N) */
x = a*c + b*d;
y = b*c - a*d;
break;
default: /* should not happen */
(void)printf("To be done: %#x: a=%g, b=%g, c=%g, d=%g\n",
cp, a, b, c, d);
break;
}/*switch*/
}
#ifdef _Imaginary_I
return x + _Imaginary_I * y;
#else
return cmplx_f(x,y); /* was x + I*y; but has problems with non-imaginary I */
#endif
}