gonzui

    1: /*      $OpenBSD: s_catan.c,v 1.2 2011/07/08 19:25:31 martynas 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: /* LINTLIBRARY */
   19: 
   20: /*                                                      catan()
   21:  *
   22:  *      Complex circular arc tangent
   23:  *
   24:  *
   25:  *
   26:  * SYNOPSIS:
   27:  *
   28:  * double complex catan();
   29:  * double complex z, w;
   30:  *
   31:  * w = catan (z);
   32:  *
   33:  *
   34:  *
   35:  * DESCRIPTION:
   36:  *
   37:  * If
   38:  *     z = x + iy,
   39:  *
   40:  * then
   41:  *          1       (    2x     )
   42:  * Re w  =  - arctan(-----------)  +  k PI
   43:  *          2       (     2    2)
   44:  *                  (1 - x  - y )
   45:  *
   46:  *               ( 2         2)
   47:  *          1    (x  +  (y+1) )
   48:  * Im w  =  - log(------------)
   49:  *          4    ( 2         2)
   50:  *               (x  +  (y-1) )
   51:  *
   52:  * Where k is an arbitrary integer.
   53:  *
   54:  * catan(z) = -i catanh(iz).
   55:  *
   56:  * ACCURACY:
   57:  *
   58:  *                      Relative error:
   59:  * arithmetic   domain     # trials      peak         rms
   60:  *    DEC       -10,+10      5900       1.3e-16     7.8e-18
   61:  *    IEEE      -10,+10     30000       2.3e-15     8.5e-17
   62:  * The check catan( ctan(z) )  =  z, with |x| and |y| < PI/2,
   63:  * had peak relative error 1.5e-16, rms relative error
   64:  * 2.9e-17.  See also clog().
   65:  */
   66: 
   67: #include <sys/cdefs.h>
   68: #include <complex.h>
   69: #include <float.h>
   70: #include <math.h>
   71: 
   72: #define MAXNUM 1.0e308
   73: 
   74: static const double DP1 = 3.14159265160560607910E0;
   75: static const double DP2 = 1.98418714791870343106E-9;
   76: static const double DP3 = 1.14423774522196636802E-17;
   77: 
   78: static double
   79: _redupi(double x)
   80: {
   81:         double t;
   82:         long i;
   83: 
   84:         t = x/M_PI;
   85:         if(t >= 0.0)
   86:                 t += 0.5;
   87:         else
   88:                 t -= 0.5;
   89: 
   90:         i = t; /* the multiple */
   91:         t = i;
   92:         t = ((x - t * DP1) - t * DP2) - t * DP3;
   93:         return (t);
   94: }
   95: 
   96: double complex
   97: catan(double complex z)
   98: {
   99:         double complex w;
  100:         double a, t, x, x2, y;
  101: 
  102:         x = creal (z);
  103:         y = cimag (z);
  104: 
  105:         if ((x == 0.0) && (y > 1.0))
  106:                 goto ovrf;
  107: 
  108:         x2 = x * x;
  109:         a = 1.0 - x2 - (y * y);
  110:         if (a == 0.0)
  111:                 goto ovrf;
  112: 
  113:         t = 0.5 * atan2 (2.0 * x, a);
  114:         w = _redupi (t);
  115: 
  116:         t = y - 1.0;
  117:         a = x2 + (t * t);
  118:         if (a == 0.0)
  119:         goto ovrf;
  120: 
  121:         t = y + 1.0;
  122:         a = (x2 + (t * t))/a;
  123:         w = w + (0.25 * log (a)) * I;
  124:         return (w);
  125: 
  126: ovrf:
  127:         /*mtherr ("catan", OVERFLOW);*/
  128:         w = MAXNUM + MAXNUM * I;
  129:         return (w);
  130: }
  131: 
  132: #if     LDBL_MANT_DIG == 53
  133: #ifdef  lint
  134: /* PROTOLIB1 */
  135: long double complex catanl(long double complex);
  136: #else   /* lint */
  137: __weak_alias(catanl, catan);
  138: #endif  /* lint */
  139: #endif  /* LDBL_MANT_DIG == 53 */