File: | src/lib/libm/src/s_ctanl.c |
Warning: | line 146, column 3 Value stored to 'd' is never read |
Press '?' to see keyboard shortcuts
Keyboard shortcuts:
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_Complex.h> |
60 | #include <float.h> |
61 | #include <math.h> |
62 | |
63 | #if LDBL_MANT_DIG64 == 64 |
64 | static const long double MACHEPL= 5.42101086242752217003726400434970855712890625E-20L; |
65 | #elif LDBL_MANT_DIG64 == 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_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_Complex |
136 | ctanl(long double complex_Complex z) |
137 | { |
138 | long double complex_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); |
Value stored to 'd' is never read | |
147 | d = ctansl(z); |
148 | } |
149 | if (d == 0.0L) { |
150 | /*mtherr( "ctan", OVERFLOW );*/ |
151 | w = LDBL_MAX1.18973149535723176502e+4932L + LDBL_MAX1.18973149535723176502e+4932L * I1.0fi; |
152 | return (w); |
153 | } |
154 | |
155 | w = sinl(2.0L * x) / d + (sinhl(2.0L * y) / d) * I1.0fi; |
156 | return (w); |
157 | } |