fix isnan
[xonotic/xonotic-data.pk3dir.git] / qcsrc / lib / warpzone / mathlib.qc
1 #include "mathlib.qh"
2 #if defined(CSQC)
3 #elif defined(MENUQC)
4 #elif defined(SVQC)
5 #endif
6
7 int fpclassify(float e)
8 {
9         if(isnan(e))
10                 return FP_NAN;
11         if(isinf(e))
12                 return FP_INFINITE;
13         if(e == 0)
14                 return FP_ZERO;
15         return FP_NORMAL;
16 }
17 bool isfinite(float e)
18 {
19         return !(isnan(e) || isinf(e));
20 }
21 bool isinf(float e)
22 {
23         return (e != 0) && (e + e == e);
24 }
25 bool isnan(float e)
26 {
27         // the sane way to detect NaN is broken because of a compiler bug
28         // (works with constants but breaks when assigned to variables)
29         // use conversion to string instead
30         //float f = e;
31         //return (e != f);
32
33         // Negative NaN ("-nan") is much more common but plain "nan" can be created by negating *some* -nans so we need to check both.
34         // DP's QCVM and GMQCC's QCVM behave differently - one needs ftos(-(0.0 / 0.0)), the other ftos(-sqrt(-1)).
35         string s = ftos(e);
36         return s == "nan" || s == "-nan";
37 }
38 bool isnormal(float e)
39 {
40         return isfinite(e);
41 }
42 bool signbit(float e)
43 {
44         return (e < 0);
45 }
46
47 float acosh(float e)
48 {
49         return log(e + sqrt(e*e - 1));
50 }
51 float asinh(float e)
52 {
53         return log(e + sqrt(e*e + 1));
54 }
55 float atanh(float e)
56 {
57         return 0.5 * log((1+e) / (1-e));
58 }
59 float cosh(float e)
60 {
61         return 0.5 * (exp(e) + exp(-e));
62 }
63 float sinh(float e)
64 {
65         return 0.5 * (exp(e) - exp(-e));
66 }
67 float tanh(float e)
68 {
69         return sinh(e) / cosh(e);
70 }
71
72 float exp(float e)
73 {
74         return pow(M_E, e);
75 }
76 float exp2(float e)
77 {
78         return pow(2, e);
79 }
80 float expm1(float e)
81 {
82         return exp(e) - 1;
83 }
84
85 vector frexp(float e)
86 {
87         vector v;
88         v.z = 0;
89         v.y = ilogb(e) + 1;
90         v.x = e / pow(2, v.y);
91         return v;
92 }
93 int ilogb(float e)
94 {
95         return floor(log2(fabs(e)));
96 }
97 float ldexp(float x, int e)
98 {
99         return x * pow(2, e);
100 }
101 float logn(float e, float base)
102 {
103         return log(e) / log(base);
104 }
105 float log10(float e)
106 {
107         return log(e) * M_LOG10E;
108 }
109 float log1p(float e)
110 {
111         return log(e + 1);
112 }
113 float log2(float e)
114 {
115         return log(e) * M_LOG2E;
116 }
117 float logb(float e)
118 {
119         return floor(log2(fabs(e)));
120 }
121 vector modf(float f)
122 {
123         return '1 0 0' * (f - trunc(f)) + '0 1 0' * trunc(f);
124 }
125
126 float scalbn(float e, int n)
127 {
128         return e * pow(2, n);
129 }
130
131 float cbrt(float e)
132 {
133         return copysign(pow(fabs(e), (1.0/3.0)), e);
134 }
135 float hypot(float e, float f)
136 {
137         return sqrt(e*e + f*f);
138 }
139
140 float erf(float e)
141 {
142         // approximation taken from wikipedia
143         float f;
144         f = e*e;
145         return copysign(sqrt(1 - exp(-f * (1.273239544735163 + 0.14001228868667 * f) / (1 + 0.14001228868667 * f))), e);
146 }
147 float erfc(float e)
148 {
149         return 1.0 - erf(e);
150 }
151 vector lgamma(float e)
152 {
153         // TODO improve accuracy
154         if(!isfinite(e))
155                 return fabs(e) * '1 0 0' + copysign(1, e) * '0 1 0';
156         if(e < 1 && e == floor(e))
157                 return nan("gamma") * '1 1 1';
158         if(e < 0.1)
159         {
160                 vector v;
161                 v = lgamma(1.0 - e);
162                 // reflection formula:
163                 // gamma(1-z) * gamma(z) = pi / sin(pi*z)
164                 // lgamma(1-z) + lgamma(z) = log(pi) - log(sin(pi*z))
165                 // sign of gamma(1-z) = sign of gamma(z) * sign of sin(pi*z)
166                 v.z = sin(M_PI * e);
167                 v.x = log(M_PI) - log(fabs(v.z)) - v.x;
168                 if(v.z < 0)
169                         v.y = -v.y;
170                 v.z = 0;
171                 return v;
172         }
173         if(e < 1.1)
174                 return lgamma(e + 1) - log(e) * '1 0 0';
175         e -= 1;
176         return (0.5 * log(2 * M_PI * e) + e * (log(e) - 1)) * '1 0 0' + '0 1 0';
177 }
178 float tgamma(float e)
179 {
180         vector v = lgamma(e);
181         return exp(v.x) * v.y;
182 }
183
184 /**
185  * Pythonic mod:
186  * TODO: %% operator?
187  *
188  *  1 %  2 ==  1
189  * -1 %  2 ==  1
190  *  1 % -2 == -1
191  * -1 % -2 == -1
192  */
193 float pymod(float e, float f)
194 {
195         return e - f * floor(e / f);
196 }
197
198 float nearbyint(float e)
199 {
200         return rint(e);
201 }
202 float trunc(float e)
203 {
204         return (e>=0) ? floor(e) : ceil(e);
205 }
206
207 float fmod(float e, float f)
208 {
209         return e - f * trunc(e / f);
210 }
211 float remainder(float e, float f)
212 {
213         return e - f * rint(e / f);
214 }
215 vector remquo(float e, float f)
216 {
217         vector v;
218         v.z = 0;
219         v.y = rint(e / f);
220         v.x = e - f * v.y;
221         return v;
222 }
223
224 float copysign(float e, float f)
225 {
226         return fabs(e) * ((f>0) ? 1 : -1);
227 }
228 /// Always use `isnan` function to compare because `float x = nan(); x == x;` gives true
229 float nan(string tag)
230 {
231         return sqrt(-1);
232 }
233 float nextafter(float e, float f)
234 {
235         // TODO very crude
236         if(e == f)
237                 return nan("nextafter");
238         if(e > f)
239                 return -nextafter(-e, -f);
240         // now we know that e < f
241         // so we need the next number > e
242         float d, a, b;
243         d = max(fabs(e), 0.00000000000000000000001);
244         a = e + d;
245         do
246         {
247                 d *= 0.5;
248                 b = a;
249                 a = e + d;
250         }
251         while(a != e);
252         return b;
253 }
254 float nexttoward(float e, float f)
255 {
256         return nextafter(e, f);
257 }
258
259 float fdim(float e, float f)
260 {
261         return max(e-f, 0);
262 }
263 float fmax(float e, float f)
264 {
265         return max(e, f);
266 }
267 float fmin(float e, float f)
268 {
269         return min(e, f);
270 }
271 float fma(float e, float f, float g)
272 {
273         return e * f + g;
274 }
275
276 int isgreater(float e, float f)
277 {
278         return e > f;
279 }
280 int isgreaterequal(float e, float f)
281 {
282         return e >= f;
283 }
284 int isless(float e, float f)
285 {
286         return e < f;
287 }
288 int islessequal(float e, float f)
289 {
290         return e <= f;
291 }
292 int islessgreater(float e, float f)
293 {
294         return e < f || e > f;
295 }
296 int isunordered(float e, float f)
297 {
298         return !(e < f || e == f || e > f);
299 }