1/*	$OpenBSD: s_ctan.c,v 1.7 2016/09/12 19:47:02 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/*							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
66static const double DP1 = 3.14159265160560607910E0;
67static const double DP2 = 1.98418714791870343106E-9;
68static const double DP3 = 1.14423774522196636802E-17;
69
70static double
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
90static double
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
134double complex
135ctan(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}
154DEF_STD(ctan);
155LDBL_MAYBE_UNUSED_CLONE(ctan);
156