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_ix(b, w) returns (I * b) / w with infinities handled 2N/A * If b and w are both finite and w is nonzero, _D_cplx_div_ix(b, w) 2N/A * delivers the complex quotient q according to the usual formula: 2N/A * let c = Re(w), and d = Im(w); then q = x + I * y where x = (b * d) 2N/A * / r and y = (b * c) / r with r = c * c + d * d. This implementa- 2N/A * tion scales to avoid premature underflow or overflow. 2N/A * If b is neither NaN nor zero and w is zero, or if b is infinite 2N/A * and w is finite and nonzero, _D_cplx_div_ix delivers an infinite 2N/A * result. If b is finite and w is infinite, _D_cplx_div_ix delivers 2N/A * If b and w are both zero or both infinite, or if either b or w is 2N/A * NaN, _D_cplx_div_ix delivers NaN + I * NaN. C99 doesn't specify 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 * scl[i].d = 2^(250*(4-i)) for i = 0, ..., 9 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 * The following is equivalent to 2N/A * c = creal(w); d = cimag(w); 2N/A c = ((
double *)&w)[0];
2N/A d = ((
double *)&w)[
1];
2N/A /* extract high-order words to estimate |b| 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 c = (
cc.i[0] < 0)? -
0.0 :
0.0;
2N/A d = (
dd.i[0] < 0)? -
0.0 :
0.0;
2N/A }
else /* w is nan */ 2N/A ((
double *)&v)[0] = b * d;
2N/A ((
double *)&v)[
1] = b * c;
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 b by 1/Re(w) - I * Im(w) */ 2N/A if (j) {
/* b is infinite */ 2N/A ((
double *)&v)[0] = (b ==
0.0)? b * c : b * d;
2N/A ((
double *)&v)[
1] = b * c;
2N/A if (
hb >=
0x7ff00000) {
/* a is inf or nan */ 2N/A ((
double *)&v)[0] = b * d;
2N/A ((
double *)&v)[
1] = b * c;
2N/A * Compute the real and imaginary parts of the quotient, 2N/A * scaling to avoid overflow or underflow. 2N/A /* compensate for scaling */ 2N/A ((
double *)&v)[0] = d;
2N/A ((
double *)&v)[
1] = c;