csqrt.c revision ddc0e0b53c661f6e439e3b7072b3ef353eadb4af
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * CDDL HEADER START
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * The contents of this file are subject to the terms of the
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Common Development and Distribution License (the "License").
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * You may not use this file except in compliance with the License.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * or http://www.opensolaris.org/os/licensing.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * See the License for the specific language governing permissions
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * and limitations under the License.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * When distributing Covered Code, include this CDDL HEADER in each
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * If applicable, add the following below this CDDL HEADER, with the
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * fields enclosed by brackets "[]" replaced with your own identifying
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * information: Portions Copyright [yyyy] [name of copyright owner]
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * CDDL HEADER END
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Use is subject to license terms.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/* INDENT OFF */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * dcomplex csqrt(dcomplex z);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Let w=r+i*s = sqrt(x+iy). Then (r + i s) = r - s + i 2sr = x + i y.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Hence x = r*r-s*s, y = 2sr.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Note that x*x+y*y = (s*s+r*r)**2. Thus, we have
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * (1) r + s = \/ x + y ,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * (2) r - s = x
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * (3) 2sr = y.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Perform (1)-(2) and (1)+(2), we obtain
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * (4) 2 r = hypot(x,y)+x,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * (5) 2*s = hypot(x,y)-x
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * where hypot(x,y) = \/ x + y .
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * In order to avoid numerical cancellation, we use formula (4) for
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * positive x, and (5) for negative x. The other component is then
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * computed by formula (3).
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * ------------------
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * (assume x and y are of medium size, i.e., no over/underflow in squaring)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * If x >=0 then
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 2 \/ x + y + x y
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * r = ---------------------, s = -------; (6)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * (note that we choose sign(s) = sign(y) to force r >=0).
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 2 \/ x + y - x y
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * s = ---------------------, r = -------; (7)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * One may use the polar coordinate of a complex number to justify the
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * following exception cases:
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * EXCEPTION CASES (conform to ISO/IEC 9899:1999(E)):
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * csqrt(+-0+ i 0 ) = 0 + i 0
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * csqrt( x + i inf ) = inf + i inf for all x (including NaN)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * csqrt( x + i NaN ) = NaN + i NaN with invalid for finite x
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * csqrt(-inf+ iy ) = 0 + i inf for finite positive-signed y
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * csqrt(+inf+ iy ) = inf + i 0 for finite positive-signed y
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * csqrt(-inf+ i NaN) = NaN +-i inf
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * csqrt(+inf+ i NaN) = inf + i NaN
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * csqrt(NaN + i y ) = NaN + i NaN for finite y
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * csqrt(NaN + i NaN) = NaN + i NaN
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/* INDENT ON */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/* INDENT OFF */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisstatic const double
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/* INDENT ON */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (ix >= 0x7ff00000 || iy >= 0x7ff00000) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* x or y is Inf or NaN */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (ix >= 0x5f300000) { /* x > 2**500 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis t = two300 * sqrt(ax + sqrt(ax * ax + y * y));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (iy < 0x20b00000) { /* y < 2**-500 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis t = twom300 * sqrt(ax + sqrt(ax * ax + y * y));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis t = sqrt(half * (ax + sqrt(ax * ax + ay * ay)));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (iy >= 0x5f300000) { /* y > 2**500 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis t = two300 * sqrt(ax + sqrt(ax * ax + y * y));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (ix < 0x20b00000) { /* x < 2**-500 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis t = twom300 * sqrt(ax + sqrt(ax * ax + y * y));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis t = sqrt(half * (ax + sqrt(ax * ax + ay * ay)));