1 /* 2 * CDDL HEADER START 3 * 4 * The contents of this file are subject to the terms of the 5 * Common Development and Distribution License, Version 1.0 only 6 * (the "License"). You may not use this file except in compliance 7 * with the License. 8 * 9 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE 10 * or http://www.opensolaris.org/os/licensing. 11 * See the License for the specific language governing permissions 12 * and limitations under the License. 13 * 14 * When distributing Covered Code, include this CDDL HEADER in each 15 * file and include the License file at usr/src/OPENSOLARIS.LICENSE. 16 * If applicable, add the following below this CDDL HEADER, with the 17 * fields enclosed by brackets "[]" replaced with your own identifying 18 * information: Portions Copyright [yyyy] [name of copyright owner] 19 * 20 * CDDL HEADER END 21 */ 22 /* 23 * Copyright 2003 Sun Microsystems, Inc. All rights reserved. 24 * Use is subject to license terms. 25 */ 26 27 #pragma ident "%Z%%M% %I% %E% SMI" 28 29 /* 30 * _F_cplx_div_ix(b, w) returns (I * b) / w with infinities handled 31 * according to C99. 32 * 33 * If b and w are both finite and w is nonzero, _F_cplx_div_ix(b, w) 34 * delivers the complex quotient q according to the usual formula: 35 * let c = Re(w), and d = Im(w); then q = x + I * y where x = (b * d) 36 * / r and y = (b * c) / r with r = c * c + d * d. This implementa- 37 * tion computes intermediate results in double precision to avoid 38 * premature underflow or overflow. 39 * 40 * If b is neither NaN nor zero and w is zero, or if b is infinite 41 * and w is finite and nonzero, _F_cplx_div_ix delivers an infinite 42 * result. If b is finite and w is infinite, _F_cplx_div_ix delivers 43 * a zero result. 44 * 45 * If b and w are both zero or both infinite, or if either b or w is 46 * NaN, _F_cplx_div_ix delivers NaN + I * NaN. C99 doesn't specify 47 * these cases. 48 * 49 * This implementation can raise spurious invalid operation, inexact, 50 * and division-by-zero exceptions. C99 allows this. 51 * 52 * Warning: Do not attempt to "optimize" this code by removing multi- 53 * plications by zero. 54 */ 55 56 #if !defined(sparc) && !defined(__sparc) 57 #error This code is for SPARC only 58 #endif 59 60 /* 61 * Return +1 if x is +Inf, -1 if x is -Inf, and 0 otherwise 62 */ 63 static int 64 testinff(float x) 65 { 66 union { 67 int i; 68 float f; 69 } xx; 70 71 xx.f = x; 72 return ((((xx.i << 1) - 0xff000000) == 0)? (1 | (xx.i >> 31)) : 0); 73 } 74 75 float _Complex 76 _F_cplx_div_ix(float b, float _Complex w) 77 { 78 float _Complex v; 79 union { 80 int i; 81 float f; 82 } cc, dd; 83 float c, d; 84 double r, x, y; 85 int i, j; 86 87 /* 88 * The following is equivalent to 89 * 90 * c = crealf(w); d = cimagf(w); 91 */ 92 c = ((float *)&w)[0]; 93 d = ((float *)&w)[1]; 94 95 r = (double)c * c + (double)d * d; 96 97 if (r == 0.0) { 98 /* w is zero; multiply b by 1/Re(w) - I * Im(w) */ 99 c = 1.0f / c; 100 j = testinff(b); 101 if (j) { /* b is infinite */ 102 b = j; 103 } 104 ((float *)&v)[0] = (b == 0.0f)? b * c : b * d; 105 ((float *)&v)[1] = b * c; 106 return (v); 107 } 108 109 r = (double)b / r; 110 x = (double)d * r; 111 y = (double)c * r; 112 113 if (x != x || y != y) { 114 /* 115 * x or y is NaN, so b and w can't both be finite and 116 * nonzero. Since we handled the case w = 0 above, the 117 * only case to check here is when w is infinite. 118 */ 119 i = testinff(c); 120 j = testinff(d); 121 if (i | j) { /* w is infinite */ 122 cc.f = c; 123 dd.f = d; 124 c = (cc.i < 0)? -0.0f : 0.0f; 125 d = (dd.i < 0)? -0.0f : 0.0f; 126 x = (double)d * b; 127 y = (double)c * b; 128 } 129 } 130 131 /* 132 * The following is equivalent to 133 * 134 * return x + I * y; 135 */ 136 ((float *)&v)[0] = (float)x; 137 ((float *)&v)[1] = (float)y; 138 return (v); 139 } 140