diff options
Diffstat (limited to 'usr/src/lib/libm/common/complex/csqrt.c')
| -rw-r--r-- | usr/src/lib/libm/common/complex/csqrt.c | 210 |
1 files changed, 210 insertions, 0 deletions
diff --git a/usr/src/lib/libm/common/complex/csqrt.c b/usr/src/lib/libm/common/complex/csqrt.c new file mode 100644 index 0000000000..1a00236677 --- /dev/null +++ b/usr/src/lib/libm/common/complex/csqrt.c @@ -0,0 +1,210 @@ +/* + * CDDL HEADER START + * + * The contents of this file are subject to the terms of the + * Common Development and Distribution License (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 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. + * + * Hence x = r*r-s*s, y = 2sr. + * + * Note that x*x+y*y = (s*s+r*r)**2. Thus, we have + * ________ + * 2 2 / 2 2 + * (1) r + s = \/ x + y , + * + * 2 2 + * (2) r - s = x + * + * (3) 2sr = y. + * + * Perform (1)-(2) and (1)+(2), we obtain + * + * 2 + * (4) 2 r = hypot(x,y)+x, + * + * 2 + * (5) 2*s = hypot(x,y)-x + * ________ + * / 2 2 + * where hypot(x,y) = \/ x + y . + * + * In order to avoid numerical cancellation, we use formula (4) for + * positive x, and (5) for negative x. The other component is then + * computed by formula (3). + * + * + * ALGORITHM + * ------------------ + * + * (assume x and y are of medium size, i.e., no over/underflow in squaring) + * + * If x >=0 then + * ________ + * / 2 2 + * 2 \/ x + y + x y + * r = ---------------------, s = -------; (6) + * 2 2 r + * + * (note that we choose sign(s) = sign(y) to force r >=0). + * Otherwise, + * ________ + * / 2 2 + * 2 \/ x + y - x y + * s = ---------------------, r = -------; (7) + * 2 2 s + * + * EXCEPTION: + * + * One may use the polar coordinate of a complex number to justify the + * following exception cases: + * + * EXCEPTION CASES (conform to ISO/IEC 9899:1999(E)): + * csqrt(+-0+ i 0 ) = 0 + i 0 + * csqrt( x + i inf ) = inf + i inf for all x (including NaN) + * csqrt( x + i NaN ) = NaN + i NaN with invalid for finite x + * csqrt(-inf+ iy ) = 0 + i inf for finite positive-signed y + * csqrt(+inf+ iy ) = inf + i 0 for finite positive-signed y + * 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); + y = D_IM(z); + hx = HI_WORD(x); + lx = LO_WORD(x); + hy = HI_WORD(y); + 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); + else if (ix <= 0x00100000) + t = half * sqrt(two * (ay + ax)); + else + t = sqrt(half * (ay + ax)); + } else if (iy >= 0x5f300000) { /* y > 2**500 */ + ax *= twom601; + y *= twom601; + 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); +} |