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_rx(v, a, w) sets *v = *a / *w with in- 31 * finities handling according to C99. 32 * 33 * On SPARC V9, _Q_cplx_div_rx(a, w) returns *a / *w with infinities 34 * handled according to C99. 35 * 36 * If a and w are both finite and w is nonzero, _Q_cplx_div_rx 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 = (a * c) / r 39 * and y = (-a * d) / r with r = c * c + d * d. This implementation 40 * scales to avoid premature underflow or overflow. 41 * 42 * If a is neither NaN nor zero and w is zero, or if a is infinite 43 * and w is finite and nonzero, _Q_cplx_div_rx delivers an infinite 44 * result. If a is finite and w is infinite, _Q_cplx_div_rx delivers 45 * a zero result. 46 * 47 * If a and w are both zero or both infinite, or if either a or w is 48 * NaN, _Q_cplx_div_rx 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_rx(const long double *pa, const long double _Complex *w) 82 { 83 long double _Complex v; 84 #else 85 void 86 _Q_cplx_div_rx(long double _Complex *v, const long double *pa, 87 const long double _Complex *w) 88 { 89 #endif 90 union { 91 int i[4]; 92 long double q; 93 } aa, cc, dd; 94 long double a, c, d, sc, sd, r; 95 int ha, hc, hd, hw, i, j; 96 97 a = *pa; 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 |a| and |w| */ 108 aa.q = a; 109 ha = aa.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 a += c + d; 126 c *= a; 127 d *= -a; 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 a by 1/Re(w) - I * Im(w) */ 134 c = 1.0l / c; 135 i = testinfl(a); 136 if (i) { /* a is infinite */ 137 a = i; 138 } 139 c *= a; 140 d = (a == 0.0l)? c : -a * d; 141 goto done; 142 } 143 144 if (ha >= 0x7fff0000) { /* a is inf or nan */ 145 c *= a; 146 d *= -a; 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 ha = (ha - 0x3fff0000) >> 16; 162 _Q_scl(&a, -ha); 163 a /= r; 164 ha -= (hw + hw); 165 166 hc = (hc - 0x3fff0000) >> 16; 167 _Q_scl(&c, -hc); 168 c *= a; 169 hc += ha; 170 171 hd = (hd - 0x3fff0000) >> 16; 172 _Q_scl(&d, -hd); 173 d *= -a; 174 hd += ha; 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] = c; 183 ((long double *)&v)[1] = d; 184 return (v); 185 #else 186 ((long double *)v)[0] = c; 187 ((long double *)v)[1] = d; 188 #endif 189 } 190