Cross Reference: _F_cplx_div.c
xref
: /
osnet-11
/
usr
/
src
/
lib
/
libc
/
sparc
/
fp
/
_F_cplx_div.c
Home
History
Annotate
Line#
Navigate
Download
Search
only in
./
2
N/A
/*
2
N/A
* CDDL HEADER START
2
N/A
*
2
N/A
* The contents of this file are subject to the terms of the
2
N/A
* Common Development and Distribution License, Version 1.0 only
2
N/A
* (the "License"). You may not use this file except in compliance
2
N/A
* with the License.
2
N/A
*
2
N/A
* You can obtain a copy of the license at
usr
/
src
/
OPENSOLARIS.LICENSE
2
N/A
* or
http://www.opensolaris.org/os/licensing
.
2
N/A
* See the License for the specific language governing permissions
2
N/A
* and limitations under the License.
2
N/A
*
2
N/A
* When distributing Covered Code, include this CDDL HEADER in each
2
N/A
* file and include the License file at
usr
/
src
/
OPENSOLARIS.LICENSE
.
2
N/A
* If applicable, add the following below this CDDL HEADER, with the
2
N/A
* fields enclosed by brackets "[]" replaced with your own identifying
2
N/A
* information: Portions Copyright [yyyy] [name of copyright owner]
2
N/A
*
2
N/A
* CDDL HEADER END
2
N/A
*/
2
N/A
/*
2
N/A
* Copyright 2003 Sun Microsystems, Inc. All rights reserved.
2
N/A
* Use is subject to license terms.
2
N/A
*/
2
N/A
2
N/A
#
pragma
ident
"%Z%%M% %I% %E% SMI"
2
N/A
2
N/A
/*
2
N/A
* _F_cplx_div(z, w) returns z / w with infinities handled according
2
N/A
* to C99.
2
N/A
*
2
N/A
* If z and w are both finite and w is nonzero, _F_cplx_div(z, w)
2
N/A
* delivers the complex quotient q according to the usual formula:
2
N/A
* let a = Re(z), b = Im(z), c = Re(w), and d = Im(w); then q = x +
2
N/A
* I * y where x = (a * c + b * d) / r and y = (b * c - a * d) / r
2
N/A
* with r = c * c + d * d. This implementation computes intermediate
2
N/A
* results in double precision to avoid premature underflow or over-
2
N/A
* flow.
2
N/A
*
2
N/A
* If z is neither NaN nor zero and w is zero, or if z is infinite
2
N/A
* and w is finite and nonzero, _F_cplx_div delivers an infinite
2
N/A
* result. If z is finite and w is infinite, _F_cplx_div delivers
2
N/A
* a zero result.
2
N/A
*
2
N/A
* If z and w are both zero or both infinite, or if either z or w is
2
N/A
* a complex NaN, _F_cplx_div delivers NaN + I * NaN. C99 doesn't
2
N/A
* specify these cases.
2
N/A
*
2
N/A
* This implementation can raise spurious invalid operation, inexact,
2
N/A
* and division-by-zero exceptions. C99 allows this.
2
N/A
*
2
N/A
* Warning: Do not attempt to "optimize" this code by removing multi-
2
N/A
* plications by zero.
2
N/A
*/
2
N/A
2
N/A
#
if
!
defined
(
sparc
) && !
defined
(
__sparc
)
2
N/A
#
error
This
code
is
for
SPARC
only
2
N/A
#
endif
2
N/A
2
N/A
static
union
{
2
N/A
int
i[
2
];
2
N/A
double
d;
2
N/A
}
inf
= {
2
N/A
0x7ff00000
, 0
2
N/A
};
2
N/A
2
N/A
/*
2
N/A
* Return +1 if x is +Inf, -1 if x is -Inf, and 0 otherwise
2
N/A
*/
2
N/A
static
int
2
N/A
testinff
(
float
x)
2
N/A
{
2
N/A
union
{
2
N/A
int
i;
2
N/A
float
f;
2
N/A
}
xx
;
2
N/A
2
N/A
xx
.f = x;
2
N/A
return
((((
xx
.i <<
1
) -
0xff000000
) == 0)? (
1
| (
xx
.i >>
31
)) : 0);
2
N/A
}
2
N/A
2
N/A
float
_Complex
2
N/A
_F_cplx_div
(
float
_Complex
z,
float
_Complex
w)
2
N/A
{
2
N/A
float
_Complex
v;
2
N/A
union
{
2
N/A
int
i;
2
N/A
float
f;
2
N/A
}
cc
,
dd
;
2
N/A
float
a, b, c, d;
2
N/A
double
r, x, y;
2
N/A
int
i, j,
recalc
;
2
N/A
2
N/A
/*
2
N/A
* The following is equivalent to
2
N/A
*
2
N/A
* a = crealf(z); b = cimagf(z);
2
N/A
* c = crealf(w); d = cimagf(w);
2
N/A
*/
2
N/A
a = ((
float
*)&z)[0];
2
N/A
b = ((
float
*)&z)[
1
];
2
N/A
c = ((
float
*)&w)[0];
2
N/A
d = ((
float
*)&w)[
1
];
2
N/A
2
N/A
r = (
double
)c * c + (
double
)d * d;
2
N/A
2
N/A
if
(r ==
0.0
) {
2
N/A
/* w is zero; multiply z by 1/Re(w) - I * Im(w) */
2
N/A
c =
1.0f
/ c;
2
N/A
i =
testinff
(a);
2
N/A
j =
testinff
(b);
2
N/A
if
(i | j) {
/* z is infinite */
2
N/A
a = i;
2
N/A
b = j;
2
N/A
}
2
N/A
((
float
*)&v)[0] = a * c + b * d;
2
N/A
((
float
*)&v)[
1
] = b * c - a * d;
2
N/A
return
(v);
2
N/A
}
2
N/A
2
N/A
r =
1.0
/ r;
2
N/A
x = ((
double
)a * c + (
double
)b * d) * r;
2
N/A
y = ((
double
)b * c - (
double
)a * d) * r;
2
N/A
2
N/A
if
(x != x && y != y) {
2
N/A
/*
2
N/A
* Both x and y are NaN, so z and w can't both be finite
2
N/A
* and nonzero. Since we handled the case w = 0 above,
2
N/A
* the only cases to check here are when one of z or w
2
N/A
* is infinite.
2
N/A
*/
2
N/A
r =
1.0
;
2
N/A
recalc
= 0;
2
N/A
i =
testinff
(a);
2
N/A
j =
testinff
(b);
2
N/A
if
(i | j) {
/* z is infinite */
2
N/A
/* "factor out" infinity */
2
N/A
a = i;
2
N/A
b = j;
2
N/A
r =
inf
.d;
2
N/A
recalc
=
1
;
2
N/A
}
2
N/A
i =
testinff
(c);
2
N/A
j =
testinff
(d);
2
N/A
if
(i | j) {
/* w is infinite */
2
N/A
/*
2
N/A
* "factor out" infinity, being careful to preserve
2
N/A
* signs of finite values
2
N/A
*/
2
N/A
cc
.f = c;
2
N/A
dd
.f = d;
2
N/A
c = i? i : ((
cc
.i < 0)? -
0.0f
:
0.0f
);
2
N/A
d = j? j : ((
dd
.i < 0)? -
0.0f
:
0.0f
);
2
N/A
r *=
0.0
;
2
N/A
recalc
=
1
;
2
N/A
}
2
N/A
if
(
recalc
) {
2
N/A
x = ((
double
)a * c + (
double
)b * d) * r;
2
N/A
y = ((
double
)b * c - (
double
)a * d) * r;
2
N/A
}
2
N/A
}
2
N/A
2
N/A
/*
2
N/A
* The following is equivalent to
2
N/A
*
2
N/A
* return x + I * y;
2
N/A
*/
2
N/A
((
float
*)&v)[0] = (
float
)x;
2
N/A
((
float
*)&v)[
1
] = (
float
)y;
2
N/A
return
(v);
2
N/A
}