Print this page
11210 libm should be cstyle(1ONBLD) clean
*** 20,37 ****
*/
/*
* Copyright 2011 Nexenta Systems, Inc. All rights reserved.
*/
/*
* Copyright 2006 Sun Microsystems, Inc. All rights reserved.
* Use is subject to license terms.
*/
#pragma weak __csqrt = csqrt
! /* INDENT OFF */
/*
* dcomplex csqrt(dcomplex z);
*
* 2 2 2
* Let w=r+i*s = sqrt(x+iy). Then (r + i s) = r - s + i 2sr = x + i y.
--- 20,38 ----
*/
/*
* Copyright 2011 Nexenta Systems, Inc. All rights reserved.
*/
+
/*
* Copyright 2006 Sun Microsystems, Inc. All rights reserved.
* Use is subject to license terms.
*/
#pragma weak __csqrt = csqrt
!
/*
* dcomplex csqrt(dcomplex z);
*
* 2 2 2
* Let w=r+i*s = sqrt(x+iy). Then (r + i s) = r - s + i 2sr = x + i y.
*** 98,125 ****
* csqrt(-inf+ i NaN) = NaN +-i inf
* csqrt(+inf+ i NaN) = inf + i NaN
* csqrt(NaN + i y ) = NaN + i NaN for finite y
* csqrt(NaN + i NaN) = NaN + i NaN
*/
- /* INDENT ON */
#include "libm.h" /* fabs/sqrt */
#include "complex_wrapper.h"
! /* INDENT OFF */
! static const double
! two300 = 2.03703597633448608627e+90,
twom300 = 4.90909346529772655310e-91,
two599 = 2.07475778444049647926e+180,
twom601 = 1.20495993255144205887e-181,
two = 2.0,
zero = 0.0,
half = 0.5;
! /* INDENT ON */
dcomplex
! csqrt(dcomplex z) {
dcomplex ans;
double x, y, t, ax, ay;
int n, ix, iy, hx, hy, lx, ly;
x = D_RE(z);
--- 99,124 ----
* csqrt(-inf+ i NaN) = NaN +-i inf
* csqrt(+inf+ i NaN) = inf + i NaN
* csqrt(NaN + i y ) = NaN + i NaN for finite y
* csqrt(NaN + i NaN) = NaN + i NaN
*/
#include "libm.h" /* fabs/sqrt */
#include "complex_wrapper.h"
! static const double two300 = 2.03703597633448608627e+90,
twom300 = 4.90909346529772655310e-91,
two599 = 2.07475778444049647926e+180,
twom601 = 1.20495993255144205887e-181,
two = 2.0,
zero = 0.0,
half = 0.5;
!
dcomplex
! csqrt(dcomplex z)
! {
dcomplex ans;
double x, y, t, ax, ay;
int n, ix, iy, hx, hy, lx, ly;
x = D_RE(z);
*** 130,184 ****
ly = LO_WORD(y);
ix = hx & 0x7fffffff;
iy = hy & 0x7fffffff;
ay = fabs(y);
ax = fabs(x);
if (ix >= 0x7ff00000 || iy >= 0x7ff00000) {
/* x or y is Inf or NaN */
! if (ISINF(iy, ly))
D_IM(ans) = D_RE(ans) = ay;
! else if (ISINF(ix, lx)) {
if (hx > 0) {
D_RE(ans) = ax;
D_IM(ans) = ay * zero;
} else {
D_RE(ans) = ay * zero;
D_IM(ans) = ax;
}
! } else
D_IM(ans) = D_RE(ans) = ax + ay;
} else if ((iy | ly) == 0) { /* y = 0 */
if (hx >= 0) {
D_RE(ans) = sqrt(ax);
D_IM(ans) = zero;
} else {
D_IM(ans) = sqrt(ax);
D_RE(ans) = zero;
}
} else if (ix >= iy) {
n = (ix - iy) >> 20;
if (n >= 30) { /* x >> y or y=0 */
t = sqrt(ax);
} else if (ix >= 0x5f300000) { /* x > 2**500 */
ax *= twom601;
y *= twom601;
t = two300 * sqrt(ax + sqrt(ax * ax + y * y));
} else if (iy < 0x20b00000) { /* y < 2**-500 */
ax *= two599;
y *= two599;
t = twom300 * sqrt(ax + sqrt(ax * ax + y * y));
! } else
t = sqrt(half * (ax + sqrt(ax * ax + ay * ay)));
if (hx >= 0) {
D_RE(ans) = t;
D_IM(ans) = ay / (t + t);
} else {
D_IM(ans) = t;
D_RE(ans) = ay / (t + t);
}
} else {
n = (iy - ix) >> 20;
if (n >= 30) { /* y >> x */
if (n >= 60)
t = sqrt(half * ay);
else if (iy >= 0x7fe00000)
t = sqrt(half * ay + half * ax);
--- 129,189 ----
ly = LO_WORD(y);
ix = hx & 0x7fffffff;
iy = hy & 0x7fffffff;
ay = fabs(y);
ax = fabs(x);
+
if (ix >= 0x7ff00000 || iy >= 0x7ff00000) {
/* x or y is Inf or NaN */
! if (ISINF(iy, ly)) {
D_IM(ans) = D_RE(ans) = ay;
! } else if (ISINF(ix, lx)) {
if (hx > 0) {
D_RE(ans) = ax;
D_IM(ans) = ay * zero;
} else {
D_RE(ans) = ay * zero;
D_IM(ans) = ax;
}
! } else {
D_IM(ans) = D_RE(ans) = ax + ay;
+ }
} else if ((iy | ly) == 0) { /* y = 0 */
if (hx >= 0) {
D_RE(ans) = sqrt(ax);
D_IM(ans) = zero;
} else {
D_IM(ans) = sqrt(ax);
D_RE(ans) = zero;
}
} else if (ix >= iy) {
n = (ix - iy) >> 20;
+
if (n >= 30) { /* x >> y or y=0 */
t = sqrt(ax);
} else if (ix >= 0x5f300000) { /* x > 2**500 */
ax *= twom601;
y *= twom601;
t = two300 * sqrt(ax + sqrt(ax * ax + y * y));
} else if (iy < 0x20b00000) { /* y < 2**-500 */
ax *= two599;
y *= two599;
t = twom300 * sqrt(ax + sqrt(ax * ax + y * y));
! } else {
t = sqrt(half * (ax + sqrt(ax * ax + ay * ay)));
+ }
+
if (hx >= 0) {
D_RE(ans) = t;
D_IM(ans) = ay / (t + t);
} else {
D_IM(ans) = t;
D_RE(ans) = ay / (t + t);
}
} else {
n = (iy - ix) >> 20;
+
if (n >= 30) { /* y >> x */
if (n >= 60)
t = sqrt(half * ay);
else if (iy >= 0x7fe00000)
t = sqrt(half * ay + half * ax);
*** 192,210 ****
t = two300 * sqrt(ax + sqrt(ax * ax + y * y));
} else if (ix < 0x20b00000) { /* x < 2**-500 */
ax *= two599;
y *= two599;
t = twom300 * sqrt(ax + sqrt(ax * ax + y * y));
! } else
t = sqrt(half * (ax + sqrt(ax * ax + ay * ay)));
if (hx >= 0) {
D_RE(ans) = t;
D_IM(ans) = ay / (t + t);
} else {
D_IM(ans) = t;
D_RE(ans) = ay / (t + t);
}
}
if (hy < 0)
D_IM(ans) = -D_IM(ans);
return (ans);
}
--- 197,219 ----
t = two300 * sqrt(ax + sqrt(ax * ax + y * y));
} else if (ix < 0x20b00000) { /* x < 2**-500 */
ax *= two599;
y *= two599;
t = twom300 * sqrt(ax + sqrt(ax * ax + y * y));
! } else {
t = sqrt(half * (ax + sqrt(ax * ax + ay * ay)));
+ }
+
if (hx >= 0) {
D_RE(ans) = t;
D_IM(ans) = ay / (t + t);
} else {
D_IM(ans) = t;
D_RE(ans) = ay / (t + t);
}
}
+
if (hy < 0)
D_IM(ans) = -D_IM(ans);
+
return (ans);
}