| OLD | NEW |
| 1 /* origin: OpenBSD /usr/src/lib/libm/src/s_catan.c */ | 1 /* origin: OpenBSD /usr/src/lib/libm/src/s_catan.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 48 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
| 59 */ | 59 */ |
| 60 | 60 |
| 61 #include "libm.h" | 61 #include "libm.h" |
| 62 | 62 |
| 63 #define MAXNUM 1.0e308 | 63 #define MAXNUM 1.0e308 |
| 64 | 64 |
| 65 static const double DP1 = 3.14159265160560607910E0; | 65 static const double DP1 = 3.14159265160560607910E0; |
| 66 static const double DP2 = 1.98418714791870343106E-9; | 66 static const double DP2 = 1.98418714791870343106E-9; |
| 67 static const double DP3 = 1.14423774522196636802E-17; | 67 static const double DP3 = 1.14423774522196636802E-17; |
| 68 | 68 |
| 69 static double _redupi(double x) | 69 static double _redupi(double x) { |
| 70 { | 70 double t; |
| 71 » double t; | 71 long i; |
| 72 » long i; | |
| 73 | 72 |
| 74 » t = x/M_PI; | 73 t = x / M_PI; |
| 75 » if (t >= 0.0) | 74 if (t >= 0.0) |
| 76 » » t += 0.5; | 75 t += 0.5; |
| 77 » else | 76 else |
| 78 » » t -= 0.5; | 77 t -= 0.5; |
| 79 | 78 |
| 80 » i = t; /* the multiple */ | 79 i = t; /* the multiple */ |
| 81 » t = i; | 80 t = i; |
| 82 » t = ((x - t * DP1) - t * DP2) - t * DP3; | 81 t = ((x - t * DP1) - t * DP2) - t * DP3; |
| 83 » return t; | 82 return t; |
| 84 } | 83 } |
| 85 | 84 |
| 86 double complex catan(double complex z) | 85 double complex catan(double complex z) { |
| 87 { | 86 double complex w; |
| 88 » double complex w; | 87 double a, t, x, x2, y; |
| 89 » double a, t, x, x2, y; | |
| 90 | 88 |
| 91 » x = creal(z); | 89 x = creal(z); |
| 92 » y = cimag(z); | 90 y = cimag(z); |
| 93 | 91 |
| 94 » if (x == 0.0 && y > 1.0) | 92 if (x == 0.0 && y > 1.0) |
| 95 » » goto ovrf; | 93 goto ovrf; |
| 96 | 94 |
| 97 » x2 = x * x; | 95 x2 = x * x; |
| 98 » a = 1.0 - x2 - (y * y); | 96 a = 1.0 - x2 - (y * y); |
| 99 » if (a == 0.0) | 97 if (a == 0.0) |
| 100 » » goto ovrf; | 98 goto ovrf; |
| 101 | 99 |
| 102 » t = 0.5 * atan2(2.0 * x, a); | 100 t = 0.5 * atan2(2.0 * x, a); |
| 103 » w = _redupi(t); | 101 w = _redupi(t); |
| 104 | 102 |
| 105 » t = y - 1.0; | 103 t = y - 1.0; |
| 106 » a = x2 + (t * t); | 104 a = x2 + (t * t); |
| 107 » if (a == 0.0) | 105 if (a == 0.0) |
| 108 » » goto ovrf; | 106 goto ovrf; |
| 109 | 107 |
| 110 » t = y + 1.0; | 108 t = y + 1.0; |
| 111 » a = (x2 + t * t)/a; | 109 a = (x2 + t * t) / a; |
| 112 » w = w + (0.25 * log(a)) * I; | 110 w = w + (0.25 * log(a)) * I; |
| 113 » return w; | 111 return w; |
| 114 | 112 |
| 115 ovrf: | 113 ovrf: |
| 116 » // FIXME | 114 // FIXME |
| 117 » w = MAXNUM + MAXNUM * I; | 115 w = MAXNUM + MAXNUM * I; |
| 118 » return w; | 116 return w; |
| 119 } | 117 } |
| OLD | NEW |