optimize isnan()
[xonotic/xonotic-data.pk3dir.git] / qcsrc / warpzonelib / mathlib.qc
1 int fpclassify(float x)
2 {
3         if(isnan(x))
4                 return FP_NAN;
5         if(isinf(x))
6                 return FP_INFINITE;
7         if(x == 0)
8                 return FP_ZERO;
9         return FP_NORMAL;
10 }
11 int isfinite(float x)
12 {
13         return !(isnan(x) || isinf(x));
14 }
15 int isinf(float x)
16 {
17         return (x != 0) && (x + x == x);
18 }
19 int isnan(float x)
20 {
21         float y;
22         y = x;
23         return (x != y);
24 }
25 int isnormal(float x)
26 {
27         return isfinite(x);
28 }
29 int signbit(float x)
30 {
31         return (x < 0);
32 }
33
34 float acosh(float x)
35 {
36         return log(x + sqrt(x*x - 1));
37 }
38 float asinh(float x)
39 {
40         return log(x + sqrt(x*x + 1));
41 }
42 float atanh(float x)
43 {
44         return 0.5 * log((1+x) / (1-x));
45 }
46 float cosh(float x)
47 {
48         return 0.5 * (exp(x) + exp(-x));
49 }
50 float sinh(float x)
51 {
52         return 0.5 * (exp(x) - exp(-x));
53 }
54 float tanh(float x)
55 {
56         return sinh(x) / cosh(x);
57 }
58
59 float exp(float x)
60 {
61         return pow(M_E, x);
62 }
63 float exp2(float x)
64 {
65         return pow(2, x);
66 }
67 float expm1(float x)
68 {
69         return exp(x) - 1;
70 }
71
72 vector frexp(float x)
73 {
74         vector v;
75         v_z = 0;
76         v_y = ilogb(x) + 1;
77         v_x = x / exp2(v_y);
78         return v;
79 }
80 int ilogb(float x)
81 {
82         return floor(log2(fabs(x)));
83 }
84 float ldexp(float x, int e)
85 {
86         return x * pow(2, e);
87 }
88 float log10(float x)
89 {
90         return log(x) * M_LOG10E;
91 }
92 float log1p(float x)
93 {
94         return log(x + 1);
95 }
96 float log2(float x)
97 {
98         return log(x) * M_LOG2E;
99 }
100 float logb(float x)
101 {
102         return floor(log2(fabs(x)));
103 }
104 vector modf(float f)
105 {
106         return '1 0 0' * (f - trunc(f)) + '0 1 0' * trunc(f);
107 }
108
109 float scalbn(float x, int n)
110 {
111         return x * pow(2, n);
112 }
113
114 float cbrt(float x)
115 {
116         return copysign(pow(fabs(x), 1.0/3.0), x);
117 }
118 float hypot(float x, float y)
119 {
120         return sqrt(x*x + y*y);
121 }
122
123 float erf(float x)
124 {
125         // approximation taken from wikipedia
126         float y;
127         y = x*x;
128         return copysign(sqrt(1 - exp(-y * (1.273239544735163 + 0.14001228868667 * y) / (1 + 0.14001228868667 * y))), x);
129 }
130 float erfc(float x)
131 {
132         return 1.0 - erf(x);
133 }
134 vector lgamma(float x)
135 {
136         // TODO improve accuracy
137         if(!isfinite(x))
138                 return fabs(x) * '1 0 0' + copysign(1, x) * '0 1 0';
139         if(x < 1 && x == floor(x))
140                 return nan("gamma") * '1 1 1';
141         if(x < 0.1)
142         {
143                 vector v;
144                 v = lgamma(1.0 - x);
145                 // reflection formula:
146                 // gamma(1-z) * gamma(z) = pi / sin(pi*z)
147                 // lgamma(1-z) + lgamma(z) = log(pi) - log(sin(pi*z))
148                 // sign of gamma(1-z) = sign of gamma(z) * sign of sin(pi*z)
149                 v_z = sin(M_PI * x);
150                 v_x = log(M_PI) - log(fabs(v_z)) - v_x;
151                 if(v_z < 0)
152                         v_y = -v_y;
153                 v_z = 0;
154                 return v;
155         }
156         if(x < 1.1)
157                 return lgamma(x + 1) - log(x) * '1 0 0';
158         x -= 1;
159         return (0.5 * log(2 * M_PI * x) + x * (log(x) - 1)) * '1 0 0' + '0 1 0';
160 }
161 float tgamma(float x)
162 {
163         vector v;
164         v = lgamma(x);
165         return exp(v_x) * v_y;
166 }
167
168 float nearbyint(float x)
169 {
170         return rint(x);
171 }
172 float trunc(float x)
173 {
174         return (x>=0) ? floor(x) : ceil(x);
175 }
176
177 float fmod(float x, float y)
178 {
179         return x - y * trunc(x / y);
180 }
181 float remainder(float x, float y)
182 {
183         return x - y * rint(x / y);
184 }
185 vector remquo(float x, float y)
186 {
187         vector v;
188         v_z = 0;
189         v_y = rint(x / y);
190         v_x = x - y * v_y;
191         return v;
192 }
193
194 float copysign(float x, float y)
195 {
196         return fabs(x) * ((y>0) ? 1 : -1);
197 }
198 float nan(string tag)
199 {
200         return sqrt(-1);
201 }
202 float nextafter(float x, float y)
203 {
204         // TODO very crude
205         if(x == y)
206                 return nan("nextafter");
207         if(x > y)
208                 return -nextafter(-x, -y);
209         // now we know that x < y
210         // so we need the next number > x
211         float d, a, b;
212         d = max(fabs(x), 0.00000000000000000000001);
213         a = x + d;
214         do
215         {
216                 d *= 0.5;
217                 b = a;
218                 a = x + d;
219         }
220         while(a != x);
221         return b;
222 }
223 float nexttoward(float x, float y)
224 {
225         return nextafter(x, y);
226 }
227
228 float fdim(float x, float y)
229 {
230         return max(x-y, 0);
231 }
232 float fmax(float x, float y)
233 {
234         return max(x, y);
235 }
236 float fmin(float x, float y)
237 {
238         return min(x, y);
239 }
240 float fma(float x, float y, float z)
241 {
242         return x * y + z;
243 }
244
245 int isgreater(float x, float y)
246 {
247         return x > y;
248 }
249 int isgreaterequal(float x, float y)
250 {
251         return x >= y;
252 }
253 int isless(float x, float y)
254 {
255         return x < y;
256 }
257 int islessequal(float x, float y)
258 {
259         return x <= y;
260 }
261 int islessgreater(float x, float y)
262 {
263         return x < y || x > y;
264 }
265 int isunordered(float x, float y)
266 {
267         return !(x < y || x == y || x > y);
268 }