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 2004 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 * _X_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, _X_cplx_div_ix de- 34 * livers the complex quotient q according to the usual formula: let 35 * c = Re(w), and d = Im(w); then q = x + I * y where x = (b * d) / r 36 * and y = (b * c) / r with r = c * c + d * d. This implementation 37 * 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, _X_cplx_div_ix delivers an infinite 41 * result. If b is finite and w is infinite, _X_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, _X_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 53 #if !defined(i386) && !defined(__i386) && !defined(__amd64) 54 #error This code is for x86 only 55 #endif 56 57 /* 58 * scl[i].e = 2^(4080*(4-i)) for i = 0, ..., 9 59 */ 60 static const union { 61 unsigned int i[3]; 62 long double e; 63 } scl[9] = { 64 { 0, 0x80000000, 0x7fbf }, 65 { 0, 0x80000000, 0x6fcf }, 66 { 0, 0x80000000, 0x5fdf }, 67 { 0, 0x80000000, 0x4fef }, 68 { 0, 0x80000000, 0x3fff }, 69 { 0, 0x80000000, 0x300f }, 70 { 0, 0x80000000, 0x201f }, 71 { 0, 0x80000000, 0x102f }, 72 { 0, 0x80000000, 0x003f } 73 }; 74 75 /* 76 * Return +1 if x is +Inf, -1 if x is -Inf, and 0 otherwise 77 */ 78 static int 79 testinfl(long double x) 80 { 81 union { 82 int i[3]; 83 long double e; 84 } xx; 85 86 xx.e = x; 87 if ((xx.i[2] & 0x7fff) != 0x7fff || ((xx.i[1] << 1) | xx.i[0]) != 0) 88 return (0); 89 return (1 | ((xx.i[2] << 16) >> 31)); 90 } 91 92 long double _Complex 93 _X_cplx_div_ix(long double b, long double _Complex w) 94 { 95 long double _Complex v; 96 union { 97 int i[3]; 98 long double e; 99 } bb, cc, dd; 100 long double c, d, sc, sd, r; 101 int eb, ec, ed, ew, i, j; 102 103 /* 104 * The following is equivalent to 105 * 106 * c = creall(*w); d = cimagl(*w); 107 */ 108 c = ((long double *)&w)[0]; 109 d = ((long double *)&w)[1]; 110 111 /* extract exponents to estimate |z| and |w| */ 112 bb.e = b; 113 eb = bb.i[2] & 0x7fff; 114 115 cc.e = c; 116 dd.e = d; 117 ec = cc.i[2] & 0x7fff; 118 ed = dd.i[2] & 0x7fff; 119 ew = (ec > ed)? ec : ed; 120 121 /* check for special cases */ 122 if (ew >= 0x7fff) { /* w is inf or nan */ 123 i = testinfl(c); 124 j = testinfl(d); 125 if (i | j) { /* w is infinite */ 126 c = ((cc.i[2] << 16) < 0)? -0.0f : 0.0f; 127 d = ((dd.i[2] << 16) < 0)? -0.0f : 0.0f; 128 } else /* w is nan */ 129 b += c + d; 130 ((long double *)&v)[0] = b * d; 131 ((long double *)&v)[1] = b * c; 132 return (v); 133 } 134 135 if (ew == 0 && (cc.i[1] | cc.i[0] | dd.i[1] | dd.i[0]) == 0) { 136 /* w is zero; multiply b by 1/Re(w) - I * Im(w) */ 137 c = 1.0f / c; 138 j = testinfl(b); 139 if (j) { /* b is infinite */ 140 b = j; 141 } 142 ((long double *)&v)[0] = (b == 0.0f)? b * c : b * d; 143 ((long double *)&v)[1] = b * c; 144 return (v); 145 } 146 147 if (eb >= 0x7fff) { /* a is inf or nan */ 148 ((long double *)&v)[0] = b * d; 149 ((long double *)&v)[1] = b * c; 150 return (v); 151 } 152 153 /* 154 * Compute the real and imaginary parts of the quotient, 155 * scaling to avoid overflow or underflow. 156 */ 157 ew = (ew - 0x3800) >> 12; 158 sc = c * scl[ew + 4].e; 159 sd = d * scl[ew + 4].e; 160 r = sc * sc + sd * sd; 161 162 eb = (eb - 0x3800) >> 12; 163 b = (b * scl[eb + 4].e) / r; 164 eb -= (ew + ew); 165 166 ec = (ec - 0x3800) >> 12; 167 c = (c * scl[ec + 4].e) * b; 168 ec += eb; 169 170 ed = (ed - 0x3800) >> 12; 171 d = (d * scl[ed + 4].e) * b; 172 ed += eb; 173 174 /* compensate for scaling */ 175 sc = scl[3].e; /* 2^4080 */ 176 if (ec < 0) { 177 ec = -ec; 178 sc = scl[5].e; /* 2^-4080 */ 179 } 180 while (ec--) 181 c *= sc; 182 183 sd = scl[3].e; 184 if (ed < 0) { 185 ed = -ed; 186 sd = scl[5].e; 187 } 188 while (ed--) 189 d *= sd; 190 191 ((long double *)&v)[0] = d; 192 ((long double *)&v)[1] = c; 193 return (v); 194 } 195