contrib/openbsd_libm/src/s_csqrtl.c

   1 /*      $OpenBSD: s_csqrtl.c,v 1.2 2011/07/20 19:28:33 martynas Exp $   */
   2
   3 /*
   4  * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
   5  *
   6  * Permission to use, copy, modify, and distribute this software for any
   7  * purpose with or without fee is hereby granted, provided that the above
   8  * copyright notice and this permission notice appear in all copies.
   9  *
  10  * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
  11  * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
  12  * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
  13  * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
  14  * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
  15  * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
  16  * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
  17  */
  18
  19 /*                                                      csqrtl()
  20  *
  21  *      Complex square root
  22  *
  23  *
  24  *
  25  * SYNOPSIS:
  26  *
  27  * long double complex csqrtl();
  28  * long double complex z, w;
  29  *
  30  * w = csqrtl( z );
  31  *
  32  *
  33  *
  34  * DESCRIPTION:
  35  *
  36  *
  37  * If z = x + iy,  r = |z|, then
  38  *
  39  *                       1/2
  40  * Re w  =  [ (r + x)/2 ]   ,
  41  *
  42  *                       1/2
  43  * Im w  =  [ (r - x)/2 ]   .
  44  *
  45  * Cancellation error in r-x or r+x is avoided by using the
  46  * identity  2 Re w Im w  =  y.
  47  *
  48  * Note that -w is also a square root of z.  The root chosen
  49  * is always in the right half plane and Im w has the same sign as y.
  50  *
  51  *
  52  *
  53  * ACCURACY:
  54  *
  55  *                      Relative error:
  56  * arithmetic   domain     # trials      peak         rms
  57  *    IEEE      -10,+10     500000      1.1e-19     3.0e-20
  58  *
  59  */
  60
  61 #include <complex.h>
  62 #include <math.h>
  63
  64 long double complex
  65 csqrtl(long double complex z)
  66 {
  67         long double complex w;
  68         long double x, y, r, t, scale;
  69
  70         x = creall(z);
  71         y = cimagl(z);
  72
  73         if (y == 0.0L) {
  74                 if (x < 0.0L) {
  75                         w = 0.0L + sqrtl(-x) * I;
  76                         return (w);
  77                 }
  78                 else {
  79                         w = sqrtl(x) + 0.0L * I;
  80                         return (w);
  81                 }
  82         }
  83
  84         if (x == 0.0L) {
  85                 r = fabsl(y);
  86                 r = sqrtl(0.5L * r);
  87                 if (y > 0.0L)
  88                         w = r + r * I;
  89                 else
  90                         w = r - r * I;
  91                 return (w);
  92         }
  93
  94         /* Rescale to avoid internal overflow or underflow.  */
  95         if ((fabsl(x) > 4.0L) || (fabsl(y) > 4.0L)) {
  96                 x *= 0.25L;
  97                 y *= 0.25L;
  98                 scale = 2.0L;
  99         }
 100         else {
 101 #if 1
 102                 x *= 7.3786976294838206464e19;  /* 2^66 */
 103                 y *= 7.3786976294838206464e19;
 104                 scale = 1.16415321826934814453125e-10;  /* 2^-33 */
 105 #else
 106                 x *= 4.0L;
 107                 y *= 4.0L;
 108                 scale = 0.5L;
 109 #endif
 110         }
 111         w = x + y * I;
 112         r = cabsl(w);
 113         if (x > 0) {
 114                 t = sqrtl(0.5L * r + 0.5L * x);
 115                 r = scale * fabsl((0.5L * y) / t);
 116                 t *= scale;
 117         }
 118         else {
 119                 r = sqrtl(0.5L * r - 0.5L * x);
 120                 t = scale * fabsl((0.5L * y) / r);
 121                 r *= scale;
 122         }
 123         if (y < 0)
 124                 w = t - r * I;
 125         else
 126                 w = t + r * I;
 127         return (w);
 128 }