1 /*	$OpenBSD: s_ctanf.c,v 1.2 2011/07/20 19:28:33 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 /*							ctanf()
19  *
20  *	Complex circular tangent
21  *
22  *
23  *
24  * SYNOPSIS:
25  *
26  * void ctanf();
27  * cmplxf z, w;
28  *
29  * ctanf( &z, &w );
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  *
49  * ACCURACY:
50  *
51  *                      Relative error:
52  * arithmetic   domain     # trials      peak         rms
53  *    IEEE      -10,+10     30000       3.3e-7       5.1e-8
54  */
55 
56 #include <complex.h>
57 #include <math.h>
58 
59 #define MACHEPF 3.0e-8
60 #define MAXNUMF 1.0e38f
61 
62 static const double DP1 = 3.140625;
63 static const double DP2 = 9.67502593994140625E-4;
64 static const double DP3 = 1.509957990978376432E-7;
65 
66 static float
67 _redupif(float xx)
68 {
69 	float x, t;
70 	long i;
71 
72 	x = xx;
73 	t = x/(float)M_PI;
74 	if(t >= 0.0)
75 		t += 0.5;
76 	else
77 		t -= 0.5;
78 
79 	i = t;	/* the multiple */
80 	t = i;
81 	t = ((x - t * DP1) - t * DP2) - t * DP3;
82 	return(t);
83 }
84 
85 /*  Taylor series expansion for cosh(2y) - cos(2x)	*/
86 
87 static float
88 _ctansf(float complex z)
89 {
90 	float f, x, x2, y, y2, rn, t, d;
91 
92 	x = fabsf(2.0f * crealf(z));
93 	y = fabsf(2.0f * cimagf(z));
94 
95 	x = _redupif(x);
96 
97 	x = x * x;
98 	y = y * y;
99 	x2 = 1.0f;
100 	y2 = 1.0f;
101 	f = 1.0f;
102 	rn = 0.0f;
103 	d = 0.0f;
104 	do {
105 		rn += 1.0f;
106 		f *= rn;
107 		rn += 1.0f;
108 		f *= rn;
109 		x2 *= x;
110 		y2 *= y;
111 		t = y2 + x2;
112 		t /= f;
113 		d += t;
114 
115 		rn += 1.0f;
116 		f *= rn;
117 		rn += 1.0f;
118 		f *= rn;
119 		x2 *= x;
120 		y2 *= y;
121 		t = y2 - x2;
122 		t /= f;
123 		d += t;
124 	}
125 	while (fabsf(t/d) > MACHEPF)
126 		;
127 	return(d);
128 }
129 
130 float complex
131 ctanf(float complex z)
132 {
133 	float complex w;
134 	float d;
135 
136 	d = cosf( 2.0f * crealf(z) ) + coshf( 2.0f * cimagf(z) );
137 
138 	if(fabsf(d) < 0.25f)
139 		d = _ctansf(z);
140 
141 	if (d == 0.0f) {
142 		/*mtherr( "ctanf", OVERFLOW );*/
143 		w = MAXNUMF + MAXNUMF * I;
144 		return (w);
145 	}
146 	w = sinf (2.0f * crealf(z)) / d + (sinhf (2.0f * cimagf(z)) / d) * I;
147 	return (w);
148 }
149