1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
|
/*
* CDDL HEADER START
*
* The contents of this file are subject to the terms of the
* Common Development and Distribution License, Version 1.0 only
* (the "License"). You may not use this file except in compliance
* with the License.
*
* You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
* or http://www.opensolaris.org/os/licensing.
* See the License for the specific language governing permissions
* and limitations under the License.
*
* When distributing Covered Code, include this CDDL HEADER in each
* file and include the License file at usr/src/OPENSOLARIS.LICENSE.
* If applicable, add the following below this CDDL HEADER, with the
* fields enclosed by brackets "[]" replaced with your own identifying
* information: Portions Copyright [yyyy] [name of copyright owner]
*
* CDDL HEADER END
*/
/*
* Copyright 2003 Sun Microsystems, Inc. All rights reserved.
* Use is subject to license terms.
*/
/*
* On SPARC V8, _Q_cplx_mul(v, z, w) sets *v = *z * *w with infinities
* handled according to C99.
*
* On SPARC V9, _Q_cplx_mul(z, w) returns *z * *w with infinities
* handled according to C99.
*
* If z and w are both finite, _Q_cplx_mul delivers the complex
* product according to the usual formula: let a = Re(z), b = Im(z),
* c = Re(w), and d = Im(w); then _Q_cplx_mul delivers x + I * y
* where x = a * c - b * d and y = a * d + b * c. Note that if both
* ac and bd overflow, then at least one of ad or bc must also over-
* flow, and vice versa, so that if one component of the product is
* NaN, the other is infinite. (Such a value is considered infinite
* according to C99.)
*
* If one of z or w is infinite and the other is either finite nonzero
* or infinite, _Q_cplx_mul delivers an infinite result. If one factor
* is infinite and the other is zero, _Q_cplx_mul delivers NaN + I * NaN.
* C99 doesn't specify the latter case.
*
* C99 also doesn't specify what should happen if either z or w is a
* complex NaN (i.e., neither finite nor infinite). This implementation
* delivers NaN + I * NaN in this case.
*
* This implementation can raise spurious underflow, overflow, invalid
* operation, and inexact exceptions. C99 allows this.
*/
#if !defined(sparc) && !defined(__sparc)
#error This code is for SPARC only
#endif
static union {
int i[4];
long double q;
} inf = {
0x7fff0000, 0, 0, 0
};
/*
* Return +1 if x is +Inf, -1 if x is -Inf, and 0 otherwise
*/
static int
testinfl(long double x)
{
union {
int i[4];
long double q;
} xx;
xx.q = x;
return (((((xx.i[0] << 1) - 0xfffe0000) | xx.i[1] | xx.i[2] | xx.i[3])
== 0)? (1 | (xx.i[0] >> 31)) : 0);
}
#ifdef __sparcv9
long double _Complex
_Q_cplx_mul(const long double _Complex *z, const long double _Complex *w)
{
long double _Complex v = 0;
#else
void
_Q_cplx_mul(long double _Complex *v, const long double _Complex *z,
const long double _Complex *w)
{
#endif
long double a, b, c, d, x, y;
int recalc, i, j;
/*
* The following is equivalent to
*
* a = creall(*z); b = cimagl(*z);
* c = creall(*w); d = cimagl(*w);
*/
a = ((long double *)z)[0];
b = ((long double *)z)[1];
c = ((long double *)w)[0];
d = ((long double *)w)[1];
x = a * c - b * d;
y = a * d + b * c;
if (x != x && y != y) {
/*
* Both x and y are NaN, so z and w can't both be finite.
* If at least one of z or w is a complex NaN, and neither
* is infinite, then we might as well deliver NaN + I * NaN.
* So the only cases to check are when one of z or w is
* infinite.
*/
recalc = 0;
i = testinfl(a);
j = testinfl(b);
if (i | j) { /* z is infinite */
/* "factor out" infinity */
a = i;
b = j;
recalc = 1;
}
i = testinfl(c);
j = testinfl(d);
if (i | j) { /* w is infinite */
/* "factor out" infinity */
c = i;
d = j;
recalc = 1;
}
if (recalc) {
x = inf.q * (a * c - b * d);
y = inf.q * (a * d + b * c);
}
}
#ifdef __sparcv9
((long double *)&v)[0] = x;
((long double *)&v)[1] = y;
return (v);
#else
((long double *)v)[0] = x;
((long double *)v)[1] = y;
#endif
}
|
