Home | History | Annotate | Download | only in src
      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 __FBSDID("$FreeBSD$");
     29 
     30 #include <complex.h>
     31 #include <float.h>
     32 #include <math.h>
     33 
     34 #include "math_private.h"
     35 
     36 /*
     37  * gcc doesn't implement complex multiplication or division correctly,
     38  * so we need to handle infinities specially. We turn on this pragma to
     39  * notify conforming c99 compilers that the fast-but-incorrect code that
     40  * gcc generates is acceptable, since the special cases have already been
     41  * handled.
     42  */
     43 #pragma	STDC CX_LIMITED_RANGE	ON
     44 
     45 /* We risk spurious overflow for components >= LDBL_MAX / (1 + sqrt(2)). */
     46 #define	THRESH	(LDBL_MAX / 2.414213562373095048801688724209698L)
     47 
     48 long double complex
     49 csqrtl(long double complex z)
     50 {
     51 	long double complex result;
     52 	long double a, b;
     53 	long double t;
     54 	int scale;
     55 
     56 	a = creall(z);
     57 	b = cimagl(z);
     58 
     59 	/* Handle special cases. */
     60 	if (z == 0)
     61 		return (cpackl(0, b));
     62 	if (isinf(b))
     63 		return (cpackl(INFINITY, b));
     64 	if (isnan(a)) {
     65 		t = (b - b) / (b - b);	/* raise invalid if b is not a NaN */
     66 		return (cpackl(a, t));	/* return NaN + NaN i */
     67 	}
     68 	if (isinf(a)) {
     69 		/*
     70 		 * csqrt(inf + NaN i)  = inf +  NaN i
     71 		 * csqrt(inf + y i)    = inf +  0 i
     72 		 * csqrt(-inf + NaN i) = NaN +- inf i
     73 		 * csqrt(-inf + y i)   = 0   +  inf i
     74 		 */
     75 		if (signbit(a))
     76 			return (cpackl(fabsl(b - b), copysignl(a, b)));
     77 		else
     78 			return (cpackl(a, copysignl(b - b, b)));
     79 	}
     80 	/*
     81 	 * The remaining special case (b is NaN) is handled just fine by
     82 	 * the normal code path below.
     83 	 */
     84 
     85 	/* Scale to avoid overflow. */
     86 	if (fabsl(a) >= THRESH || fabsl(b) >= THRESH) {
     87 		a *= 0.25;
     88 		b *= 0.25;
     89 		scale = 1;
     90 	} else {
     91 		scale = 0;
     92 	}
     93 
     94 	/* Algorithm 312, CACM vol 10, Oct 1967. */
     95 	if (a >= 0) {
     96 		t = sqrtl((a + hypotl(a, b)) * 0.5);
     97 		result = cpackl(t, b / (2 * t));
     98 	} else {
     99 		t = sqrtl((-a + hypotl(a, b)) * 0.5);
    100 		result = cpackl(fabsl(b) / (2 * t), copysignl(t, b));
    101 	}
    102 
    103 	/* Rescale. */
    104 	if (scale)
    105 		return (result * 2);
    106 	else
    107 		return (result);
    108 }
    109