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 * _F_cplx_div_rx(a, w) returns a / w with infinities handled according 2N/A * If a and w are both finite and w is nonzero, _F_cplx_div_rx(a, 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 = (a * c) 2N/A * / r and y = (-a * d) / r with r = c * c + d * d. This implementa- 2N/A * tion computes intermediate results in double precision to avoid 2N/A * premature underflow or overflow. 2N/A * If a is neither NaN nor zero and w is zero, or if a is infinite 2N/A * and w is finite and nonzero, _F_cplx_div_rx delivers an infinite 2N/A * result. If a is finite and w is infinite, _F_cplx_div_rx delivers 2N/A * If a and w are both zero or both infinite, or if either a or w is 2N/A * NaN, _F_cplx_div_rx 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 * 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 = (
double)c * c + (
double)d * d;
2N/A /* w is zero; multiply a by 1/Re(w) - I * Im(w) */ 2N/A if (i) {
/* a is infinite */ 2N/A ((
float *)&v)[0] = a * c;
2N/A ((
float *)&v)[
1] = (a ==
0.0f)? a * c : -a * d;
2N/A if (x != x || y != y) {
2N/A * x or y is NaN, so a 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 * The following is equivalent to 2N/A ((
float *)&v)[0] = (
float)x;
2N/A ((
float *)&v)[
1] = (
float)y;