OLD | NEW |
1 /* origin: OpenBSD /usr/src/lib/libm/src/s_catanl.c */ | 1 /* origin: OpenBSD /usr/src/lib/libm/src/s_catanl.c */ |
2 /* | 2 /* |
3 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> | 3 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> |
4 * | 4 * |
5 * Permission to use, copy, modify, and distribute this software for any | 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 | 6 * purpose with or without fee is hereby granted, provided that the above |
7 * copyright notice and this permission notice appear in all copies. | 7 * copyright notice and this permission notice appear in all copies. |
8 * | 8 * |
9 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES | 9 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES |
10 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF | 10 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF |
(...skipping 44 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
55 * The check catan( ctan(z) ) = z, with |x| and |y| < PI/2, | 55 * The check catan( ctan(z) ) = z, with |x| and |y| < PI/2, |
56 * had peak relative error 1.5e-16, rms relative error | 56 * had peak relative error 1.5e-16, rms relative error |
57 * 2.9e-17. See also clog(). | 57 * 2.9e-17. See also clog(). |
58 */ | 58 */ |
59 | 59 |
60 #include <complex.h> | 60 #include <complex.h> |
61 #include <float.h> | 61 #include <float.h> |
62 #include "libm.h" | 62 #include "libm.h" |
63 | 63 |
64 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 | 64 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 |
65 long double complex catanl(long double complex z) | 65 long double complex catanl(long double complex z) { |
66 { | 66 return catan(z); |
67 » return catan(z); | |
68 } | 67 } |
69 #else | 68 #else |
70 static const long double PIL = 3.141592653589793238462643383279502884197169L; | 69 static const long double PIL = 3.141592653589793238462643383279502884197169L; |
71 static const long double DP1 = 3.14159265358979323829596852490908531763125L; | 70 static const long double DP1 = 3.14159265358979323829596852490908531763125L; |
72 static const long double DP2 = 1.6667485837041756656403424829301998703007e-19L; | 71 static const long double DP2 = 1.6667485837041756656403424829301998703007e-19L; |
73 static const long double DP3 = 1.8830410776607851167459095484560349402753e-39L; | 72 static const long double DP3 = 1.8830410776607851167459095484560349402753e-39L; |
74 | 73 |
75 static long double redupil(long double x) | 74 static long double redupil(long double x) { |
76 { | 75 long double t; |
77 » long double t; | 76 long i; |
78 » long i; | |
79 | 77 |
80 » t = x / PIL; | 78 t = x / PIL; |
81 » if (t >= 0.0L) | 79 if (t >= 0.0L) |
82 » » t += 0.5L; | 80 t += 0.5L; |
83 » else | 81 else |
84 » » t -= 0.5L; | 82 t -= 0.5L; |
85 | 83 |
86 » i = t; /* the multiple */ | 84 i = t; /* the multiple */ |
87 » t = i; | 85 t = i; |
88 » t = ((x - t * DP1) - t * DP2) - t * DP3; | 86 t = ((x - t * DP1) - t * DP2) - t * DP3; |
89 » return t; | 87 return t; |
90 } | 88 } |
91 | 89 |
92 long double complex catanl(long double complex z) | 90 long double complex catanl(long double complex z) { |
93 { | 91 long double complex w; |
94 » long double complex w; | 92 long double a, t, x, x2, y; |
95 » long double a, t, x, x2, y; | |
96 | 93 |
97 » x = creall(z); | 94 x = creall(z); |
98 » y = cimagl(z); | 95 y = cimagl(z); |
99 | 96 |
100 » if ((x == 0.0L) && (y > 1.0L)) | 97 if ((x == 0.0L) && (y > 1.0L)) |
101 » » goto ovrf; | 98 goto ovrf; |
102 | 99 |
103 » x2 = x * x; | 100 x2 = x * x; |
104 » a = 1.0L - x2 - (y * y); | 101 a = 1.0L - x2 - (y * y); |
105 » if (a == 0.0L) | 102 if (a == 0.0L) |
106 » » goto ovrf; | 103 goto ovrf; |
107 | 104 |
108 » t = atan2l(2.0L * x, a) * 0.5L; | 105 t = atan2l(2.0L * x, a) * 0.5L; |
109 » w = redupil(t); | 106 w = redupil(t); |
110 | 107 |
111 » t = y - 1.0L; | 108 t = y - 1.0L; |
112 » a = x2 + (t * t); | 109 a = x2 + (t * t); |
113 » if (a == 0.0L) | 110 if (a == 0.0L) |
114 » » goto ovrf; | 111 goto ovrf; |
115 | 112 |
116 » t = y + 1.0L; | 113 t = y + 1.0L; |
117 » a = (x2 + (t * t)) / a; | 114 a = (x2 + (t * t)) / a; |
118 » w = w + (0.25L * logl(a)) * I; | 115 w = w + (0.25L * logl(a)) * I; |
119 » return w; | 116 return w; |
120 | 117 |
121 ovrf: | 118 ovrf: |
122 » // FIXME | 119 // FIXME |
123 » w = LDBL_MAX + LDBL_MAX * I; | 120 w = LDBL_MAX + LDBL_MAX * I; |
124 » return w; | 121 return w; |
125 } | 122 } |
126 #endif | 123 #endif |
OLD | NEW |