gonzui

    1: /*      $OpenBSD: s_catanf.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: /*                                                      catanf()
   19:  *
   20:  *      Complex circular arc tangent
   21:  *
   22:  *
   23:  *
   24:  * SYNOPSIS:
   25:  *
   26:  * float complex catanf();
   27:  * float complex z, w;
   28:  *
   29:  * w = catanf( z );
   30:  *
   31:  *
   32:  *
   33:  * DESCRIPTION:
   34:  *
   35:  * If
   36:  *     z = x + iy,
   37:  *
   38:  * then
   39:  *          1       (    2x     )
   40:  * Re w  =  - arctan(-----------)  +  k PI
   41:  *          2       (     2    2)
   42:  *                  (1 - x  - y )
   43:  *
   44:  *               ( 2         2)
   45:  *          1    (x  +  (y+1) )
   46:  * Im w  =  - log(------------)
   47:  *          4    ( 2         2)
   48:  *               (x  +  (y-1) )
   49:  *
   50:  * Where k is an arbitrary integer.
   51:  *
   52:  *
   53:  *
   54:  * ACCURACY:
   55:  *
   56:  *                      Relative error:
   57:  * arithmetic   domain     # trials      peak         rms
   58:  *    IEEE      -10,+10     30000        2.3e-6      5.2e-8
   59:  *
   60:  */
   61: 
   62: #include <complex.h>
   63: #include <math.h>
   64: 
   65: #define MAXNUMF 1.0e38F
   66: 
   67: static const double DP1 = 3.140625;
   68: static const double DP2 = 9.67502593994140625E-4;
   69: static const double DP3 = 1.509957990978376432E-7;
   70: 
   71: static float
   72: _redupif(float xx)
   73: {
   74:         float x, t;
   75:         long i;
   76: 
   77:         x = xx;
   78:         t = x/(float)M_PI;
   79:         if(t >= 0.0)
   80:                 t += 0.5;
   81:         else
   82:                 t -= 0.5;
   83: 
   84:         i = t; /* the multiple */
   85:         t = i;
   86:         t = ((x - t * DP1) - t * DP2) - t * DP3;
   87:         return(t);
   88: }
   89: 
   90: float complex
   91: catanf(float complex z)
   92: {
   93:         float complex w;
   94:         float a, t, x, x2, y;
   95: 
   96:         x = crealf(z);
   97:         y = cimagf(z);
   98: 
   99:         if((x == 0.0f) && (y > 1.0f))
  100:                 goto ovrf;
  101: 
  102:         x2 = x * x;
  103:         a = 1.0f - x2 - (y * y);
  104:         if (a == 0.0f)
  105:                 goto ovrf;
  106: 
  107:         t = 0.5f * atan2f(2.0f * x, a);
  108:         w = _redupif(t);
  109: 
  110:         t = y - 1.0f;
  111:         a = x2 + (t * t);
  112:         if(a == 0.0f)
  113:                 goto ovrf;
  114: 
  115:         t = y + 1.0f;
  116:         a = (x2 + (t * t))/a;
  117:         w = w + (0.25f * logf (a)) * I;
  118:         return (w);
  119: 
  120: ovrf:
  121:         /*mtherr( "catanf", OVERFLOW );*/
  122:         w = MAXNUMF + MAXNUMF * I;
  123:         return (w);
  124: }