gonzui

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