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 /*
28 * _D_cplx_div(z, w) returns z / w with infinities handled according
29 * to C99.
30 *
31 * If z and w are both finite and w is nonzero, _D_cplx_div(z, w)
32 * delivers the complex quotient q according to the usual formula:
33 * let a = Re(z), b = Im(z), c = Re(w), and d = Im(w); then q = x +
34 * I * y where x = (a * c + b * d) / r and y = (b * c - a * d) / r
35 * with r = c * c + d * d. This implementation computes intermediate
36 * results in extended precision to avoid premature underflow or over-
37 * flow.
38 *
39 * If z is neither NaN nor zero and w is zero, or if z is infinite
40 * and w is finite and nonzero, _D_cplx_div delivers an infinite
41 * result. If z is finite and w is infinite, _D_cplx_div delivers
42 * a zero result.
43 *
44 * If z and w are both zero or both infinite, or if either z or w is
45 * a complex NaN, _D_cplx_div delivers NaN + I * NaN. C99 doesn't
46 * specify these cases.
47 *
48 * This implementation can raise spurious invalid operation, inexact,
49 * and division-by-zero exceptions. C99 allows this.
50 *
51 * Warning: Do not attempt to "optimize" this code by removing multi-
52 * plications by zero.
53 */
54
55 #if !defined(i386) && !defined(__i386) && !defined(__amd64)
56 #error This code is for x86 only
57 #endif
58
59 static union {
60 int i;
61 float f;
62 } inf = {
63 0x7f800000
64 };
65
66 /*
67 * Return +1 if x is +Inf, -1 if x is -Inf, and 0 otherwise
68 */
69 static int
testinf(double x)70 testinf(double x)
71 {
72 union {
73 int i[2];
74 double d;
75 } xx;
76
77 xx.d = x;
78 return (((((xx.i[1] << 1) - 0xffe00000) | xx.i[0]) == 0)?
79 (1 | (xx.i[1] >> 31)) : 0);
80 }
81
82 double _Complex
_D_cplx_div(double _Complex z,double _Complex w)83 _D_cplx_div(double _Complex z, double _Complex w)
84 {
85 double _Complex v;
86 union {
87 int i[2];
88 double d;
89 } cc, dd;
90 double a, b, c, d;
91 long double r, x, y;
92 int i, j, recalc;
93
94 /*
95 * The following is equivalent to
96 *
97 * a = creal(z); b = cimag(z);
98 * c = creal(w); d = cimag(w);
99 */
100 /* LINTED alignment */
101 a = ((double *)&z)[0];
102 /* LINTED alignment */
103 b = ((double *)&z)[1];
104 /* LINTED alignment */
105 c = ((double *)&w)[0];
106 /* LINTED alignment */
107 d = ((double *)&w)[1];
108
109 r = (long double)c * c + (long double)d * d;
110
111 if (r == 0.0f) {
112 /* w is zero; multiply z by 1/Re(w) - I * Im(w) */
113 c = 1.0f / c;
114 i = testinf(a);
115 j = testinf(b);
116 if (i | j) { /* z is infinite */
117 a = i;
118 b = j;
119 }
120 /* LINTED alignment */
121 ((double *)&v)[0] = a * c + b * d;
122 /* LINTED alignment */
123 ((double *)&v)[1] = b * c - a * d;
124 return (v);
125 }
126
127 r = 1.0f / r;
128 x = ((long double)a * c + (long double)b * d) * r;
129 y = ((long double)b * c - (long double)a * d) * r;
130
131 if (x != x && y != y) {
132 /*
133 * Both x and y are NaN, so z and w can't both be finite
134 * and nonzero. Since we handled the case w = 0 above,
135 * the only cases to check here are when one of z or w
136 * is infinite.
137 */
138 r = 1.0f;
139 recalc = 0;
140 i = testinf(a);
141 j = testinf(b);
142 if (i | j) { /* z is infinite */
143 /* "factor out" infinity */
144 a = i;
145 b = j;
146 r = inf.f;
147 recalc = 1;
148 }
149 i = testinf(c);
150 j = testinf(d);
151 if (i | j) { /* w is infinite */
152 /*
153 * "factor out" infinity, being careful to preserve
154 * signs of finite values
155 */
156 cc.d = c;
157 dd.d = d;
158 c = i? i : ((cc.i[1] < 0)? -0.0f : 0.0f);
159 d = j? j : ((dd.i[1] < 0)? -0.0f : 0.0f);
160 r *= 0.0f;
161 recalc = 1;
162 }
163 if (recalc) {
164 x = ((long double)a * c + (long double)b * d) * r;
165 y = ((long double)b * c - (long double)a * d) * r;
166 }
167 }
168
169 /*
170 * The following is equivalent to
171 *
172 * return x + I * y;
173 */
174 /* LINTED alignment */
175 ((double *)&v)[0] = (double)x;
176 /* LINTED alignment */
177 ((double *)&v)[1] = (double)y;
178 return (v);
179 }
180