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 #pragma ident "%Z%%M% %I% %E% SMI" 28 29 /* 30 * On SPARC V8, _Q_cplx_div(v, z, w) sets *v = *z / *w with infin- 31 * ities handling according to C99. 32 * 33 * On SPARC V9, _Q_cplx_div(z, w) returns *z / *w with infinities 34 * handled according to C99. 35 * 36 * If z and w are both finite and w is nonzero, _Q_cplx_div delivers 37 * the complex quotient q according to the usual formula: let a = 38 * Re(z), b = Im(z), c = Re(w), and d = Im(w); then q = x + I * y 39 * where x = (a * c + b * d) / r and y = (b * c - a * d) / r with 40 * r = c * c + d * d. This implementation scales to avoid premature 41 * underflow or overflow. 42 * 43 * If z is neither NaN nor zero and w is zero, or if z is infinite 44 * and w is finite and nonzero, _Q_cplx_div delivers an infinite 45 * result. If z is finite and w is infinite, _Q_cplx_div delivers 46 * a zero result. 47 * 48 * If z and w are both zero or both infinite, or if either z or w is 49 * a complex NaN, _Q_cplx_div delivers NaN + I * NaN. C99 doesn't 50 * specify these cases. 51 * 52 * This implementation can raise spurious underflow, overflow, in- 53 * valid operation, inexact, and division-by-zero exceptions. C99 54 * allows this. 55 */ 56 57 #if !defined(sparc) && !defined(__sparc) 58 #error This code is for SPARC only 59 #endif 60 61 static union { 62 int i[4]; 63 long double q; 64 } inf = { 65 0x7fff0000, 0, 0, 0 66 }; 67 68 /* 69 * Return +1 if x is +Inf, -1 if x is -Inf, and 0 otherwise 70 */ 71 static int 72 testinfl(long double x) 73 { 74 union { 75 int i[4]; 76 long double q; 77 } xx; 78 79 xx.q = x; 80 return (((((xx.i[0] << 1) - 0xfffe0000) | xx.i[1] | xx.i[2] | xx.i[3]) 81 == 0)? (1 | (xx.i[0] >> 31)) : 0); 82 } 83 84 #ifdef __sparcv9 85 long double _Complex 86 _Q_cplx_div(const long double _Complex *z, const long double _Complex *w) 87 { 88 long double _Complex v; 89 #else 90 void 91 _Q_cplx_div(long double _Complex *v, const long double _Complex *z, 92 const long double _Complex *w) 93 { 94 #endif 95 union { 96 int i[4]; 97 long double q; 98 } aa, bb, cc, dd, ss; 99 long double a, b, c, d, r; 100 int ha, hb, hc, hd, hz, hw, hs, i, j; 101 102 /* 103 * The following is equivalent to 104 * 105 * a = creall(*z); b = cimagl(*z); 106 * c = creall(*w); d = cimagl(*w); 107 */ 108 a = ((long double *)z)[0]; 109 b = ((long double *)z)[1]; 110 c = ((long double *)w)[0]; 111 d = ((long double *)w)[1]; 112 113 /* extract high-order words to estimate |z| and |w| */ 114 aa.q = a; 115 bb.q = b; 116 ha = aa.i[0] & ~0x80000000; 117 hb = bb.i[0] & ~0x80000000; 118 hz = (ha > hb)? ha : hb; 119 120 cc.q = c; 121 dd.q = d; 122 hc = cc.i[0] & ~0x80000000; 123 hd = dd.i[0] & ~0x80000000; 124 hw = (hc > hd)? hc : hd; 125 126 /* check for special cases */ 127 if (hw >= 0x7fff0000) { /* w is inf or nan */ 128 r = 0.0l; 129 i = testinfl(c); 130 j = testinfl(d); 131 if (i | j) { /* w is infinite */ 132 /* 133 * "factor out" infinity, being careful to preserve 134 * signs of finite values 135 */ 136 c = i? i : ((cc.i[0] < 0)? -0.0l : 0.0l); 137 d = j? j : ((dd.i[0] < 0)? -0.0l : 0.0l); 138 if (hz >= 0x7ffe0000) { 139 /* scale to avoid overflow below */ 140 c *= 0.5l; 141 d *= 0.5l; 142 } 143 } 144 goto done; 145 } 146 147 if (hw == 0 && (cc.i[1] | cc.i[2] | cc.i[3] | 148 dd.i[1] | dd.i[2] | dd.i[3]) == 0) { 149 /* w is zero; multiply z by 1/Re(w) - I * Im(w) */ 150 r = 1.0l; 151 c = 1.0l / c; 152 i = testinfl(a); 153 j = testinfl(b); 154 if (i | j) { /* z is infinite */ 155 a = i; 156 b = j; 157 } 158 goto done; 159 } 160 161 if (hz >= 0x7fff0000) { /* z is inf or nan */ 162 r = 1.0l; 163 i = testinfl(a); 164 j = testinfl(b); 165 if (i | j) { /* z is infinite */ 166 a = i; 167 b = j; 168 r = inf.q; 169 } 170 goto done; 171 } 172 173 /* 174 * Scale c and d to compute 1/|w|^2 and the real and imaginary 175 * parts of the quotient. 176 */ 177 hs = (((hw >> 2) - hw) + 0x6ffd7fff) & 0xffff0000; 178 if (hz < 0x00ea0000) { /* |z| < 2^-16149 */ 179 if (((hw - 0x3e380000) | (0x40e90000 - hw)) >= 0) 180 hs = (((0x40e90000 - hw) >> 1) & 0xffff0000) 181 + 0x3fff0000; 182 } 183 ss.i[0] = hs; 184 ss.i[1] = ss.i[2] = ss.i[3] = 0; 185 186 c *= ss.q; 187 d *= ss.q; 188 r = 1.0l / (c * c + d * d); 189 190 c *= ss.q; 191 d *= ss.q; 192 193 done: 194 #ifdef __sparcv9 195 ((long double *)&v)[0] = (a * c + b * d) * r; 196 ((long double *)&v)[1] = (b * c - a * d) * r; 197 return (v); 198 #else 199 ((long double *)v)[0] = (a * c + b * d) * r; 200 ((long double *)v)[1] = (b * c - a * d) * r; 201 #endif 202 } 203