newlib/libm/complex/csqrt.c

   1 /* $NetBSD: csqrt.c,v 1.1 2007/08/20 16:01:37 drochner Exp $ */
   2
   3 /*-
   4  * Copyright (c) 2007 The NetBSD Foundation, Inc.
   5  * All rights reserved.
   6  *
   7  * This code is derived from software written by Stephen L. Moshier.
   8  * It is redistributed by the NetBSD Foundation by permission of the author.
   9  *
  10  * Redistribution and use in source and binary forms, with or without
  11  * modification, are permitted provided that the following conditions
  12  * are met:
  13  * 1. Redistributions of source code must retain the above copyright
  14  *    notice, this list of conditions and the following disclaimer.
  15  * 2. Redistributions in binary form must reproduce the above copyright
  16  *    notice, this list of conditions and the following disclaimer in the
  17  *    documentation and/or other materials provided with the distribution.
  18  *
  19  * THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC. AND CONTRIBUTORS
  20  * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
  21  * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
  22  * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR CONTRIBUTORS
  23  * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
  24  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
  25  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  26  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
  27  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  28  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  29  * POSSIBILITY OF SUCH DAMAGE.
  30  *
  31  * imported and modified include for newlib 2010/10/03
  32  * Marco Atzeri <marco_atzeri@yahoo.it>
  33  */
  34
  35 /*
  36 FUNCTION
  37         <<csqrt>>, <<csqrtf>>---complex square root
  38
  39 INDEX
  40         csqrt
  41 INDEX
  42         csqrtf
  43
  44 SYNOPSIS
  45        #include <complex.h>
  46        double complex csqrt(double complex <[z]>);
  47        float complex csqrtf(float complex <[z]>);
  48
  49
  50 DESCRIPTION
  51         These functions compute the complex square root of <[z]>, with
  52         a branch cut along the negative real axis.
  53
  54         <<csqrtf>> is identical to <<csqrt>>, except that it performs
  55         its calculations on <<floats complex>>.
  56
  57 RETURNS
  58         The csqrt functions return the complex square root value, in
  59         the range of the right halfplane (including the imaginary axis).
  60
  61 PORTABILITY
  62         <<csqrt>> and <<csqrtf>> are ISO C99
  63
  64 QUICKREF
  65         <<csqrt>> and <<csqrtf>> are ISO C99
  66
  67 */
  68
  69
  70 #include <complex.h>
  71 #include <math.h>
  72
  73 double complex
  74 csqrt(double complex z)
  75 {
  76         double complex w;
  77         double x, y, r, t, scale;
  78
  79         x = creal (z);
  80         y = cimag (z);
  81
  82         if (y == 0.0) {
  83                 if (x == 0.0) {
  84                         w = 0.0 + y * I;
  85                 } else {
  86                         r = fabs(x);
  87                         r = sqrt(r);
  88                         if (x < 0.0) {
  89                                 w = 0.0 + r * I;
  90                         } else {
  91                                 w = r + y * I;
  92                         }
  93                 }
  94                 return w;
  95         }
  96         if (x == 0.0) {
  97                 r = fabs(y);
  98                 r = sqrt(0.5 * r);
  99                 if (y > 0)
 100                         w = r + r * I;
 101                 else
 102                         w = r - r * I;
 103                 return w;
 104         }
 105         /* Rescale to avoid internal overflow or underflow.  */
 106         if ((fabs(x) > 4.0) || (fabs(y) > 4.0)) {
 107                 x *= 0.25;
 108                 y *= 0.25;
 109                 scale = 2.0;
 110         } else {
 111 #if 1
 112                 x *= 1.8014398509481984e16;  /* 2^54 */
 113                 y *= 1.8014398509481984e16;
 114                 scale = 7.450580596923828125e-9; /* 2^-27 */
 115 #else
 116                 x *= 4.0;
 117                 y *= 4.0;
 118                 scale = 0.5;
 119 #endif
 120         }
 121         w = x + y * I;
 122         r = cabs(w);
 123         if (x > 0) {
 124                 t = sqrt(0.5 * r + 0.5 * x);
 125                 r = scale * fabs((0.5 * y) / t );
 126                 t *= scale;
 127         } else {
 128                 r = sqrt(0.5 * r - 0.5 * x);
 129                 t = scale * fabs((0.5 * y) / r);
 130                 r *= scale;
 131         }
 132         if (y < 0)
 133                 w = t - r * I;
 134         else
 135                 w = t + r * I;
 136         return w;
 137 }