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