| 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 |