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