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 * _D_cplx_div_rx(a, w) returns a / w with infinities handled according
31 * to C99.
32 *
33 * If a and w are both finite and w is nonzero, _D_cplx_div_rx(a, w)
34 * delivers the complex quotient q according to the usual formula:
35 * let c = Re(w), and d = Im(w); then q = x + I * y where x = (a * c)
36 * / r and y = (-a * d) / r with r = c * c + d * d. This implementa-
37 * tion scales to avoid premature underflow or overflow.
38 *
39 * If a is neither NaN nor zero and w is zero, or if a is infinite
40 * and w is finite and nonzero, _D_cplx_div_rx delivers an infinite
41 * result. If a is finite and w is infinite, _D_cplx_div_rx delivers
42 * a zero result.
43 *
44 * If a and w are both zero or both infinite, or if either a or w is
45 * NaN, _D_cplx_div_rx delivers NaN + I * NaN. C99 doesn't specify
46 * these cases.
47 *
48 * This implementation can raise spurious underflow, overflow, in-
49 * valid operation, inexact, and division-by-zero exceptions. C99
50 * allows this.
51 *
52 * Warning: Do not attempt to "optimize" this code by removing multi-
53 * plications by zero.
54 */
55
56 #if !defined(sparc) && !defined(__sparc)
57 #error This code is for SPARC only
58 #endif
59
60 /*
61 * scl[i].d = 2^(250*(4-i)) for i = 0, ..., 9
62 */
63 static const union {
64 int i[2];
65 double d;
66 } scl[9] = {
67 { 0x7e700000, 0 },
68 { 0x6ed00000, 0 },
69 { 0x5f300000, 0 },
70 { 0x4f900000, 0 },
71 { 0x3ff00000, 0 },
72 { 0x30500000, 0 },
73 { 0x20b00000, 0 },
74 { 0x11100000, 0 },
75 { 0x01700000, 0 }
76 };
77
78 /*
79 * Return +1 if x is +Inf, -1 if x is -Inf, and 0 otherwise
80 */
81 static int
testinf(double x)82 testinf(double x)
83 {
84 union {
85 int i[2];
86 double d;
87 } xx;
88
89 xx.d = x;
90 return (((((xx.i[0] << 1) - 0xffe00000) | xx.i[1]) == 0)?
91 (1 | (xx.i[0] >> 31)) : 0);
92 }
93
94 double _Complex
_D_cplx_div_rx(double a,double _Complex w)95 _D_cplx_div_rx(double a, double _Complex w)
96 {
97 double _Complex v;
98 union {
99 int i[2];
100 double d;
101 } aa, cc, dd;
102 double c, d, sc, sd, r;
103 int ha, hc, hd, hw, i, j;
104
105 /*
106 * The following is equivalent to
107 *
108 * c = creal(w); d = cimag(w);
109 */
110 c = ((double *)&w)[0];
111 d = ((double *)&w)[1];
112
113 /* extract high-order words to estimate |a| and |w| */
114 aa.d = a;
115 ha = aa.i[0] & ~0x80000000;
116
117 cc.d = c;
118 dd.d = d;
119 hc = cc.i[0] & ~0x80000000;
120 hd = dd.i[0] & ~0x80000000;
121 hw = (hc > hd)? hc : hd;
122
123 /* check for special cases */
124 if (hw >= 0x7ff00000) { /* w is inf or nan */
125 i = testinf(c);
126 j = testinf(d);
127 if (i | j) { /* w is infinite */
128 c = (cc.i[0] < 0)? -0.0 : 0.0;
129 d = (dd.i[0] < 0)? -0.0 : 0.0;
130 } else /* w is nan */
131 a *= c * d;
132 ((double *)&v)[0] = a * c;
133 ((double *)&v)[1] = -a * d;
134 return (v);
135 }
136
137 if (hw < 0x00100000) {
138 /*
139 * This nonsense is needed to work around some SPARC
140 * implementations of nonstandard mode; if both parts
141 * of w are subnormal, multiply them by one to force
142 * them to be flushed to zero when nonstandard mode
143 * is enabled. Sheesh.
144 */
145 cc.d = c = c * 1.0;
146 dd.d = d = d * 1.0;
147 hc = cc.i[0] & ~0x80000000;
148 hd = dd.i[0] & ~0x80000000;
149 hw = (hc > hd)? hc : hd;
150 }
151
152 if (hw == 0 && (cc.i[1] | dd.i[1]) == 0) {
153 /* w is zero; multiply a by 1/Re(w) - I * Im(w) */
154 c = 1.0 / c;
155 i = testinf(a);
156 if (i) { /* a is infinite */
157 a = i;
158 }
159 ((double *)&v)[0] = a * c;
160 ((double *)&v)[1] = (a == 0.0)? a * c : -a * d;
161 return (v);
162 }
163
164 if (ha >= 0x7ff00000) { /* a is inf or nan */
165 ((double *)&v)[0] = a * c;
166 ((double *)&v)[1] = -a * d;
167 return (v);
168 }
169
170 /*
171 * Compute the real and imaginary parts of the quotient,
172 * scaling to avoid overflow or underflow.
173 */
174 hw = (hw - 0x38000000) >> 28;
175 sc = c * scl[hw + 4].d;
176 sd = d * scl[hw + 4].d;
177 r = sc * sc + sd * sd;
178
179 ha = (ha - 0x38000000) >> 28;
180 a = (a * scl[ha + 4].d) / r;
181 ha -= (hw + hw);
182
183 hc = (hc - 0x38000000) >> 28;
184 c = (c * scl[hc + 4].d) * a;
185 hc += ha;
186
187 hd = (hd - 0x38000000) >> 28;
188 d = -(d * scl[hd + 4].d) * a;
189 hd += ha;
190
191 /* compensate for scaling */
192 sc = scl[3].d; /* 2^250 */
193 if (hc < 0) {
194 hc = -hc;
195 sc = scl[5].d; /* 2^-250 */
196 }
197 while (hc--)
198 c *= sc;
199
200 sd = scl[3].d;
201 if (hd < 0) {
202 hd = -hd;
203 sd = scl[5].d;
204 }
205 while (hd--)
206 d *= sd;
207
208 ((double *)&v)[0] = c;
209 ((double *)&v)[1] = d;
210 return (v);
211 }
212