]>
Commit | Line | Data |
---|---|---|
320054e8 DG |
1 | /* origin: OpenBSD /usr/src/lib/libm/src/s_catanl.c */ |
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 | * Complex circular arc tangent | |
19 | * | |
20 | * | |
21 | * SYNOPSIS: | |
22 | * | |
23 | * long double complex catanl(); | |
24 | * long double complex z, w; | |
25 | * | |
26 | * w = catanl( z ); | |
27 | * | |
28 | * | |
29 | * DESCRIPTION: | |
30 | * | |
31 | * If | |
32 | * z = x + iy, | |
33 | * | |
34 | * then | |
35 | * 1 ( 2x ) | |
36 | * Re w = - arctan(-----------) + k PI | |
37 | * 2 ( 2 2) | |
38 | * (1 - x - y ) | |
39 | * | |
40 | * ( 2 2) | |
41 | * 1 (x + (y+1) ) | |
42 | * Im w = - log(------------) | |
43 | * 4 ( 2 2) | |
44 | * (x + (y-1) ) | |
45 | * | |
46 | * Where k is an arbitrary integer. | |
47 | * | |
48 | * | |
49 | * ACCURACY: | |
50 | * | |
51 | * Relative error: | |
52 | * arithmetic domain # trials peak rms | |
53 | * DEC -10,+10 5900 1.3e-16 7.8e-18 | |
54 | * IEEE -10,+10 30000 2.3e-15 8.5e-17 | |
55 | * The check catan( ctan(z) ) = z, with |x| and |y| < PI/2, | |
56 | * had peak relative error 1.5e-16, rms relative error | |
57 | * 2.9e-17. See also clog(). | |
58 | */ | |
59 | ||
60 | #include <complex.h> | |
61 | #include <float.h> | |
d4db3fa2 | 62 | #include "complex_impl.h" |
320054e8 DG |
63 | |
64 | #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 | |
65 | long double complex catanl(long double complex z) | |
66 | { | |
67 | return catan(z); | |
68 | } | |
69 | #else | |
70 | static const long double PIL = 3.141592653589793238462643383279502884197169L; | |
71 | static const long double DP1 = 3.14159265358979323829596852490908531763125L; | |
72 | static const long double DP2 = 1.6667485837041756656403424829301998703007e-19L; | |
73 | static const long double DP3 = 1.8830410776607851167459095484560349402753e-39L; | |
74 | ||
75 | static long double 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 | long double complex catanl(long double complex z) | |
93 | { | |
94 | long double complex w; | |
95 | long double a, t, x, x2, y; | |
96 | ||
97 | x = creall(z); | |
98 | y = cimagl(z); | |
99 | ||
320054e8 DG |
100 | x2 = x * x; |
101 | a = 1.0L - x2 - (y * y); | |
320054e8 DG |
102 | |
103 | t = atan2l(2.0L * x, a) * 0.5L; | |
104 | w = redupil(t); | |
105 | ||
106 | t = y - 1.0L; | |
107 | a = x2 + (t * t); | |
320054e8 DG |
108 | |
109 | t = y + 1.0L; | |
110 | a = (x2 + (t * t)) / a; | |
575e1579 | 111 | w = CMPLXF(w, 0.25L * logl(a)); |
320054e8 DG |
112 | return w; |
113 | } | |
114 | #endif |