gonzui

    1: /*      $OpenBSD: s_ctan.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: /*                                                      ctan()
   21:  *
   22:  *      Complex circular tangent
   23:  *
   24:  *
   25:  *
   26:  * SYNOPSIS:
   27:  *
   28:  * double complex ctan();
   29:  * double complex z, w;
   30:  *
   31:  * w = ctan (z);
   32:  *
   33:  *
   34:  *
   35:  * DESCRIPTION:
   36:  *
   37:  * If
   38:  *     z = x + iy,
   39:  *
   40:  * then
   41:  *
   42:  *           sin 2x  +  i sinh 2y
   43:  *     w  =  --------------------.
   44:  *            cos 2x  +  cosh 2y
   45:  *
   46:  * On the real axis the denominator is zero at odd multiples
   47:  * of PI/2.  The denominator is evaluated by its Taylor
   48:  * series near these points.
   49:  *
   50:  * ctan(z) = -i ctanh(iz).
   51:  *
   52:  * ACCURACY:
   53:  *
   54:  *                      Relative error:
   55:  * arithmetic   domain     # trials      peak         rms
   56:  *    DEC       -10,+10      5200       7.1e-17     1.6e-17
   57:  *    IEEE      -10,+10     30000       7.2e-16     1.2e-16
   58:  * Also tested by ctan * ccot = 1 and catan(ctan(z))  =  z.
   59:  */
   60: 
   61: #include <sys/cdefs.h>
   62: #include <complex.h>
   63: #include <float.h>
   64: #include <math.h>
   65: 
   66: #define MACHEP 1.1e-16
   67: #define MAXNUM 1.0e308
   68: 
   69: static const double DP1 = 3.14159265160560607910E0;
   70: static const double DP2 = 1.98418714791870343106E-9;
   71: static const double DP3 = 1.14423774522196636802E-17;
   72: 
   73: static double
   74: _redupi(double x)
   75: {
   76:         double t;
   77:         long i;
   78: 
   79:         t = x/M_PI;
   80:         if (t >= 0.0)
   81:                 t += 0.5;
   82:         else
   83:                 t -= 0.5;
   84: 
   85:         i = t; /* the multiple */
   86:         t = i;
   87:         t = ((x - t * DP1) - t * DP2) - t * DP3;
   88:         return (t);
   89: }
   90: 
   91: /*  Taylor series expansion for cosh(2y) - cos(2x)      */
   92: 
   93: static double
   94: _ctans(double complex z)
   95: {
   96:         double f, x, x2, y, y2, rn, t;
   97:         double d;
   98: 
   99:         x = fabs (2.0 * creal (z));
  100:         y = fabs (2.0 * cimag(z));
  101: 
  102:         x = _redupi(x);
  103: 
  104:         x = x * x;
  105:         y = y * y;
  106:         x2 = 1.0;
  107:         y2 = 1.0;
  108:         f = 1.0;
  109:         rn = 0.0;
  110:         d = 0.0;
  111:         do {
  112:                 rn += 1.0;
  113:                 f *= rn;
  114:                 rn += 1.0;
  115:                 f *= rn;
  116:                 x2 *= x;
  117:                 y2 *= y;
  118:                 t = y2 + x2;
  119:                 t /= f;
  120:                 d += t;
  121: 
  122:                 rn += 1.0;
  123:                 f *= rn;
  124:                 rn += 1.0;
  125:                 f *= rn;
  126:                 x2 *= x;
  127:                 y2 *= y;
  128:                 t = y2 - x2;
  129:                 t /= f;
  130:                 d += t;
  131:         }
  132:         while (fabs(t/d) > MACHEP)
  133:                 ;
  134:         return (d);
  135: }
  136: 
  137: double complex
  138: ctan(double complex z)
  139: {
  140:         double complex w;
  141:         double d;
  142: 
  143:         d = cos (2.0 * creal (z)) + cosh (2.0 * cimag (z));
  144: 
  145:         if (fabs(d) < 0.25)
  146:                 d = _ctans (z);
  147: 
  148:         if (d == 0.0) {
  149:                 /*mtherr ("ctan", OVERFLOW);*/
  150:                 w = MAXNUM + MAXNUM * I;
  151:                 return (w);
  152:         }
  153: 
  154:         w = sin (2.0 * creal(z)) / d + (sinh (2.0 * cimag(z)) / d) * I;
  155:         return (w);
  156: }
  157: 
  158: #if     LDBL_MANT_DIG == 53
  159: #ifdef  lint
  160: /* PROTOLIB1 */
  161: long double complex ctanl(long double complex);
  162: #else   /* lint */
  163: __weak_alias(ctanl, ctan);
  164: #endif  /* lint */
  165: #endif  /* LDBL_MANT_DIG == 53 */