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_mul(z, w) returns z * w with infinities handled according
31 * to C99.
32 *
33 * If z and w are both finite, _X_cplx_mul(z, w) delivers the complex
34 * product according to the usual formula: let a = Re(z), b = Im(z),
35 * c = Re(w), and d = Im(w); then _X_cplx_mul(z, w) delivers x + I * y
36 * where x = a * c - b * d and y = a * d + b * c. Note that if both
37 * ac and bd overflow, then at least one of ad or bc must also over-
38 * flow, and vice versa, so that if one component of the product is
39 * NaN, the other is infinite. (Such a value is considered infinite
40 * according to C99.)
41 *
42 * If one of z or w is infinite and the other is either finite nonzero
43 * or infinite, _X_cplx_mul delivers an infinite result. If one factor
44 * is infinite and the other is zero, _X_cplx_mul delivers NaN + I * NaN.
45 * C99 doesn't specify the latter case.
46 *
47 * C99 also doesn't specify what should happen if either z or w is a
48 * complex NaN (i.e., neither finite nor infinite). This implementation
49 * delivers NaN + I * NaN in this case.
50 *
51 * This implementation can raise spurious underflow, overflow, invalid
52 * operation, and inexact exceptions. C99 allows this.
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
testinfl(long double x)70 testinfl(long double x)
71 {
72 union {
73 int i[3];
74 long double e;
75 } xx;
76
77 xx.e = x;
78 if ((xx.i[2] & 0x7fff) != 0x7fff || ((xx.i[1] << 1) | xx.i[0]) != 0)
79 return (0);
80 return (1 | ((xx.i[2] << 16) >> 31));
81 }
82
83 long double _Complex
_X_cplx_mul(long double _Complex z,long double _Complex w)84 _X_cplx_mul(long double _Complex z, long double _Complex w)
85 {
86 long double _Complex v;
87 long double a, b, c, d, x, y;
88 int recalc, i, j;
89
90 /*
91 * The following is equivalent to
92 *
93 * a = creall(z); b = cimagl(z);
94 * c = creall(w); d = cimagl(w);
95 */
96 a = ((long double *)&z)[0];
97 b = ((long double *)&z)[1];
98 c = ((long double *)&w)[0];
99 d = ((long double *)&w)[1];
100
101 x = a * c - b * d;
102 y = a * d + b * c;
103
104 if (x != x && y != y) {
105 /*
106 * Both x and y are NaN, so z and w can't both be finite.
107 * If at least one of z or w is a complex NaN, and neither
108 * is infinite, then we might as well deliver NaN + I * NaN.
109 * So the only cases to check are when one of z or w is
110 * infinite.
111 */
112 recalc = 0;
113 i = testinfl(a);
114 j = testinfl(b);
115 if (i | j) { /* z is infinite */
116 /* "factor out" infinity */
117 a = i;
118 b = j;
119 recalc = 1;
120 }
121 i = testinfl(c);
122 j = testinfl(d);
123 if (i | j) { /* w is infinite */
124 /* "factor out" infinity */
125 c = i;
126 d = j;
127 recalc = 1;
128 }
129 if (recalc) {
130 x = inf.f * (a * c - b * d);
131 y = inf.f * (a * d + b * c);
132 }
133 }
134
135 /*
136 * The following is equivalent to
137 *
138 * return x + I * y;
139 */
140 ((long double *)&v)[0] = x;
141 ((long double *)&v)[1] = y;
142 return (v);
143 }
144