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