contrib/openbsd_libm/src/s_ctanl.c

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