Cross Reference: _Q_cplx_div.c
xref
: /
osnet-11
/
usr
/
src
/
lib
/
libc
/
sparc
/
fp
/
_Q_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
* On SPARC V8, _Q_cplx_div(v, z, w) sets *v = *z / *w with infin-
2
N/A
* ities handling according to C99.
2
N/A
*
2
N/A
* On SPARC V9, _Q_cplx_div(z, w) returns *z / *w with infinities
2
N/A
* handled according to C99.
2
N/A
*
2
N/A
* If z and w are both finite and w is nonzero, _Q_cplx_div delivers
2
N/A
* the complex quotient q according to the usual formula: let a =
2
N/A
* Re(z), b = Im(z), c = Re(w), and d = Im(w); then q = x + I * y
2
N/A
* where x = (a * c + b * d) / r and y = (b * c - a * d) / r with
2
N/A
* r = c * c + d * d. This implementation scales to avoid premature
2
N/A
* underflow or overflow.
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, _Q_cplx_div delivers an infinite
2
N/A
* result. If z is finite and w is infinite, _Q_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, _Q_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 underflow, overflow, in-
2
N/A
* valid operation, inexact, and division-by-zero exceptions. C99
2
N/A
* allows this.
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[
4
];
2
N/A
long
double
q;
2
N/A
}
inf
= {
2
N/A
0x7fff0000
, 0, 0, 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
testinfl
(
long
double
x)
2
N/A
{
2
N/A
union
{
2
N/A
int
i[
4
];
2
N/A
long
double
q;
2
N/A
}
xx
;
2
N/A
2
N/A
xx
.q = x;
2
N/A
return
(((((
xx
.i[0] <<
1
) -
0xfffe0000
) |
xx
.i[
1
] |
xx
.i[
2
] |
xx
.i[
3
])
2
N/A
== 0)? (
1
| (
xx
.i[0] >>
31
)) : 0);
2
N/A
}
2
N/A
2
N/A
#
ifdef
__sparcv9
2
N/A
long
double
_Complex
2
N/A
_Q_cplx_div
(
const
long
double
_Complex
*z,
const
long
double
_Complex
*w)
2
N/A
{
2
N/A
long
double
_Complex
v;
2
N/A
#
else
2
N/A
void
2
N/A
_Q_cplx_div
(
long
double
_Complex
*v,
const
long
double
_Complex
*z,
2
N/A
const
long
double
_Complex
*w)
2
N/A
{
2
N/A
#
endif
2
N/A
union
{
2
N/A
int
i[
4
];
2
N/A
long
double
q;
2
N/A
}
aa
,
bb
,
cc
,
dd
,
ss
;
2
N/A
long
double
a, b, c, d, r;
2
N/A
int
ha
,
hb
,
hc
,
hd
,
hz
,
hw
,
hs
, i, j;
2
N/A
2
N/A
/*
2
N/A
* The following is equivalent to
2
N/A
*
2
N/A
* a = creall(*z); b = cimagl(*z);
2
N/A
* c = creall(*w); d = cimagl(*w);
2
N/A
*/
2
N/A
a = ((
long
double
*)z)[0];
2
N/A
b = ((
long
double
*)z)[
1
];
2
N/A
c = ((
long
double
*)w)[0];
2
N/A
d = ((
long
double
*)w)[
1
];
2
N/A
2
N/A
/* extract high-order words to estimate |z| and |w| */
2
N/A
aa
.q = a;
2
N/A
bb
.q = b;
2
N/A
ha
=
aa
.i[0] & ~
0x80000000
;
2
N/A
hb
=
bb
.i[0] & ~
0x80000000
;
2
N/A
hz
= (
ha
>
hb
)?
ha
:
hb
;
2
N/A
2
N/A
cc
.q = c;
2
N/A
dd
.q = d;
2
N/A
hc
=
cc
.i[0] & ~
0x80000000
;
2
N/A
hd
=
dd
.i[0] & ~
0x80000000
;
2
N/A
hw
= (
hc
>
hd
)?
hc
:
hd
;
2
N/A
2
N/A
/* check for special cases */
2
N/A
if
(
hw
>=
0x7fff0000
) {
/* w is inf or nan */
2
N/A
r =
0.0l
;
2
N/A
i =
testinfl
(c);
2
N/A
j =
testinfl
(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
c = i? i : ((
cc
.i[0] < 0)? -
0.0l
:
0.0l
);
2
N/A
d = j? j : ((
dd
.i[0] < 0)? -
0.0l
:
0.0l
);
2
N/A
if
(
hz
>=
0x7ffe0000
) {
2
N/A
/* scale to avoid overflow below */
2
N/A
c *=
0.5l
;
2
N/A
d *=
0.5l
;
2
N/A
}
2
N/A
}
2
N/A
goto
done
;
2
N/A
}
2
N/A
2
N/A
if
(
hw
== 0 && (
cc
.i[
1
] |
cc
.i[
2
] |
cc
.i[
3
] |
2
N/A
dd
.i[
1
] |
dd
.i[
2
] |
dd
.i[
3
]) == 0) {
2
N/A
/* w is zero; multiply z by 1/Re(w) - I * Im(w) */
2
N/A
r =
1.0l
;
2
N/A
c =
1.0l
/ c;
2
N/A
i =
testinfl
(a);
2
N/A
j =
testinfl
(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
goto
done
;
2
N/A
}
2
N/A
2
N/A
if
(
hz
>=
0x7fff0000
) {
/* z is inf or nan */
2
N/A
r =
1.0l
;
2
N/A
i =
testinfl
(a);
2
N/A
j =
testinfl
(b);
2
N/A
if
(i | j) {
/* z is infinite */
2
N/A
a = i;
2
N/A
b = j;
2
N/A
r =
inf
.q;
2
N/A
}
2
N/A
goto
done
;
2
N/A
}
2
N/A
2
N/A
/*
2
N/A
* Scale c and d to compute 1/|w|^2 and the real and imaginary
2
N/A
* parts of the quotient.
2
N/A
*/
2
N/A
hs
= (((
hw
>>
2
) -
hw
) +
0x6ffd7fff
) &
0xffff0000
;
2
N/A
if
(
hz
<
0x00ea0000
) {
/* |z| < 2^-16149 */
2
N/A
if
(((
hw
-
0x3e380000
) | (
0x40e90000
-
hw
)) >= 0)
2
N/A
hs
= (((
0x40e90000
-
hw
) >>
1
) &
0xffff0000
)
2
N/A
+
0x3fff0000
;
2
N/A
}
2
N/A
ss
.i[0] =
hs
;
2
N/A
ss
.i[
1
] =
ss
.i[
2
] =
ss
.i[
3
] = 0;
2
N/A
2
N/A
c *=
ss
.q;
2
N/A
d *=
ss
.q;
2
N/A
r =
1.0l
/ (c * c + d * d);
2
N/A
2
N/A
c *=
ss
.q;
2
N/A
d *=
ss
.q;
2
N/A
2
N/A
done
:
2
N/A
#
ifdef
__sparcv9
2
N/A
((
long
double
*)&v)[0] = (a * c + b * d) * r;
2
N/A
((
long
double
*)&v)[
1
] = (b * c - a * d) * r;
2
N/A
return
(v);
2
N/A
#
else
2
N/A
((
long
double
*)v)[0] = (a * c + b * d) * r;
2
N/A
((
long
double
*)v)[
1
] = (b * c - a * d) * r;
2
N/A
#
endif
2
N/A
}