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