Home | History | Annotate | Line # | Download | only in complex
csqrtl.c revision 1.2.2.2
      1  1.2.2.2  martin /*-
      2  1.2.2.2  martin  * Copyright (c) 2007-2008 David Schultz <das (at) FreeBSD.ORG>
      3  1.2.2.2  martin  * All rights reserved.
      4  1.2.2.2  martin  *
      5  1.2.2.2  martin  * Redistribution and use in source and binary forms, with or without
      6  1.2.2.2  martin  * modification, are permitted provided that the following conditions
      7  1.2.2.2  martin  * are met:
      8  1.2.2.2  martin  * 1. Redistributions of source code must retain the above copyright
      9  1.2.2.2  martin  *    notice, this list of conditions and the following disclaimer.
     10  1.2.2.2  martin  * 2. Redistributions in binary form must reproduce the above copyright
     11  1.2.2.2  martin  *    notice, this list of conditions and the following disclaimer in the
     12  1.2.2.2  martin  *    documentation and/or other materials provided with the distribution.
     13  1.2.2.2  martin  *
     14  1.2.2.2  martin  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
     15  1.2.2.2  martin  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
     16  1.2.2.2  martin  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
     17  1.2.2.2  martin  * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
     18  1.2.2.2  martin  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
     19  1.2.2.2  martin  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
     20  1.2.2.2  martin  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
     21  1.2.2.2  martin  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
     22  1.2.2.2  martin  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
     23  1.2.2.2  martin  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
     24  1.2.2.2  martin  * SUCH DAMAGE.
     25  1.2.2.2  martin  */
     26  1.2.2.2  martin 
     27  1.2.2.2  martin #include <sys/cdefs.h>
     28  1.2.2.2  martin #if 0
     29  1.2.2.2  martin __FBSDID("$FreeBSD: head/lib/msun/src/s_csqrtl.c 181402 2008-08-08 00:15:16Z das $");
     30  1.2.2.2  martin #else
     31  1.2.2.2  martin __RCSID("$NetBSD: csqrtl.c,v 1.2.2.2 2014/10/13 19:34:58 martin Exp $");
     32  1.2.2.2  martin #endif
     33  1.2.2.2  martin 
     34  1.2.2.2  martin #include "../src/namespace.h"
     35  1.2.2.2  martin #include <complex.h>
     36  1.2.2.2  martin #include <float.h>
     37  1.2.2.2  martin #include <math.h>
     38  1.2.2.2  martin #include <stdbool.h>
     39  1.2.2.2  martin 
     40  1.2.2.2  martin /*
     41  1.2.2.2  martin  * gcc doesn't implement complex multiplication or division correctly,
     42  1.2.2.2  martin  * so we need to handle infinities specially. We turn on this pragma to
     43  1.2.2.2  martin  * notify conforming c99 compilers that the fast-but-incorrect code that
     44  1.2.2.2  martin  * gcc generates is acceptable, since the special cases have already been
     45  1.2.2.2  martin  * handled.
     46  1.2.2.2  martin  */
     47  1.2.2.2  martin // #pragma	STDC CX_LIMITED_RANGE	ON
     48  1.2.2.2  martin 
     49  1.2.2.2  martin /* We risk spurious overflow for components >= LDBL_MAX / (1 + sqrt(2)). */
     50  1.2.2.2  martin #define	THRESH	(LDBL_MAX / 2.414213562373095048801688724209698L)
     51  1.2.2.2  martin 
     52  1.2.2.2  martin #define cpackl(r, i) ((r) + (i) * I)
     53  1.2.2.2  martin 
     54  1.2.2.2  martin long double complex
     55  1.2.2.2  martin csqrtl(long double complex z)
     56  1.2.2.2  martin {
     57  1.2.2.2  martin 	long double complex result;
     58  1.2.2.2  martin 	long double a, b;
     59  1.2.2.2  martin 	long double t;
     60  1.2.2.2  martin 	bool scale;
     61  1.2.2.2  martin 
     62  1.2.2.2  martin 	a = creall(z);
     63  1.2.2.2  martin 	b = cimagl(z);
     64  1.2.2.2  martin 
     65  1.2.2.2  martin 	/* Handle special cases. */
     66  1.2.2.2  martin 	if (z == 0.0L)
     67  1.2.2.2  martin 		return (cpackl(0.0L, b));
     68  1.2.2.2  martin 	if (isinf(b))
     69  1.2.2.2  martin 		return (cpackl(INFINITY, b));
     70  1.2.2.2  martin 	if (isnan(a)) {
     71  1.2.2.2  martin 		t = (b - b) / (b - b);	/* raise invalid if b is not a NaN */
     72  1.2.2.2  martin 		return (cpackl(a, t));	/* return NaN + NaN i */
     73  1.2.2.2  martin 	}
     74  1.2.2.2  martin 	if (isinf(a)) {
     75  1.2.2.2  martin 		/*
     76  1.2.2.2  martin 		 * csqrt(inf + NaN i)  = inf +  NaN i
     77  1.2.2.2  martin 		 * csqrt(inf + y i)    = inf +  0 i
     78  1.2.2.2  martin 		 * csqrt(-inf + NaN i) = NaN +- inf i
     79  1.2.2.2  martin 		 * csqrt(-inf + y i)   = 0   +  inf i
     80  1.2.2.2  martin 		 */
     81  1.2.2.2  martin 		if (signbit(a))
     82  1.2.2.2  martin 			return (cpackl(fabsl(b - b), copysignl(a, b)));
     83  1.2.2.2  martin 		else
     84  1.2.2.2  martin 			return (cpackl(a, copysignl(b - b, b)));
     85  1.2.2.2  martin 	}
     86  1.2.2.2  martin 	/*
     87  1.2.2.2  martin 	 * The remaining special case (b is NaN) is handled just fine by
     88  1.2.2.2  martin 	 * the normal code path below.
     89  1.2.2.2  martin 	 */
     90  1.2.2.2  martin 
     91  1.2.2.2  martin 	/* Scale to avoid overflow. */
     92  1.2.2.2  martin 	if (fabsl(a) >= THRESH || fabsl(b) >= THRESH) {
     93  1.2.2.2  martin 		a *= 0.25L;
     94  1.2.2.2  martin 		b *= 0.25L;
     95  1.2.2.2  martin 		scale = true;
     96  1.2.2.2  martin 	} else {
     97  1.2.2.2  martin 		scale = false;
     98  1.2.2.2  martin 	}
     99  1.2.2.2  martin 
    100  1.2.2.2  martin 	/* Algorithm 312, CACM vol 10, Oct 1967. */
    101  1.2.2.2  martin 	if (a >= 0L) {
    102  1.2.2.2  martin 		t = sqrtl((a + hypotl(a, b)) * 0.5L);
    103  1.2.2.2  martin 		result = cpackl(t, b / (2.0L * t));
    104  1.2.2.2  martin 	} else {
    105  1.2.2.2  martin 		t = sqrtl((-a + hypotl(a, b)) * 0.5L);
    106  1.2.2.2  martin 		result = cpackl(fabsl(b) / (2.0L * t), copysignl(t, b));
    107  1.2.2.2  martin 	}
    108  1.2.2.2  martin 
    109  1.2.2.2  martin 	/* Rescale. */
    110  1.2.2.2  martin 	if (scale)
    111  1.2.2.2  martin 		return (result * 2.0L);
    112  1.2.2.2  martin 	else
    113  1.2.2.2  martin 		return (result);
    114  1.2.2.2  martin }
    115