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_ix(v, b, w) sets *v = (I * *b / *w) with 31 * infinities handling according to C99. 32 * 33 * On SPARC V9, _Q_cplx_div_ix(b, w) returns (I * *b) / *w with in- 34 * finities handled according to C99. 35 * 36 * If b and w are both finite and w is nonzero, _Q_cplx_div_ix de- 37 * livers the complex quotient q according to the usual formula: let 38 * c = Re(w), and d = Im(w); then q = x + I * y where x = (b * d) / r 39 * and y = (-b * d) / r with r = c * c + d * d. This implementation 40 * scales to avoid premature underflow or overflow. 41 * 42 * If b is neither NaN nor zero and w is zero, or if b is infinite 43 * and w is finite and nonzero, _Q_cplx_div_ix delivers an infinite 44 * result. If b is finite and w is infinite, _Q_cplx_div_ix delivers 45 * a zero result. 46 * 47 * If b and w are both zero or both infinite, or if either b or w is 48 * NaN, _Q_cplx_div_ix delivers NaN + I * NaN. C99 doesn't specify 49 * these cases. 50 * 51 * This implementation can raise spurious underflow, overflow, in- 52 * valid operation, inexact, and division-by-zero exceptions. C99 53 * allows this. 54 */ 55 56 #if !defined(sparc) && !defined(__sparc) 57 #error This code is for SPARC only 58 #endif 59 60 extern void _Q_scl(long double *, int); 61 extern void _Q_scle(long double *, int); 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[4]; 71 long double q; 72 } xx; 73 74 xx.q = x; 75 return (((((xx.i[0] << 1) - 0xfffe0000) | xx.i[1] | xx.i[2] | xx.i[3]) 76 == 0)? (1 | (xx.i[0] >> 31)) : 0); 77 } 78 79 #ifdef __sparcv9 80 long double _Complex 81 _Q_cplx_div_ix(const long double *pb, const long double _Complex *w) 82 { 83 long double _Complex v; 84 #else 85 void 86 _Q_cplx_div_ix(long double _Complex *v, const long double *pb, 87 const long double _Complex *w) 88 { 89 #endif 90 union { 91 int i[4]; 92 long double q; 93 } bb, cc, dd; 94 long double b, c, d, sc, sd, r; 95 int hb, hc, hd, hw, i, j; 96 97 b = *pb; 98 99 /* 100 * The following is equivalent to 101 * 102 * c = creall(*w); d = cimagl(*w); 103 */ 104 c = ((long double *)w)[0]; 105 d = ((long double *)w)[1]; 106 107 /* extract high-order words to estimate |b| and |w| */ 108 bb.q = b; 109 hb = bb.i[0] & ~0x80000000; 110 111 cc.q = c; 112 dd.q = d; 113 hc = cc.i[0] & ~0x80000000; 114 hd = dd.i[0] & ~0x80000000; 115 hw = (hc > hd)? hc : hd; 116 117 /* check for special cases */ 118 if (hw >= 0x7fff0000) { /* w is inf or nan */ 119 i = testinfl(c); 120 j = testinfl(d); 121 if (i | j) { /* w is infinite */ 122 c = (cc.i[0] < 0)? -0.0l : 0.0l; 123 d = (dd.i[0] < 0)? -0.0l : 0.0l; 124 } else /* w is nan */ 125 b += c + d; 126 c *= b; 127 d *= b; 128 goto done; 129 } 130 131 if (hw == 0 && (cc.i[1] | cc.i[2] | cc.i[3] | 132 dd.i[1] | dd.i[2] | dd.i[3]) == 0) { 133 /* w is zero; multiply b by 1/Re(w) - I * Im(w) */ 134 c = 1.0l / c; 135 j = testinfl(b); 136 if (j) { /* b is infinite */ 137 b = j; 138 } 139 c *= b; 140 d = (b == 0.0l)? c : b * d; 141 goto done; 142 } 143 144 if (hb >= 0x7fff0000) { /* a is inf or nan */ 145 c *= b; 146 d *= b; 147 goto done; 148 } 149 150 /* 151 * Compute the real and imaginary parts of the quotient, 152 * scaling to avoid overflow or underflow. 153 */ 154 hw = (hw - 0x3fff0000) >> 16; 155 sc = c; 156 sd = d; 157 _Q_scl(&sc, -hw); 158 _Q_scl(&sd, -hw); 159 r = sc * sc + sd * sd; 160 161 hb = (hb - 0x3fff0000) >> 16; 162 _Q_scl(&b, -hb); 163 b /= r; 164 hb -= (hw + hw); 165 166 hc = (hc - 0x3fff0000) >> 16; 167 _Q_scl(&c, -hc); 168 c *= b; 169 hc += hb; 170 171 hd = (hd - 0x3fff0000) >> 16; 172 _Q_scl(&d, -hd); 173 d *= b; 174 hd += hb; 175 176 /* compensate for scaling */ 177 _Q_scle(&c, hc); 178 _Q_scle(&d, hd); 179 180 done: 181 #ifdef __sparcv9 182 ((long double *)&v)[0] = d; 183 ((long double *)&v)[1] = c; 184 return (v); 185 #else 186 ((long double *)v)[0] = d; 187 ((long double *)v)[1] = c; 188 #endif 189 } 190