Branch data Line data Source code
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 <float.h>
60 : : #include <openlibm_complex.h>
61 : : #include <openlibm_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 : 0 : redupil(long double x)
76 : : {
77 : : long double t;
78 : : long i;
79 : :
80 : 0 : t = x / PIL;
81 [ # # ]: 0 : if (t >= 0.0L)
82 : 0 : t += 0.5L;
83 : : else
84 : 0 : t -= 0.5L;
85 : :
86 : 0 : i = t; /* the multiple */
87 : 0 : t = i;
88 : 0 : t = ((x - t * DP1) - t * DP2) - t * DP3;
89 : 0 : return (t);
90 : : }
91 : :
92 : : static long double
93 : 0 : ctansl(long double complex z)
94 : : {
95 : : long double f, x, x2, y, y2, rn, t;
96 : : long double d;
97 : :
98 : 0 : x = fabsl(2.0L * creall(z));
99 : 0 : y = fabsl(2.0L * cimagl(z));
100 : :
101 : 0 : x = redupil(x);
102 : :
103 : 0 : x = x * x;
104 : 0 : y = y * y;
105 : 0 : x2 = 1.0L;
106 : 0 : y2 = 1.0L;
107 : 0 : f = 1.0L;
108 : 0 : rn = 0.0L;
109 : 0 : d = 0.0L;
110 : : do {
111 : 0 : rn += 1.0L;
112 : 0 : f *= rn;
113 : 0 : rn += 1.0L;
114 : 0 : f *= rn;
115 : 0 : x2 *= x;
116 : 0 : y2 *= y;
117 : 0 : t = y2 + x2;
118 : 0 : t /= f;
119 : 0 : d += t;
120 : :
121 : 0 : rn += 1.0L;
122 : 0 : f *= rn;
123 : 0 : rn += 1.0L;
124 : 0 : f *= rn;
125 : 0 : x2 *= x;
126 : 0 : y2 *= y;
127 : 0 : t = y2 - x2;
128 : 0 : t /= f;
129 : 0 : d += t;
130 : : }
131 [ # # ]: 0 : while (fabsl(t/d) > MACHEPL);
132 : 0 : return(d);
133 : : }
134 : :
135 : : long double complex
136 : 0 : ctanl(long double complex z)
137 : : {
138 : : long double complex w;
139 : : long double d, x, y;
140 : :
141 : 0 : x = creall(z);
142 : 0 : y = cimagl(z);
143 : 0 : d = cosl(2.0L * x) + coshl(2.0L * y);
144 : :
145 [ # # ]: 0 : if (fabsl(d) < 0.25L) {
146 : 0 : d = fabsl(d);
147 : 0 : d = ctansl(z);
148 : : }
149 [ # # ]: 0 : if (d == 0.0L) {
150 : : /*mtherr( "ctan", OVERFLOW );*/
151 : 0 : w = LDBL_MAX + LDBL_MAX * I;
152 : 0 : return (w);
153 : : }
154 : :
155 : 0 : w = sinl(2.0L * x) / d + (sinhl(2.0L * y) / d) * I;
156 : 0 : return (w);
157 : : }
|
