1 /* $OpenBSD: s_ctan.c,v 1.6 2013/07/03 04:46:36 espie 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 /* ctan()
19 *
20 * Complex circular tangent
21 *
22 *
23 *
24 * SYNOPSIS:
25 *
26 * double complex ctan();
27 * double complex z, w;
28 *
29 * w = ctan (z);
30 *
31 *
32 *
33 * DESCRIPTION:
34 *
35 * If
36 * z = x + iy,
37 *
38 * then
39 *
40 * sin 2x + i sinh 2y
41 * w = --------------------.
42 * cos 2x + cosh 2y
43 *
44 * On the real axis the denominator is zero at odd multiples
45 * of PI/2. The denominator is evaluated by its Taylor
46 * series near these points.
47 *
48 * ctan(z) = -i ctanh(iz).
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 #define MACHEP 1.1e-16
64 #define MAXNUM 1.0e308
65
66 static const double DP1 = 3.14159265160560607910E0;
67 static const double DP2 = 1.98418714791870343106E-9;
68 static const double DP3 = 1.14423774522196636802E-17;
69
70 static double
_redupi(double x)71 _redupi(double x)
72 {
73 double t;
74 long i;
75
76 t = x/M_PI;
77 if (t >= 0.0)
78 t += 0.5;
79 else
80 t -= 0.5;
81
82 i = t; /* the multiple */
83 t = i;
84 t = ((x - t * DP1) - t * DP2) - t * DP3;
85 return (t);
86 }
87
88 /* Taylor series expansion for cosh(2y) - cos(2x) */
89
90 static double
_ctans(double complex z)91 _ctans(double complex z)
92 {
93 double f, x, x2, y, y2, rn, t;
94 double d;
95
96 x = fabs (2.0 * creal (z));
97 y = fabs (2.0 * cimag(z));
98
99 x = _redupi(x);
100
101 x = x * x;
102 y = y * y;
103 x2 = 1.0;
104 y2 = 1.0;
105 f = 1.0;
106 rn = 0.0;
107 d = 0.0;
108 do {
109 rn += 1.0;
110 f *= rn;
111 rn += 1.0;
112 f *= rn;
113 x2 *= x;
114 y2 *= y;
115 t = y2 + x2;
116 t /= f;
117 d += t;
118
119 rn += 1.0;
120 f *= rn;
121 rn += 1.0;
122 f *= rn;
123 x2 *= x;
124 y2 *= y;
125 t = y2 - x2;
126 t /= f;
127 d += t;
128 }
129 while (fabs(t/d) > MACHEP)
130 ;
131 return (d);
132 }
133
134 double complex
ctan(double complex z)135 ctan(double complex z)
136 {
137 double complex w;
138 double d;
139
140 d = cos (2.0 * creal (z)) + cosh (2.0 * cimag (z));
141
142 if (fabs(d) < 0.25)
143 d = _ctans (z);
144
145 if (d == 0.0) {
146 /*mtherr ("ctan", OVERFLOW);*/
147 w = MAXNUM + MAXNUM * I;
148 return (w);
149 }
150
151 w = sin (2.0 * creal(z)) / d + (sinh (2.0 * cimag(z)) / d) * I;
152 return (w);
153 }
154
155 #if LDBL_MANT_DIG == DBL_MANT_DIG
156 __strong_alias(ctanl, ctan);
157 #endif /* LDBL_MANT_DIG == DBL_MANT_DIG */
158