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 |