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(z, w) returns z / w with infinities handled according 29 * to C99. 30 * 31 * If z and w are both finite and w is nonzero, _X_cplx_div delivers 32 * the complex quotient q according to the usual formula: let a = 33 * Re(z), b = Im(z), c = Re(w), and d = Im(w); then q = x + I * y 34 * where x = (a * c + b * d) / r and y = (b * c - a * d) / r with 35 * r = c * c + d * d. This implementation scales to avoid premature 36 * underflow or overflow. 37 * 38 * If z is neither NaN nor zero and w is zero, or if z is infinite 39 * and w is finite and nonzero, _X_cplx_div delivers an infinite 40 * result. If z is finite and w is infinite, _X_cplx_div delivers 41 * a zero result. 42 * 43 * If z and w are both zero or both infinite, or if either z or w is 44 * a complex NaN, _X_cplx_div delivers NaN + I * NaN. C99 doesn't 45 * specify these cases. 46 * 47 * This implementation can raise spurious underflow, overflow, in- 48 * valid operation, inexact, and division-by-zero exceptions. C99 49 * allows this. 50 */ 51 52 #if !defined(i386) && !defined(__i386) && !defined(__amd64) 53 #error This code is for x86 only 54 #endif 55 56 static union { 57 int i; 58 float f; 59 } inf = { 60 0x7f800000 61 }; 62 63 /* 64 * Return +1 if x is +Inf, -1 if x is -Inf, and 0 otherwise 65 */ 66 static int 67 testinfl(long double x) 68 { 69 union { 70 int i[3]; 71 long double e; 72 } xx; 73 74 xx.e = x; 75 if ((xx.i[2] & 0x7fff) != 0x7fff || ((xx.i[1] << 1) | xx.i[0]) != 0) 76 return (0); 77 return (1 | ((xx.i[2] << 16) >> 31)); 78 } 79 80 long double _Complex 81 _X_cplx_div(long double _Complex z, long double _Complex w) 82 { 83 long double _Complex v; 84 union { 85 int i[3]; 86 long double e; 87 } aa, bb, cc, dd, ss; 88 long double a, b, c, d, r; 89 int ea, eb, ec, ed, ez, ew, es, i, j; 90 91 /* 92 * The following is equivalent to 93 * 94 * a = creall(*z); b = cimagl(*z); 95 * c = creall(*w); d = cimagl(*w); 96 */ 97 a = ((long double *)&z)[0]; 98 b = ((long double *)&z)[1]; 99 c = ((long double *)&w)[0]; 100 d = ((long double *)&w)[1]; 101 102 /* extract exponents to estimate |z| and |w| */ 103 aa.e = a; 104 bb.e = b; 105 ea = aa.i[2] & 0x7fff; 106 eb = bb.i[2] & 0x7fff; 107 ez = (ea > eb)? ea : eb; 108 109 cc.e = c; 110 dd.e = d; 111 ec = cc.i[2] & 0x7fff; 112 ed = dd.i[2] & 0x7fff; 113 ew = (ec > ed)? ec : ed; 114 115 /* check for special cases */ 116 if (ew >= 0x7fff) { /* w is inf or nan */ 117 r = 0.0f; 118 i = testinfl(c); 119 j = testinfl(d); 120 if (i | j) { /* w is infinite */ 121 /* 122 * "factor out" infinity, being careful to preserve 123 * signs of finite values 124 */ 125 c = i? i : (((cc.i[2] << 16) < 0)? -0.0f : 0.0f); 126 d = j? j : (((dd.i[2] << 16) < 0)? -0.0f : 0.0f); 127 if (ez >= 0x7ffe) { 128 /* scale to avoid overflow below */ 129 c *= 0.5f; 130 d *= 0.5f; 131 } 132 } 133 ((long double *)&v)[0] = (a * c + b * d) * r; 134 ((long double *)&v)[1] = (b * c - a * d) * r; 135 return (v); 136 } 137 138 if (ew == 0 && (cc.i[1] | cc.i[0] | dd.i[1] | dd.i[0]) == 0) { 139 /* w is zero; multiply z by 1/Re(w) - I * Im(w) */ 140 c = 1.0f / c; 141 i = testinfl(a); 142 j = testinfl(b); 143 if (i | j) { /* z is infinite */ 144 a = i; 145 b = j; 146 } 147 ((long double *)&v)[0] = a * c + b * d; 148 ((long double *)&v)[1] = b * c - a * d; 149 return (v); 150 } 151 152 if (ez >= 0x7fff) { /* z is inf or nan */ 153 i = testinfl(a); 154 j = testinfl(b); 155 if (i | j) { /* z is infinite */ 156 a = i; 157 b = j; 158 r = inf.f; 159 } 160 ((long double *)&v)[0] = a * c + b * d; 161 ((long double *)&v)[1] = b * c - a * d; 162 return (v); 163 } 164 165 /* 166 * Scale c and d to compute 1/|w|^2 and the real and imaginary 167 * parts of the quotient. 168 */ 169 es = ((ew >> 2) - ew) + 0x6ffd; 170 if (ez < 0x0086) { /* |z| < 2^-16249 */ 171 if (((ew - 0x3efe) | (0x4083 - ew)) >= 0) 172 es = ((0x4083 - ew) >> 1) + 0x3fff; 173 } 174 ss.i[2] = es; 175 ss.i[1] = 0x80000000; 176 ss.i[0] = 0; 177 178 c *= ss.e; 179 d *= ss.e; 180 r = 1.0f / (c * c + d * d); 181 182 c *= ss.e; 183 d *= ss.e; 184 185 ((long double *)&v)[0] = (a * c + b * d) * r; 186 ((long double *)&v)[1] = (b * c - a * d) * r; 187 return (v); 188 } 189