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