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 2004 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 * _F_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, _F_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 computes intermediate results in extended precision to avoid 2N/A * 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, _F_cplx_div_ix delivers an infinite 2N/A * result. If b is finite and w is infinite, _F_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, _F_cplx_div_ix delivers NaN + I * NaN. C99 doesn't specify 2N/A * This implementation can raise spurious invalid operation, inexact, 2N/A * and division-by-zero exceptions. C99 allows this. 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 <<
1) -
0xff000000) == 0)? (
1 | (
xx.i >>
31)) : 0);
2N/A long double r, x, y;
2N/A * The following is equivalent to 2N/A * c = crealf(w); d = cimagf(w); 2N/A c = ((
float *)&w)[0];
2N/A d = ((
float *)&w)[
1];
2N/A r = (
long double)c * c + (
long double)d * d;
2N/A /* w is zero; multiply b by 1/Re(w) - I * Im(w) */ 2N/A if (j) {
/* b is infinite */ 2N/A ((
float *)&v)[0] = (b ==
0.0f)? b * c : b * d;
2N/A ((
float *)&v)[
1] = b * c;
2N/A r = (
long double)b / r;
2N/A x = (
long double)d * r;
2N/A y = (
long double)c * r;
2N/A if (x != x || y != y) {
2N/A * x or y is NaN, so b and w can't both be finite and 2N/A * nonzero. Since we handled the case w = 0 above, the 2N/A * only case to check here is when w is infinite. 2N/A if (i | j) {
/* w is infinite */ 2N/A c = (
cc.i < 0)? -
0.0f :
0.0f;
2N/A d = (
dd.i < 0)? -
0.0f :
0.0f;
2N/A x = (
long double)d * b;
2N/A y = (
long double)c * b;
2N/A * The following is equivalent to 2N/A ((
float *)&v)[0] = (
float)x;
2N/A ((
float *)&v)[
1] = (
float)y;