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