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