OLD | NEW |
1 /* origin: OpenBSD /usr/src/lib/libm/src/s_catanf.c */ | 1 /* origin: OpenBSD /usr/src/lib/libm/src/s_catanf.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 43 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
54 */ | 54 */ |
55 | 55 |
56 #include "libm.h" | 56 #include "libm.h" |
57 | 57 |
58 #define MAXNUMF 1.0e38F | 58 #define MAXNUMF 1.0e38F |
59 | 59 |
60 static const double DP1 = 3.140625; | 60 static const double DP1 = 3.140625; |
61 static const double DP2 = 9.67502593994140625E-4; | 61 static const double DP2 = 9.67502593994140625E-4; |
62 static const double DP3 = 1.509957990978376432E-7; | 62 static const double DP3 = 1.509957990978376432E-7; |
63 | 63 |
64 static float _redupif(float xx) | 64 static float _redupif(float xx) { |
65 { | 65 float x, t; |
66 » float x, t; | 66 long i; |
67 » long i; | |
68 | 67 |
69 » x = xx; | 68 x = xx; |
70 » t = x/(float)M_PI; | 69 t = x / (float)M_PI; |
71 » if (t >= 0.0f) | 70 if (t >= 0.0f) |
72 » » t += 0.5f; | 71 t += 0.5f; |
73 » else | 72 else |
74 » » t -= 0.5f; | 73 t -= 0.5f; |
75 | 74 |
76 » i = t; /* the multiple */ | 75 i = t; /* the multiple */ |
77 » t = i; | 76 t = i; |
78 » t = ((x - t * DP1) - t * DP2) - t * DP3; | 77 t = ((x - t * DP1) - t * DP2) - t * DP3; |
79 » return t; | 78 return t; |
80 } | 79 } |
81 | 80 |
82 float complex catanf(float complex z) | 81 float complex catanf(float complex z) { |
83 { | 82 float complex w; |
84 » float complex w; | 83 float a, t, x, x2, y; |
85 » float a, t, x, x2, y; | |
86 | 84 |
87 » x = crealf(z); | 85 x = crealf(z); |
88 » y = cimagf(z); | 86 y = cimagf(z); |
89 | 87 |
90 » if ((x == 0.0f) && (y > 1.0f)) | 88 if ((x == 0.0f) && (y > 1.0f)) |
91 » » goto ovrf; | 89 goto ovrf; |
92 | 90 |
93 » x2 = x * x; | 91 x2 = x * x; |
94 » a = 1.0f - x2 - (y * y); | 92 a = 1.0f - x2 - (y * y); |
95 » if (a == 0.0f) | 93 if (a == 0.0f) |
96 » » goto ovrf; | 94 goto ovrf; |
97 | 95 |
98 » t = 0.5f * atan2f(2.0f * x, a); | 96 t = 0.5f * atan2f(2.0f * x, a); |
99 » w = _redupif(t); | 97 w = _redupif(t); |
100 | 98 |
101 » t = y - 1.0f; | 99 t = y - 1.0f; |
102 » a = x2 + (t * t); | 100 a = x2 + (t * t); |
103 » if (a == 0.0f) | 101 if (a == 0.0f) |
104 » » goto ovrf; | 102 goto ovrf; |
105 | 103 |
106 » t = y + 1.0f; | 104 t = y + 1.0f; |
107 » a = (x2 + (t * t))/a; | 105 a = (x2 + (t * t)) / a; |
108 » w = w + (0.25f * logf (a)) * I; | 106 w = w + (0.25f * logf(a)) * I; |
109 » return w; | 107 return w; |
110 | 108 |
111 ovrf: | 109 ovrf: |
112 » // FIXME | 110 // FIXME |
113 » w = MAXNUMF + MAXNUMF * I; | 111 w = MAXNUMF + MAXNUMF * I; |
114 » return w; | 112 return w; |
115 } | 113 } |
OLD | NEW |