1*7c478bd9Sstevel@tonic-gate /* 2*7c478bd9Sstevel@tonic-gate * CDDL HEADER START 3*7c478bd9Sstevel@tonic-gate * 4*7c478bd9Sstevel@tonic-gate * The contents of this file are subject to the terms of the 5*7c478bd9Sstevel@tonic-gate * Common Development and Distribution License, Version 1.0 only 6*7c478bd9Sstevel@tonic-gate * (the "License"). You may not use this file except in compliance 7*7c478bd9Sstevel@tonic-gate * with the License. 8*7c478bd9Sstevel@tonic-gate * 9*7c478bd9Sstevel@tonic-gate * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE 10*7c478bd9Sstevel@tonic-gate * or http://www.opensolaris.org/os/licensing. 11*7c478bd9Sstevel@tonic-gate * See the License for the specific language governing permissions 12*7c478bd9Sstevel@tonic-gate * and limitations under the License. 13*7c478bd9Sstevel@tonic-gate * 14*7c478bd9Sstevel@tonic-gate * When distributing Covered Code, include this CDDL HEADER in each 15*7c478bd9Sstevel@tonic-gate * file and include the License file at usr/src/OPENSOLARIS.LICENSE. 16*7c478bd9Sstevel@tonic-gate * If applicable, add the following below this CDDL HEADER, with the 17*7c478bd9Sstevel@tonic-gate * fields enclosed by brackets "[]" replaced with your own identifying 18*7c478bd9Sstevel@tonic-gate * information: Portions Copyright [yyyy] [name of copyright owner] 19*7c478bd9Sstevel@tonic-gate * 20*7c478bd9Sstevel@tonic-gate * CDDL HEADER END 21*7c478bd9Sstevel@tonic-gate */ 22*7c478bd9Sstevel@tonic-gate /* 23*7c478bd9Sstevel@tonic-gate * Copyright 2003 Sun Microsystems, Inc. All rights reserved. 24*7c478bd9Sstevel@tonic-gate * Use is subject to license terms. 25*7c478bd9Sstevel@tonic-gate */ 26*7c478bd9Sstevel@tonic-gate 27*7c478bd9Sstevel@tonic-gate #pragma ident "%Z%%M% %I% %E% SMI" 28*7c478bd9Sstevel@tonic-gate 29*7c478bd9Sstevel@tonic-gate /* 30*7c478bd9Sstevel@tonic-gate * On SPARC V8, _Q_cplx_div_rx(v, a, w) sets *v = *a / *w with in- 31*7c478bd9Sstevel@tonic-gate * finities handling according to C99. 32*7c478bd9Sstevel@tonic-gate * 33*7c478bd9Sstevel@tonic-gate * On SPARC V9, _Q_cplx_div_rx(a, w) returns *a / *w with infinities 34*7c478bd9Sstevel@tonic-gate * handled according to C99. 35*7c478bd9Sstevel@tonic-gate * 36*7c478bd9Sstevel@tonic-gate * If a and w are both finite and w is nonzero, _Q_cplx_div_rx de- 37*7c478bd9Sstevel@tonic-gate * livers the complex quotient q according to the usual formula: let 38*7c478bd9Sstevel@tonic-gate * c = Re(w), and d = Im(w); then q = x + I * y where x = (a * c) / r 39*7c478bd9Sstevel@tonic-gate * and y = (-a * d) / r with r = c * c + d * d. This implementation 40*7c478bd9Sstevel@tonic-gate * scales to avoid premature underflow or overflow. 41*7c478bd9Sstevel@tonic-gate * 42*7c478bd9Sstevel@tonic-gate * If a is neither NaN nor zero and w is zero, or if a is infinite 43*7c478bd9Sstevel@tonic-gate * and w is finite and nonzero, _Q_cplx_div_rx delivers an infinite 44*7c478bd9Sstevel@tonic-gate * result. If a is finite and w is infinite, _Q_cplx_div_rx delivers 45*7c478bd9Sstevel@tonic-gate * a zero result. 46*7c478bd9Sstevel@tonic-gate * 47*7c478bd9Sstevel@tonic-gate * If a and w are both zero or both infinite, or if either a or w is 48*7c478bd9Sstevel@tonic-gate * NaN, _Q_cplx_div_rx delivers NaN + I * NaN. C99 doesn't specify 49*7c478bd9Sstevel@tonic-gate * these cases. 50*7c478bd9Sstevel@tonic-gate * 51*7c478bd9Sstevel@tonic-gate * This implementation can raise spurious underflow, overflow, in- 52*7c478bd9Sstevel@tonic-gate * valid operation, inexact, and division-by-zero exceptions. C99 53*7c478bd9Sstevel@tonic-gate * allows this. 54*7c478bd9Sstevel@tonic-gate */ 55*7c478bd9Sstevel@tonic-gate 56*7c478bd9Sstevel@tonic-gate #if !defined(sparc) && !defined(__sparc) 57*7c478bd9Sstevel@tonic-gate #error This code is for SPARC only 58*7c478bd9Sstevel@tonic-gate #endif 59*7c478bd9Sstevel@tonic-gate 60*7c478bd9Sstevel@tonic-gate extern void _Q_scl(long double *, int); 61*7c478bd9Sstevel@tonic-gate extern void _Q_scle(long double *, int); 62*7c478bd9Sstevel@tonic-gate 63*7c478bd9Sstevel@tonic-gate /* 64*7c478bd9Sstevel@tonic-gate * Return +1 if x is +Inf, -1 if x is -Inf, and 0 otherwise 65*7c478bd9Sstevel@tonic-gate */ 66*7c478bd9Sstevel@tonic-gate static int 67*7c478bd9Sstevel@tonic-gate testinfl(long double x) 68*7c478bd9Sstevel@tonic-gate { 69*7c478bd9Sstevel@tonic-gate union { 70*7c478bd9Sstevel@tonic-gate int i[4]; 71*7c478bd9Sstevel@tonic-gate long double q; 72*7c478bd9Sstevel@tonic-gate } xx; 73*7c478bd9Sstevel@tonic-gate 74*7c478bd9Sstevel@tonic-gate xx.q = x; 75*7c478bd9Sstevel@tonic-gate return (((((xx.i[0] << 1) - 0xfffe0000) | xx.i[1] | xx.i[2] | xx.i[3]) 76*7c478bd9Sstevel@tonic-gate == 0)? (1 | (xx.i[0] >> 31)) : 0); 77*7c478bd9Sstevel@tonic-gate } 78*7c478bd9Sstevel@tonic-gate 79*7c478bd9Sstevel@tonic-gate #ifdef __sparcv9 80*7c478bd9Sstevel@tonic-gate long double _Complex 81*7c478bd9Sstevel@tonic-gate _Q_cplx_div_rx(const long double *pa, const long double _Complex *w) 82*7c478bd9Sstevel@tonic-gate { 83*7c478bd9Sstevel@tonic-gate long double _Complex v; 84*7c478bd9Sstevel@tonic-gate #else 85*7c478bd9Sstevel@tonic-gate void 86*7c478bd9Sstevel@tonic-gate _Q_cplx_div_rx(long double _Complex *v, const long double *pa, 87*7c478bd9Sstevel@tonic-gate const long double _Complex *w) 88*7c478bd9Sstevel@tonic-gate { 89*7c478bd9Sstevel@tonic-gate #endif 90*7c478bd9Sstevel@tonic-gate union { 91*7c478bd9Sstevel@tonic-gate int i[4]; 92*7c478bd9Sstevel@tonic-gate long double q; 93*7c478bd9Sstevel@tonic-gate } aa, cc, dd; 94*7c478bd9Sstevel@tonic-gate long double a, c, d, sc, sd, r; 95*7c478bd9Sstevel@tonic-gate int ha, hc, hd, hw, i, j; 96*7c478bd9Sstevel@tonic-gate 97*7c478bd9Sstevel@tonic-gate a = *pa; 98*7c478bd9Sstevel@tonic-gate 99*7c478bd9Sstevel@tonic-gate /* 100*7c478bd9Sstevel@tonic-gate * The following is equivalent to 101*7c478bd9Sstevel@tonic-gate * 102*7c478bd9Sstevel@tonic-gate * c = creall(*w); d = cimagl(*w); 103*7c478bd9Sstevel@tonic-gate */ 104*7c478bd9Sstevel@tonic-gate c = ((long double *)w)[0]; 105*7c478bd9Sstevel@tonic-gate d = ((long double *)w)[1]; 106*7c478bd9Sstevel@tonic-gate 107*7c478bd9Sstevel@tonic-gate /* extract high-order words to estimate |a| and |w| */ 108*7c478bd9Sstevel@tonic-gate aa.q = a; 109*7c478bd9Sstevel@tonic-gate ha = aa.i[0] & ~0x80000000; 110*7c478bd9Sstevel@tonic-gate 111*7c478bd9Sstevel@tonic-gate cc.q = c; 112*7c478bd9Sstevel@tonic-gate dd.q = d; 113*7c478bd9Sstevel@tonic-gate hc = cc.i[0] & ~0x80000000; 114*7c478bd9Sstevel@tonic-gate hd = dd.i[0] & ~0x80000000; 115*7c478bd9Sstevel@tonic-gate hw = (hc > hd)? hc : hd; 116*7c478bd9Sstevel@tonic-gate 117*7c478bd9Sstevel@tonic-gate /* check for special cases */ 118*7c478bd9Sstevel@tonic-gate if (hw >= 0x7fff0000) { /* w is inf or nan */ 119*7c478bd9Sstevel@tonic-gate i = testinfl(c); 120*7c478bd9Sstevel@tonic-gate j = testinfl(d); 121*7c478bd9Sstevel@tonic-gate if (i | j) { /* w is infinite */ 122*7c478bd9Sstevel@tonic-gate c = (cc.i[0] < 0)? -0.0l : 0.0l; 123*7c478bd9Sstevel@tonic-gate d = (dd.i[0] < 0)? -0.0l : 0.0l; 124*7c478bd9Sstevel@tonic-gate } else /* w is nan */ 125*7c478bd9Sstevel@tonic-gate a += c + d; 126*7c478bd9Sstevel@tonic-gate c *= a; 127*7c478bd9Sstevel@tonic-gate d *= -a; 128*7c478bd9Sstevel@tonic-gate goto done; 129*7c478bd9Sstevel@tonic-gate } 130*7c478bd9Sstevel@tonic-gate 131*7c478bd9Sstevel@tonic-gate if (hw == 0 && (cc.i[1] | cc.i[2] | cc.i[3] | 132*7c478bd9Sstevel@tonic-gate dd.i[1] | dd.i[2] | dd.i[3]) == 0) { 133*7c478bd9Sstevel@tonic-gate /* w is zero; multiply a by 1/Re(w) - I * Im(w) */ 134*7c478bd9Sstevel@tonic-gate c = 1.0l / c; 135*7c478bd9Sstevel@tonic-gate i = testinfl(a); 136*7c478bd9Sstevel@tonic-gate if (i) { /* a is infinite */ 137*7c478bd9Sstevel@tonic-gate a = i; 138*7c478bd9Sstevel@tonic-gate } 139*7c478bd9Sstevel@tonic-gate c *= a; 140*7c478bd9Sstevel@tonic-gate d = (a == 0.0l)? c : -a * d; 141*7c478bd9Sstevel@tonic-gate goto done; 142*7c478bd9Sstevel@tonic-gate } 143*7c478bd9Sstevel@tonic-gate 144*7c478bd9Sstevel@tonic-gate if (ha >= 0x7fff0000) { /* a is inf or nan */ 145*7c478bd9Sstevel@tonic-gate c *= a; 146*7c478bd9Sstevel@tonic-gate d *= -a; 147*7c478bd9Sstevel@tonic-gate goto done; 148*7c478bd9Sstevel@tonic-gate } 149*7c478bd9Sstevel@tonic-gate 150*7c478bd9Sstevel@tonic-gate /* 151*7c478bd9Sstevel@tonic-gate * Compute the real and imaginary parts of the quotient, 152*7c478bd9Sstevel@tonic-gate * scaling to avoid overflow or underflow. 153*7c478bd9Sstevel@tonic-gate */ 154*7c478bd9Sstevel@tonic-gate hw = (hw - 0x3fff0000) >> 16; 155*7c478bd9Sstevel@tonic-gate sc = c; 156*7c478bd9Sstevel@tonic-gate sd = d; 157*7c478bd9Sstevel@tonic-gate _Q_scl(&sc, -hw); 158*7c478bd9Sstevel@tonic-gate _Q_scl(&sd, -hw); 159*7c478bd9Sstevel@tonic-gate r = sc * sc + sd * sd; 160*7c478bd9Sstevel@tonic-gate 161*7c478bd9Sstevel@tonic-gate ha = (ha - 0x3fff0000) >> 16; 162*7c478bd9Sstevel@tonic-gate _Q_scl(&a, -ha); 163*7c478bd9Sstevel@tonic-gate a /= r; 164*7c478bd9Sstevel@tonic-gate ha -= (hw + hw); 165*7c478bd9Sstevel@tonic-gate 166*7c478bd9Sstevel@tonic-gate hc = (hc - 0x3fff0000) >> 16; 167*7c478bd9Sstevel@tonic-gate _Q_scl(&c, -hc); 168*7c478bd9Sstevel@tonic-gate c *= a; 169*7c478bd9Sstevel@tonic-gate hc += ha; 170*7c478bd9Sstevel@tonic-gate 171*7c478bd9Sstevel@tonic-gate hd = (hd - 0x3fff0000) >> 16; 172*7c478bd9Sstevel@tonic-gate _Q_scl(&d, -hd); 173*7c478bd9Sstevel@tonic-gate d *= -a; 174*7c478bd9Sstevel@tonic-gate hd += ha; 175*7c478bd9Sstevel@tonic-gate 176*7c478bd9Sstevel@tonic-gate /* compensate for scaling */ 177*7c478bd9Sstevel@tonic-gate _Q_scle(&c, hc); 178*7c478bd9Sstevel@tonic-gate _Q_scle(&d, hd); 179*7c478bd9Sstevel@tonic-gate 180*7c478bd9Sstevel@tonic-gate done: 181*7c478bd9Sstevel@tonic-gate #ifdef __sparcv9 182*7c478bd9Sstevel@tonic-gate ((long double *)&v)[0] = c; 183*7c478bd9Sstevel@tonic-gate ((long double *)&v)[1] = d; 184*7c478bd9Sstevel@tonic-gate return (v); 185*7c478bd9Sstevel@tonic-gate #else 186*7c478bd9Sstevel@tonic-gate ((long double *)v)[0] = c; 187*7c478bd9Sstevel@tonic-gate ((long double *)v)[1] = d; 188*7c478bd9Sstevel@tonic-gate #endif 189*7c478bd9Sstevel@tonic-gate } 190