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 * _D_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, _D_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 scales to avoid premature underflow or overflow. 38 * 39 * If b is neither NaN nor zero and w is zero, or if b is infinite 40 * and w is finite and nonzero, _D_cplx_div_ix delivers an infinite 41 * result. If b is finite and w is infinite, _D_cplx_div_ix delivers 42 * a zero result. 43 * 44 * If b and w are both zero or both infinite, or if either b or w is 45 * NaN, _D_cplx_div_ix delivers NaN + I * NaN. C99 doesn't specify 46 * these cases. 47 * 48 * This implementation can raise spurious underflow, overflow, in- 49 * valid operation, inexact, and division-by-zero exceptions. C99 50 * 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 * scl[i].d = 2^(250*(4-i)) for i = 0, ..., 9 62 */ 63 static const union { 64 int i[2]; 65 double d; 66 } scl[9] = { 67 { 0x7e700000, 0 }, 68 { 0x6ed00000, 0 }, 69 { 0x5f300000, 0 }, 70 { 0x4f900000, 0 }, 71 { 0x3ff00000, 0 }, 72 { 0x30500000, 0 }, 73 { 0x20b00000, 0 }, 74 { 0x11100000, 0 }, 75 { 0x01700000, 0 } 76 }; 77 78 /* 79 * Return +1 if x is +Inf, -1 if x is -Inf, and 0 otherwise 80 */ 81 static int 82 testinf(double x) 83 { 84 union { 85 int i[2]; 86 double d; 87 } xx; 88 89 xx.d = x; 90 return (((((xx.i[0] << 1) - 0xffe00000) | xx.i[1]) == 0)? 91 (1 | (xx.i[0] >> 31)) : 0); 92 } 93 94 double _Complex 95 _D_cplx_div_ix(double b, double _Complex w) 96 { 97 double _Complex v; 98 union { 99 int i[2]; 100 double d; 101 } bb, cc, dd; 102 double c, d, sc, sd, r; 103 int hb, hc, hd, hw, i, j; 104 105 /* 106 * The following is equivalent to 107 * 108 * c = creal(w); d = cimag(w); 109 */ 110 c = ((double *)&w)[0]; 111 d = ((double *)&w)[1]; 112 113 /* extract high-order words to estimate |b| and |w| */ 114 bb.d = b; 115 hb = bb.i[0] & ~0x80000000; 116 117 cc.d = c; 118 dd.d = d; 119 hc = cc.i[0] & ~0x80000000; 120 hd = dd.i[0] & ~0x80000000; 121 hw = (hc > hd)? hc : hd; 122 123 /* check for special cases */ 124 if (hw >= 0x7ff00000) { /* w is inf or nan */ 125 i = testinf(c); 126 j = testinf(d); 127 if (i | j) { /* w is infinite */ 128 c = (cc.i[0] < 0)? -0.0 : 0.0; 129 d = (dd.i[0] < 0)? -0.0 : 0.0; 130 } else /* w is nan */ 131 b *= c * d; 132 ((double *)&v)[0] = b * d; 133 ((double *)&v)[1] = b * c; 134 return (v); 135 } 136 137 if (hw < 0x00100000) { 138 /* 139 * This nonsense is needed to work around some SPARC 140 * implementations of nonstandard mode; if both parts 141 * of w are subnormal, multiply them by one to force 142 * them to be flushed to zero when nonstandard mode 143 * is enabled. Sheesh. 144 */ 145 cc.d = c = c * 1.0; 146 dd.d = d = d * 1.0; 147 hc = cc.i[0] & ~0x80000000; 148 hd = dd.i[0] & ~0x80000000; 149 hw = (hc > hd)? hc : hd; 150 } 151 152 if (hw == 0 && (cc.i[1] | dd.i[1]) == 0) { 153 /* w is zero; multiply b by 1/Re(w) - I * Im(w) */ 154 c = 1.0 / c; 155 j = testinf(b); 156 if (j) { /* b is infinite */ 157 b = j; 158 } 159 ((double *)&v)[0] = (b == 0.0)? b * c : b * d; 160 ((double *)&v)[1] = b * c; 161 return (v); 162 } 163 164 if (hb >= 0x7ff00000) { /* a is inf or nan */ 165 ((double *)&v)[0] = b * d; 166 ((double *)&v)[1] = b * c; 167 return (v); 168 } 169 170 /* 171 * Compute the real and imaginary parts of the quotient, 172 * scaling to avoid overflow or underflow. 173 */ 174 hw = (hw - 0x38000000) >> 28; 175 sc = c * scl[hw + 4].d; 176 sd = d * scl[hw + 4].d; 177 r = sc * sc + sd * sd; 178 179 hb = (hb - 0x38000000) >> 28; 180 b = (b * scl[hb + 4].d) / r; 181 hb -= (hw + hw); 182 183 hc = (hc - 0x38000000) >> 28; 184 c = (c * scl[hc + 4].d) * b; 185 hc += hb; 186 187 hd = (hd - 0x38000000) >> 28; 188 d = (d * scl[hd + 4].d) * b; 189 hd += hb; 190 191 /* compensate for scaling */ 192 sc = scl[3].d; /* 2^250 */ 193 if (hc < 0) { 194 hc = -hc; 195 sc = scl[5].d; /* 2^-250 */ 196 } 197 while (hc--) 198 c *= sc; 199 200 sd = scl[3].d; 201 if (hd < 0) { 202 hd = -hd; 203 sd = scl[5].d; 204 } 205 while (hd--) 206 d *= sd; 207 208 ((double *)&v)[0] = d; 209 ((double *)&v)[1] = c; 210 return (v); 211 } 212