2N/A * The contents of this file are subject to the terms of the 2N/A * Common Development and Distribution License, Version 1.0 only 2N/A * (the "License"). You may not use this file except in compliance 2N/A * See the License for the specific language governing permissions 2N/A * and limitations under the License. 2N/A * When distributing Covered Code, include this CDDL HEADER in each 2N/A * If applicable, add the following below this CDDL HEADER, with the 2N/A * fields enclosed by brackets "[]" replaced with your own identifying 2N/A * information: Portions Copyright [yyyy] [name of copyright owner] 2N/A * Copyright 2003 Sun Microsystems, Inc. All rights reserved. 2N/A * Use is subject to license terms. 2N/A#
pragma ident "%Z%%M% %I% %E% SMI" 2N/A * _D_cplx_div(z, w) returns z / w with infinities handled according 2N/A * If z and w are both finite and w is nonzero, _D_cplx_div(z, w) 2N/A * delivers the complex quotient q according to the usual formula: 2N/A * let a = Re(z), b = Im(z), c = Re(w), and d = Im(w); then q = x + 2N/A * I * y where x = (a * c + b * d) / r and y = (b * c - a * d) / r 2N/A * with r = c * c + d * d. This implementation scales to avoid 2N/A * premature underflow or overflow. 2N/A * If z is neither NaN nor zero and w is zero, or if z is infinite 2N/A * and w is finite and nonzero, _D_cplx_div delivers an infinite 2N/A * result. If z is finite and w is infinite, _D_cplx_div delivers 2N/A * If z and w are both zero or both infinite, or if either z or w is 2N/A * a complex NaN, _D_cplx_div delivers NaN + I * NaN. C99 doesn't 2N/A * specify these cases. 2N/A * This implementation can raise spurious underflow, overflow, in- 2N/A * valid operation, inexact, and division-by-zero exceptions. C99 2N/A * Warning: Do not attempt to "optimize" this code by removing multi- 2N/A * plications by zero. 2N/A * Return +1 if x is +Inf, -1 if x is -Inf, and 0 otherwise 2N/A return (((((
xx.i[0] <<
1) -
0xffe00000) |
xx.i[
1]) == 0)?
2N/A (
1 | (
xx.i[0] >>
31)) : 0);
2N/A double a, b, c, d, r;
2N/A * The following is equivalent to 2N/A * a = creal(z); b = cimag(z); 2N/A * c = creal(w); d = cimag(w); 2N/A a = ((
double *)&z)[0];
2N/A b = ((
double *)&z)[
1];
2N/A c = ((
double *)&w)[0];
2N/A d = ((
double *)&w)[
1];
2N/A /* extract high-order words to estimate |z| and |w| */ 2N/A /* check for special cases */ 2N/A if (
hw >=
0x7ff00000) {
/* w is inf or nan */ 2N/A if (i | j) {
/* w is infinite */ 2N/A * "factor out" infinity, being careful to preserve 2N/A * signs of finite values 2N/A c = i? i : ((
cc.i[0] < 0)? -
0.0 :
0.0);
2N/A d = j? j : ((
dd.i[0] < 0)? -
0.0 :
0.0);
2N/A /* scale to avoid overflow below */ 2N/A ((
double *)&v)[0] = (a * c + b * d) * r;
2N/A ((
double *)&v)[
1] = (b * c - a * d) * r;
2N/A * This nonsense is needed to work around some SPARC 2N/A * implementations of nonstandard mode; if both parts 2N/A * of w are subnormal, multiply them by one to force 2N/A * them to be flushed to zero when nonstandard mode 2N/A * is enabled. Sheesh. 2N/A /* w is zero; multiply z by 1/Re(w) - I * Im(w) */ 2N/A if (i | j) {
/* z is infinite */ 2N/A ((
double *)&v)[0] = a * c + b * d;
2N/A ((
double *)&v)[
1] = b * c - a * d;
2N/A if (
hz >=
0x7ff00000) {
/* z is inf or nan */ 2N/A if (i | j) {
/* z is infinite */ 2N/A ((
double *)&v)[0] = (a * c + b * d) * r;
2N/A ((
double *)&v)[
1] = (b * c - a * d) * r;
2N/A * Scale c and d to compute 1/|w|^2 and the real and imaginary 2N/A * parts of the quotient. 2N/A * Note that for any s, if we let c' = sc, d' = sd, c'' = sc', 2N/A * and d'' = sd', then 2N/A * (ac'' + bd'') / (c'^2 + d'^2) = (ac + bd) / (c^2 + d^2) 2N/A * and similarly for the imaginary part of the quotient. We want 2N/A * to choose s such that (i) r := 1/(c'^2 + d'^2) can be computed 2N/A * without overflow or harmful underflow, and (ii) (ac'' + bd'') 2N/A * and (bc'' - ad'') can be computed without spurious overflow or 2N/A * harmful underflow. To avoid unnecessary rounding, we restrict 2N/A * s to a power of two. 2N/A * To satisfy (i), we need to choose s such that max(|c'|,|d'|) 2N/A * is not too far from one. To satisfy (ii), we need to choose 2N/A * s such that max(|c''|,|d''|) is also not too far from one. 2N/A * There is some leeway in our choice, but to keep the logic 2N/A * from getting overly complicated, we simply attempt to roughly 2N/A * balance these constraints by choosing s so as to make r about 2N/A * the same size as max(|c''|,|d''|). This corresponds to choos- 2N/A * ing s to be a power of two near |w|^(-3/4). 2N/A * Regarding overflow, observe that if max(|c''|,|d''|) <= 1/2, 2N/A * then the computation of (ac'' + bd'') and (bc'' - ad'') can- 2N/A * not overflow; otherwise, the computation of either of these 2N/A * values can only incur overflow if the true result would be 2N/A * within a factor of two of the overflow threshold. In other 2N/A * words, if we bias the choice of s such that at least one of 2N/A * max(|c''|,|d''|) <= 1/2 or r >= 2 2N/A * always holds, then no undeserved overflow can occur. 2N/A * To cope with underflow, note that if r < 2^-53, then any 2N/A * intermediate results that underflow are insignificant; either 2N/A * they will be added to normal results, rendering the under- 2N/A * flow no worse than ordinary roundoff, or they will contribute 2N/A * to a final result that is smaller than the smallest subnormal 2N/A * number. Therefore, we need only modify the preceding logic 2N/A * when z is very small and w is not too far from one. In that 2N/A * case, we can reduce the effect of any intermediate underflow 2N/A * to no worse than ordinary roundoff error by choosing s so as 2N/A * to make max(|c''|,|d''|) large enough that at least one of 2N/A * (ac'' + bd'') or (bc'' - ad'') is normal. 2N/A hs = (((
hw >>
2) -
hw) +
0x6fd7ffff) &
0xfff00000;
2N/A if (
hz <
0x07200000) {
/* |z| < 2^-909 */ 2N/A if (((
hw -
0x32800000) | (
0x47100000 -
hw)) >= 0)
2N/A hs = (((
0x47100000 -
hw) >>
1) &
0xfff00000)
2N/A r =
1.0 / (c * c + d * d);
2N/A ((
double *)&v)[0] = (a * c + b * d) * r;
2N/A ((
double *)&v)[
1] = (b * c - a * d) * r;