65e53fb1ca7362e12272a6a6b8421b4dfcf82dde
[divverent/nexuiz.git] / data / qcsrc / common / mathlib.qc
1 /*
2 Copyright (c) 2009 Rudolf Polzer
3
4 Permission is hereby granted, free of charge, to any person obtaining a copy of
5 this software and associated documentation files (the "Software"), to deal in
6 the Software without restriction, including without limitation the rights to
7 use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
8 of the Software, and to permit persons to whom the Software is furnished to do
9 so, subject to the following conditions:
10
11 The above copyright notice and this permission notice shall be included in all
12 copies or substantial portions of the Software.
13
14 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
15 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
16 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
17 AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
18 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
19 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
20 SOFTWARE.
21 */
22
23 int fpclassify(float x)
24 {
25         if(isnan(x))
26                 return FP_NAN;
27         if(isinf(x))
28                 return FP_INFINITE;
29         if(x == 0)
30                 return FP_ZERO;
31         return FP_NORMAL;
32 }
33 int isfinite(float x)
34 {
35         return !(isnan(x) || isinf(x));
36 }
37 int isinf(float x)
38 {
39         return (x != 0) && (x + x == x);
40 }
41 int isnan(float x)
42 {
43         return isunordered(x, x);
44 }
45 int isnormal(float x)
46 {
47         return isfinite(x);
48 }
49 int signbit(float x)
50 {
51         return (x < 0);
52 }
53
54 float acosh(float x)
55 {
56         return log(x + sqrt(x*x - 1));
57 }
58 float asinh(float x)
59 {
60         return log(x + sqrt(x*x + 1));
61 }
62 float atanh(float x)
63 {
64         return 0.5 * log((1+x) / (1-x));
65 }
66 float cosh(float x)
67 {
68         return 0.5 * (exp(x) + exp(-x));
69 }
70 float sinh(float x)
71 {
72         return 0.5 * (exp(x) - exp(-x));
73 }
74 float tanh(float x)
75 {
76         return sinh(x) / cosh(x);
77 }
78
79 float exp(float x)
80 {
81         return pow(M_E, x);
82 }
83 float exp2(float x)
84 {
85         return pow(2, x);
86 }
87 float expm1(float x)
88 {
89         return exp(x) - 1;
90 }
91
92 vector frexp(float x)
93 {
94         vector v;
95         v_z = 0;
96         v_y = ilogb(x) + 1;
97         v_x = x / exp2(v_y);
98         return v;
99 }
100 int ilogb(float x)
101 {
102         return floor(log2(fabs(x)));
103 }
104 float ldexp(float x, int e)
105 {
106         return x * pow(2, e);
107 }
108 float log(float x)
109 {
110         // TODO improve speed
111         float i;
112         float r, r0;
113         if(x <= 0)
114                 return nan("log");
115         if(!isfinite(x))
116                 return x;
117         if(x >= 8)
118                 return -log(1 / x); // faster
119         if(x < 0.0009765625)
120                 return 2 * log(sqrt(x)); // faster
121         r = 1;
122         r0 = 0;
123         for(i = 1; fabs(r - r0) >= 0.00001; ++i)
124         {
125                 // Newton iteration on exp(r) = x:
126                 //   r <- r - (exp(r) - x) / (exp(r))
127                 //   r <- r - 1 + x / exp(r)
128                 r0 = r;
129                 r = r0 - 1 + x / exp(r0);
130         }
131         return r;
132 }
133 float log10(float x)
134 {
135         return log(x) * M_LOG10E;
136 }
137 float log1p(float x)
138 {
139         return log(x + 1);
140 }
141 float log2(float x)
142 {
143         return log(x) * M_LOG2E;
144 }
145 float logb(float x)
146 {
147         return floor(log2(fabs(x)));
148 }
149 vector modf(float f)
150 {
151         return '1 0 0' * (f - trunc(f)) + '0 1 0' * trunc(f);
152 }
153
154 float scalbn(float x, int n)
155 {
156         return x * pow(2, n);
157 }
158
159 float cbrt(float x)
160 {
161         return copysign(pow(fabs(x), 1.0/3.0), x);
162 }
163 float hypot(float x, float y)
164 {
165         return sqrt(x*x + y*y);
166 }
167
168 float erf(float x)
169 {
170         // approximation taken from wikipedia
171         float y;
172         y = x*x;
173         return copysign(sqrt(1 - exp(-y * (1.273239544735163 + 0.14001228868667 * y) / (1 + 0.14001228868667 * y))), x);
174 }
175 float erfc(float x)
176 {
177         return 1.0 - erf(x);
178 }
179 vector lgamma(float x)
180 {
181         // TODO improve accuracy
182         if(!isfinite(x))
183                 return fabs(x) * '1 0 0' + copysign(1, x) * '0 1 0';
184         if(x < 1 && x == floor(x))
185                 return nan("gamma") * '1 1 1';
186         if(x < 0.1)
187         {
188                 vector v;
189                 v = lgamma(1.0 - x);
190                 // reflection formula:
191                 // gamma(1-z) * gamma(z) = pi / sin(pi*z)
192                 // lgamma(1-z) + lgamma(z) = log(pi) - log(sin(pi*z))
193                 // sign of gamma(1-z) = sign of gamma(z) * sign of sin(pi*z)
194                 v_z = sin(M_PI * x);
195                 v_x = log(M_PI) - log(fabs(v_z)) - v_x;
196                 if(v_z < 0)
197                         v_y = -v_y;
198                 v_z = 0;
199                 return v;
200         }
201         if(x < 1.1)
202                 return lgamma(x + 1) - log(x) * '1 0 0';
203         x -= 1;
204         return (0.5 * log(2 * M_PI * x) + x * (log(x) - 1)) * '1 0 0' + '0 1 0';
205 }
206 float tgamma(float x)
207 {
208         vector v;
209         v = lgamma(x);
210         return exp(v_x) * v_y;
211 }
212
213 float nearbyint(float x)
214 {
215         return rint(x);
216 }
217 float trunc(float x)
218 {
219         return (x>=0) ? floor(x) : ceil(x);
220 }
221
222 float fmod(float x, float y)
223 {
224         return x - y * trunc(x / y);
225 }
226 float remainder(float x, float y)
227 {
228         return x - y * rint(x / y);
229 }
230 vector remquo(float x, float y)
231 {
232         vector v;
233         v_z = 0;
234         v_y = rint(x / y);
235         v_x = x - y * v_y;
236         return v;
237 }
238
239 float copysign(float x, float y)
240 {
241         return fabs(x) * ((y>0) ? 1 : -1);
242 }
243 float nan(string tag)
244 {
245         return sqrt(-1);
246 }
247 float nextafter(float x, float y)
248 {
249         // TODO very crude
250         if(x == y)
251                 return nan("nextafter");
252         if(x > y)
253                 return -nextafter(-x, -y);
254         // now we know that x < y
255         // so we need the next number > x
256         float d, a, b;
257         d = max(fabs(x), 0.00000000000000000000001);
258         a = x + d;
259         do
260         {
261                 d *= 0.5;
262                 b = a;
263                 a = x + d;
264         }
265         while(a != x);
266         return b;
267 }
268 float nexttoward(float x, float y)
269 {
270         return nextafter(x, y);
271 }
272
273 float fdim(float x, float y)
274 {
275         return max(x-y, 0);
276 }
277 float fmax(float x, float y)
278 {
279         return max(x, y);
280 }
281 float fmin(float x, float y)
282 {
283         return min(x, y);
284 }
285 float fma(float x, float y, float z)
286 {
287         return x * y + z;
288 }
289
290 int isgreater(float x, float y)
291 {
292         return x > y;
293 }
294 int isgreaterequal(float x, float y)
295 {
296         return x >= y;
297 }
298 int isless(float x, float y)
299 {
300         return x < y;
301 }
302 int islessequal(float x, float y)
303 {
304         return x <= y;
305 }
306 int islessgreater(float x, float y)
307 {
308         return x < y || x > y;
309 }
310 int isunordered(float x, float y)
311 {
312         return !(x < y || x == y || x > y);
313 }