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 2004 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 /*
28*7c478bd9Sstevel@tonic-gate * _X_cplx_div_ix(b, w) returns (I * b) / w with infinities handled
29*7c478bd9Sstevel@tonic-gate * according to C99.
30*7c478bd9Sstevel@tonic-gate *
31*7c478bd9Sstevel@tonic-gate * If b and w are both finite and w is nonzero, _X_cplx_div_ix de-
32*7c478bd9Sstevel@tonic-gate * livers the complex quotient q according to the usual formula: let
33*7c478bd9Sstevel@tonic-gate * c = Re(w), and d = Im(w); then q = x + I * y where x = (b * d) / r
34*7c478bd9Sstevel@tonic-gate * and y = (b * c) / r with r = c * c + d * d. This implementation
35*7c478bd9Sstevel@tonic-gate * scales to avoid premature underflow or overflow.
36*7c478bd9Sstevel@tonic-gate *
37*7c478bd9Sstevel@tonic-gate * If b is neither NaN nor zero and w is zero, or if b is infinite
38*7c478bd9Sstevel@tonic-gate * and w is finite and nonzero, _X_cplx_div_ix delivers an infinite
39*7c478bd9Sstevel@tonic-gate * result. If b is finite and w is infinite, _X_cplx_div_ix delivers
40*7c478bd9Sstevel@tonic-gate * a zero result.
41*7c478bd9Sstevel@tonic-gate *
42*7c478bd9Sstevel@tonic-gate * If b and w are both zero or both infinite, or if either b or w is
43*7c478bd9Sstevel@tonic-gate * NaN, _X_cplx_div_ix delivers NaN + I * NaN. C99 doesn't specify
44*7c478bd9Sstevel@tonic-gate * these cases.
45*7c478bd9Sstevel@tonic-gate *
46*7c478bd9Sstevel@tonic-gate * This implementation can raise spurious underflow, overflow, in-
47*7c478bd9Sstevel@tonic-gate * valid operation, inexact, and division-by-zero exceptions. C99
48*7c478bd9Sstevel@tonic-gate * allows this.
49*7c478bd9Sstevel@tonic-gate */
50*7c478bd9Sstevel@tonic-gate
51*7c478bd9Sstevel@tonic-gate #if !defined(i386) && !defined(__i386) && !defined(__amd64)
52*7c478bd9Sstevel@tonic-gate #error This code is for x86 only
53*7c478bd9Sstevel@tonic-gate #endif
54*7c478bd9Sstevel@tonic-gate
55*7c478bd9Sstevel@tonic-gate /*
56*7c478bd9Sstevel@tonic-gate * scl[i].e = 2^(4080*(4-i)) for i = 0, ..., 9
57*7c478bd9Sstevel@tonic-gate */
58*7c478bd9Sstevel@tonic-gate static const union {
59*7c478bd9Sstevel@tonic-gate unsigned int i[3];
60*7c478bd9Sstevel@tonic-gate long double e;
61*7c478bd9Sstevel@tonic-gate } scl[9] = {
62*7c478bd9Sstevel@tonic-gate { 0, 0x80000000, 0x7fbf },
63*7c478bd9Sstevel@tonic-gate { 0, 0x80000000, 0x6fcf },
64*7c478bd9Sstevel@tonic-gate { 0, 0x80000000, 0x5fdf },
65*7c478bd9Sstevel@tonic-gate { 0, 0x80000000, 0x4fef },
66*7c478bd9Sstevel@tonic-gate { 0, 0x80000000, 0x3fff },
67*7c478bd9Sstevel@tonic-gate { 0, 0x80000000, 0x300f },
68*7c478bd9Sstevel@tonic-gate { 0, 0x80000000, 0x201f },
69*7c478bd9Sstevel@tonic-gate { 0, 0x80000000, 0x102f },
70*7c478bd9Sstevel@tonic-gate { 0, 0x80000000, 0x003f }
71*7c478bd9Sstevel@tonic-gate };
72*7c478bd9Sstevel@tonic-gate
73*7c478bd9Sstevel@tonic-gate /*
74*7c478bd9Sstevel@tonic-gate * Return +1 if x is +Inf, -1 if x is -Inf, and 0 otherwise
75*7c478bd9Sstevel@tonic-gate */
76*7c478bd9Sstevel@tonic-gate static int
testinfl(long double x)77*7c478bd9Sstevel@tonic-gate testinfl(long double x)
78*7c478bd9Sstevel@tonic-gate {
79*7c478bd9Sstevel@tonic-gate union {
80*7c478bd9Sstevel@tonic-gate int i[3];
81*7c478bd9Sstevel@tonic-gate long double e;
82*7c478bd9Sstevel@tonic-gate } xx;
83*7c478bd9Sstevel@tonic-gate
84*7c478bd9Sstevel@tonic-gate xx.e = x;
85*7c478bd9Sstevel@tonic-gate if ((xx.i[2] & 0x7fff) != 0x7fff || ((xx.i[1] << 1) | xx.i[0]) != 0)
86*7c478bd9Sstevel@tonic-gate return (0);
87*7c478bd9Sstevel@tonic-gate return (1 | ((xx.i[2] << 16) >> 31));
88*7c478bd9Sstevel@tonic-gate }
89*7c478bd9Sstevel@tonic-gate
90*7c478bd9Sstevel@tonic-gate long double _Complex
_X_cplx_div_ix(long double b,long double _Complex w)91*7c478bd9Sstevel@tonic-gate _X_cplx_div_ix(long double b, long double _Complex w)
92*7c478bd9Sstevel@tonic-gate {
93*7c478bd9Sstevel@tonic-gate long double _Complex v;
94*7c478bd9Sstevel@tonic-gate union {
95*7c478bd9Sstevel@tonic-gate int i[3];
96*7c478bd9Sstevel@tonic-gate long double e;
97*7c478bd9Sstevel@tonic-gate } bb, cc, dd;
98*7c478bd9Sstevel@tonic-gate long double c, d, sc, sd, r;
99*7c478bd9Sstevel@tonic-gate int eb, ec, ed, ew, i, j;
100*7c478bd9Sstevel@tonic-gate
101*7c478bd9Sstevel@tonic-gate /*
102*7c478bd9Sstevel@tonic-gate * The following is equivalent to
103*7c478bd9Sstevel@tonic-gate *
104*7c478bd9Sstevel@tonic-gate * c = creall(*w); d = cimagl(*w);
105*7c478bd9Sstevel@tonic-gate */
106*7c478bd9Sstevel@tonic-gate c = ((long double *)&w)[0];
107*7c478bd9Sstevel@tonic-gate d = ((long double *)&w)[1];
108*7c478bd9Sstevel@tonic-gate
109*7c478bd9Sstevel@tonic-gate /* extract exponents to estimate |z| and |w| */
110*7c478bd9Sstevel@tonic-gate bb.e = b;
111*7c478bd9Sstevel@tonic-gate eb = bb.i[2] & 0x7fff;
112*7c478bd9Sstevel@tonic-gate
113*7c478bd9Sstevel@tonic-gate cc.e = c;
114*7c478bd9Sstevel@tonic-gate dd.e = d;
115*7c478bd9Sstevel@tonic-gate ec = cc.i[2] & 0x7fff;
116*7c478bd9Sstevel@tonic-gate ed = dd.i[2] & 0x7fff;
117*7c478bd9Sstevel@tonic-gate ew = (ec > ed)? ec : ed;
118*7c478bd9Sstevel@tonic-gate
119*7c478bd9Sstevel@tonic-gate /* check for special cases */
120*7c478bd9Sstevel@tonic-gate if (ew >= 0x7fff) { /* w is inf or nan */
121*7c478bd9Sstevel@tonic-gate i = testinfl(c);
122*7c478bd9Sstevel@tonic-gate j = testinfl(d);
123*7c478bd9Sstevel@tonic-gate if (i | j) { /* w is infinite */
124*7c478bd9Sstevel@tonic-gate c = ((cc.i[2] << 16) < 0)? -0.0f : 0.0f;
125*7c478bd9Sstevel@tonic-gate d = ((dd.i[2] << 16) < 0)? -0.0f : 0.0f;
126*7c478bd9Sstevel@tonic-gate } else /* w is nan */
127*7c478bd9Sstevel@tonic-gate b += c + d;
128*7c478bd9Sstevel@tonic-gate ((long double *)&v)[0] = b * d;
129*7c478bd9Sstevel@tonic-gate ((long double *)&v)[1] = b * c;
130*7c478bd9Sstevel@tonic-gate return (v);
131*7c478bd9Sstevel@tonic-gate }
132*7c478bd9Sstevel@tonic-gate
133*7c478bd9Sstevel@tonic-gate if (ew == 0 && (cc.i[1] | cc.i[0] | dd.i[1] | dd.i[0]) == 0) {
134*7c478bd9Sstevel@tonic-gate /* w is zero; multiply b by 1/Re(w) - I * Im(w) */
135*7c478bd9Sstevel@tonic-gate c = 1.0f / c;
136*7c478bd9Sstevel@tonic-gate j = testinfl(b);
137*7c478bd9Sstevel@tonic-gate if (j) { /* b is infinite */
138*7c478bd9Sstevel@tonic-gate b = j;
139*7c478bd9Sstevel@tonic-gate }
140*7c478bd9Sstevel@tonic-gate ((long double *)&v)[0] = (b == 0.0f)? b * c : b * d;
141*7c478bd9Sstevel@tonic-gate ((long double *)&v)[1] = b * c;
142*7c478bd9Sstevel@tonic-gate return (v);
143*7c478bd9Sstevel@tonic-gate }
144*7c478bd9Sstevel@tonic-gate
145*7c478bd9Sstevel@tonic-gate if (eb >= 0x7fff) { /* a is inf or nan */
146*7c478bd9Sstevel@tonic-gate ((long double *)&v)[0] = b * d;
147*7c478bd9Sstevel@tonic-gate ((long double *)&v)[1] = b * c;
148*7c478bd9Sstevel@tonic-gate return (v);
149*7c478bd9Sstevel@tonic-gate }
150*7c478bd9Sstevel@tonic-gate
151*7c478bd9Sstevel@tonic-gate /*
152*7c478bd9Sstevel@tonic-gate * Compute the real and imaginary parts of the quotient,
153*7c478bd9Sstevel@tonic-gate * scaling to avoid overflow or underflow.
154*7c478bd9Sstevel@tonic-gate */
155*7c478bd9Sstevel@tonic-gate ew = (ew - 0x3800) >> 12;
156*7c478bd9Sstevel@tonic-gate sc = c * scl[ew + 4].e;
157*7c478bd9Sstevel@tonic-gate sd = d * scl[ew + 4].e;
158*7c478bd9Sstevel@tonic-gate r = sc * sc + sd * sd;
159*7c478bd9Sstevel@tonic-gate
160*7c478bd9Sstevel@tonic-gate eb = (eb - 0x3800) >> 12;
161*7c478bd9Sstevel@tonic-gate b = (b * scl[eb + 4].e) / r;
162*7c478bd9Sstevel@tonic-gate eb -= (ew + ew);
163*7c478bd9Sstevel@tonic-gate
164*7c478bd9Sstevel@tonic-gate ec = (ec - 0x3800) >> 12;
165*7c478bd9Sstevel@tonic-gate c = (c * scl[ec + 4].e) * b;
166*7c478bd9Sstevel@tonic-gate ec += eb;
167*7c478bd9Sstevel@tonic-gate
168*7c478bd9Sstevel@tonic-gate ed = (ed - 0x3800) >> 12;
169*7c478bd9Sstevel@tonic-gate d = (d * scl[ed + 4].e) * b;
170*7c478bd9Sstevel@tonic-gate ed += eb;
171*7c478bd9Sstevel@tonic-gate
172*7c478bd9Sstevel@tonic-gate /* compensate for scaling */
173*7c478bd9Sstevel@tonic-gate sc = scl[3].e; /* 2^4080 */
174*7c478bd9Sstevel@tonic-gate if (ec < 0) {
175*7c478bd9Sstevel@tonic-gate ec = -ec;
176*7c478bd9Sstevel@tonic-gate sc = scl[5].e; /* 2^-4080 */
177*7c478bd9Sstevel@tonic-gate }
178*7c478bd9Sstevel@tonic-gate while (ec--)
179*7c478bd9Sstevel@tonic-gate c *= sc;
180*7c478bd9Sstevel@tonic-gate
181*7c478bd9Sstevel@tonic-gate sd = scl[3].e;
182*7c478bd9Sstevel@tonic-gate if (ed < 0) {
183*7c478bd9Sstevel@tonic-gate ed = -ed;
184*7c478bd9Sstevel@tonic-gate sd = scl[5].e;
185*7c478bd9Sstevel@tonic-gate }
186*7c478bd9Sstevel@tonic-gate while (ed--)
187*7c478bd9Sstevel@tonic-gate d *= sd;
188*7c478bd9Sstevel@tonic-gate
189*7c478bd9Sstevel@tonic-gate ((long double *)&v)[0] = d;
190*7c478bd9Sstevel@tonic-gate ((long double *)&v)[1] = c;
191*7c478bd9Sstevel@tonic-gate return (v);
192*7c478bd9Sstevel@tonic-gate }
193