/*
* 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
* 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.
*/
/* 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
* ------------------
*
*
* 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:
*
* 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 "complex_wrapper.h"
/* INDENT OFF */
static const double
/* INDENT ON */
x = D_RE(z);
y = D_IM(z);
/* x or y is Inf or NaN */
if (hx > 0) {
} else {
}
} else
if (hx >= 0) {
} else {
}
if (n >= 30) { /* x >> y or y=0 */
y *= twom601;
y *= two599;
} else
if (hx >= 0) {
} else {
}
} else {
if (n >= 30) { /* y >> x */
if (n >= 60)
else if (iy >= 0x7fe00000)
else if (ix <= 0x00100000)
else
y *= twom601;
y *= two599;
} else
if (hx >= 0) {
} else {
}
}
if (hy < 0)
return (ans);
}