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 /* FLT_MAX_EXP */ #include /* copysign(), fabs(), fmax(), logb(), isinf(), * nextafter(), scalbnl(), isnan(), isfinite() */ #include /* 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 /* FLT_MAX_EXP */ #include /* copysign(), fabs(), fmax(), logb(), isinf(), * nextafter(), scalbnl(), isnan(), isfinite() */ #include /* 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 /* FLT_MAX */ #include /* fabs, isinf, isnan, isfinite */ #include /* 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 }