Fix setinfo.
[xonotic/darkplaces.git] / mathlib.c
1 /*
2 Copyright (C) 1996-1997 Id Software, Inc.
3
4 This program is free software; you can redistribute it and/or
5 modify it under the terms of the GNU General Public License
6 as published by the Free Software Foundation; either version 2
7 of the License, or (at your option) any later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
12
13 See the GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public License
16 along with this program; if not, write to the Free Software
17 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
18
19 */
20 // mathlib.c -- math primitives
21
22 #include "quakedef.h"
23
24 #include <math.h>
25
26 vec3_t vec3_origin = {0,0,0};
27 float ixtable[4096];
28
29 /*-----------------------------------------------------------------*/
30
31 float m_bytenormals[NUMVERTEXNORMALS][3] =
32 {
33 {-0.525731f, 0.000000f, 0.850651f}, {-0.442863f, 0.238856f, 0.864188f},
34 {-0.295242f, 0.000000f, 0.955423f}, {-0.309017f, 0.500000f, 0.809017f},
35 {-0.162460f, 0.262866f, 0.951056f}, {0.000000f, 0.000000f, 1.000000f},
36 {0.000000f, 0.850651f, 0.525731f}, {-0.147621f, 0.716567f, 0.681718f},
37 {0.147621f, 0.716567f, 0.681718f}, {0.000000f, 0.525731f, 0.850651f},
38 {0.309017f, 0.500000f, 0.809017f}, {0.525731f, 0.000000f, 0.850651f},
39 {0.295242f, 0.000000f, 0.955423f}, {0.442863f, 0.238856f, 0.864188f},
40 {0.162460f, 0.262866f, 0.951056f}, {-0.681718f, 0.147621f, 0.716567f},
41 {-0.809017f, 0.309017f, 0.500000f}, {-0.587785f, 0.425325f, 0.688191f},
42 {-0.850651f, 0.525731f, 0.000000f}, {-0.864188f, 0.442863f, 0.238856f},
43 {-0.716567f, 0.681718f, 0.147621f}, {-0.688191f, 0.587785f, 0.425325f},
44 {-0.500000f, 0.809017f, 0.309017f}, {-0.238856f, 0.864188f, 0.442863f},
45 {-0.425325f, 0.688191f, 0.587785f}, {-0.716567f, 0.681718f, -0.147621f},
46 {-0.500000f, 0.809017f, -0.309017f}, {-0.525731f, 0.850651f, 0.000000f},
47 {0.000000f, 0.850651f, -0.525731f}, {-0.238856f, 0.864188f, -0.442863f},
48 {0.000000f, 0.955423f, -0.295242f}, {-0.262866f, 0.951056f, -0.162460f},
49 {0.000000f, 1.000000f, 0.000000f}, {0.000000f, 0.955423f, 0.295242f},
50 {-0.262866f, 0.951056f, 0.162460f}, {0.238856f, 0.864188f, 0.442863f},
51 {0.262866f, 0.951056f, 0.162460f}, {0.500000f, 0.809017f, 0.309017f},
52 {0.238856f, 0.864188f, -0.442863f}, {0.262866f, 0.951056f, -0.162460f},
53 {0.500000f, 0.809017f, -0.309017f}, {0.850651f, 0.525731f, 0.000000f},
54 {0.716567f, 0.681718f, 0.147621f}, {0.716567f, 0.681718f, -0.147621f},
55 {0.525731f, 0.850651f, 0.000000f}, {0.425325f, 0.688191f, 0.587785f},
56 {0.864188f, 0.442863f, 0.238856f}, {0.688191f, 0.587785f, 0.425325f},
57 {0.809017f, 0.309017f, 0.500000f}, {0.681718f, 0.147621f, 0.716567f},
58 {0.587785f, 0.425325f, 0.688191f}, {0.955423f, 0.295242f, 0.000000f},
59 {1.000000f, 0.000000f, 0.000000f}, {0.951056f, 0.162460f, 0.262866f},
60 {0.850651f, -0.525731f, 0.000000f}, {0.955423f, -0.295242f, 0.000000f},
61 {0.864188f, -0.442863f, 0.238856f}, {0.951056f, -0.162460f, 0.262866f},
62 {0.809017f, -0.309017f, 0.500000f}, {0.681718f, -0.147621f, 0.716567f},
63 {0.850651f, 0.000000f, 0.525731f}, {0.864188f, 0.442863f, -0.238856f},
64 {0.809017f, 0.309017f, -0.500000f}, {0.951056f, 0.162460f, -0.262866f},
65 {0.525731f, 0.000000f, -0.850651f}, {0.681718f, 0.147621f, -0.716567f},
66 {0.681718f, -0.147621f, -0.716567f}, {0.850651f, 0.000000f, -0.525731f},
67 {0.809017f, -0.309017f, -0.500000f}, {0.864188f, -0.442863f, -0.238856f},
68 {0.951056f, -0.162460f, -0.262866f}, {0.147621f, 0.716567f, -0.681718f},
69 {0.309017f, 0.500000f, -0.809017f}, {0.425325f, 0.688191f, -0.587785f},
70 {0.442863f, 0.238856f, -0.864188f}, {0.587785f, 0.425325f, -0.688191f},
71 {0.688191f, 0.587785f, -0.425325f}, {-0.147621f, 0.716567f, -0.681718f},
72 {-0.309017f, 0.500000f, -0.809017f}, {0.000000f, 0.525731f, -0.850651f},
73 {-0.525731f, 0.000000f, -0.850651f}, {-0.442863f, 0.238856f, -0.864188f},
74 {-0.295242f, 0.000000f, -0.955423f}, {-0.162460f, 0.262866f, -0.951056f},
75 {0.000000f, 0.000000f, -1.000000f}, {0.295242f, 0.000000f, -0.955423f},
76 {0.162460f, 0.262866f, -0.951056f}, {-0.442863f, -0.238856f, -0.864188f},
77 {-0.309017f, -0.500000f, -0.809017f}, {-0.162460f, -0.262866f, -0.951056f},
78 {0.000000f, -0.850651f, -0.525731f}, {-0.147621f, -0.716567f, -0.681718f},
79 {0.147621f, -0.716567f, -0.681718f}, {0.000000f, -0.525731f, -0.850651f},
80 {0.309017f, -0.500000f, -0.809017f}, {0.442863f, -0.238856f, -0.864188f},
81 {0.162460f, -0.262866f, -0.951056f}, {0.238856f, -0.864188f, -0.442863f},
82 {0.500000f, -0.809017f, -0.309017f}, {0.425325f, -0.688191f, -0.587785f},
83 {0.716567f, -0.681718f, -0.147621f}, {0.688191f, -0.587785f, -0.425325f},
84 {0.587785f, -0.425325f, -0.688191f}, {0.000000f, -0.955423f, -0.295242f},
85 {0.000000f, -1.000000f, 0.000000f}, {0.262866f, -0.951056f, -0.162460f},
86 {0.000000f, -0.850651f, 0.525731f}, {0.000000f, -0.955423f, 0.295242f},
87 {0.238856f, -0.864188f, 0.442863f}, {0.262866f, -0.951056f, 0.162460f},
88 {0.500000f, -0.809017f, 0.309017f}, {0.716567f, -0.681718f, 0.147621f},
89 {0.525731f, -0.850651f, 0.000000f}, {-0.238856f, -0.864188f, -0.442863f},
90 {-0.500000f, -0.809017f, -0.309017f}, {-0.262866f, -0.951056f, -0.162460f},
91 {-0.850651f, -0.525731f, 0.000000f}, {-0.716567f, -0.681718f, -0.147621f},
92 {-0.716567f, -0.681718f, 0.147621f}, {-0.525731f, -0.850651f, 0.000000f},
93 {-0.500000f, -0.809017f, 0.309017f}, {-0.238856f, -0.864188f, 0.442863f},
94 {-0.262866f, -0.951056f, 0.162460f}, {-0.864188f, -0.442863f, 0.238856f},
95 {-0.809017f, -0.309017f, 0.500000f}, {-0.688191f, -0.587785f, 0.425325f},
96 {-0.681718f, -0.147621f, 0.716567f}, {-0.442863f, -0.238856f, 0.864188f},
97 {-0.587785f, -0.425325f, 0.688191f}, {-0.309017f, -0.500000f, 0.809017f},
98 {-0.147621f, -0.716567f, 0.681718f}, {-0.425325f, -0.688191f, 0.587785f},
99 {-0.162460f, -0.262866f, 0.951056f}, {0.442863f, -0.238856f, 0.864188f},
100 {0.162460f, -0.262866f, 0.951056f}, {0.309017f, -0.500000f, 0.809017f},
101 {0.147621f, -0.716567f, 0.681718f}, {0.000000f, -0.525731f, 0.850651f},
102 {0.425325f, -0.688191f, 0.587785f}, {0.587785f, -0.425325f, 0.688191f},
103 {0.688191f, -0.587785f, 0.425325f}, {-0.955423f, 0.295242f, 0.000000f},
104 {-0.951056f, 0.162460f, 0.262866f}, {-1.000000f, 0.000000f, 0.000000f},
105 {-0.850651f, 0.000000f, 0.525731f}, {-0.955423f, -0.295242f, 0.000000f},
106 {-0.951056f, -0.162460f, 0.262866f}, {-0.864188f, 0.442863f, -0.238856f},
107 {-0.951056f, 0.162460f, -0.262866f}, {-0.809017f, 0.309017f, -0.500000f},
108 {-0.864188f, -0.442863f, -0.238856f}, {-0.951056f, -0.162460f, -0.262866f},
109 {-0.809017f, -0.309017f, -0.500000f}, {-0.681718f, 0.147621f, -0.716567f},
110 {-0.681718f, -0.147621f, -0.716567f}, {-0.850651f, 0.000000f, -0.525731f},
111 {-0.688191f, 0.587785f, -0.425325f}, {-0.587785f, 0.425325f, -0.688191f},
112 {-0.425325f, 0.688191f, -0.587785f}, {-0.425325f, -0.688191f, -0.587785f},
113 {-0.587785f, -0.425325f, -0.688191f}, {-0.688191f, -0.587785f, -0.425325f},
114 };
115
116 #if 0
117 unsigned char NormalToByte(const vec3_t n)
118 {
119         int i, best;
120         float bestdistance, distance;
121
122         best = 0;
123         bestdistance = DotProduct (n, m_bytenormals[0]);
124         for (i = 1;i < NUMVERTEXNORMALS;i++)
125         {
126                 distance = DotProduct (n, m_bytenormals[i]);
127                 if (distance > bestdistance)
128                 {
129                         bestdistance = distance;
130                         best = i;
131                 }
132         }
133         return best;
134 }
135
136 // note: uses byte partly to force unsigned for the validity check
137 void ByteToNormal(unsigned char num, vec3_t n)
138 {
139         if (num < NUMVERTEXNORMALS)
140                 VectorCopy(m_bytenormals[num], n);
141         else
142                 VectorClear(n); // FIXME: complain?
143 }
144
145 // assumes "src" is normalized
146 void PerpendicularVector( vec3_t dst, const vec3_t src )
147 {
148         // LordHavoc: optimized to death and beyond
149         int pos;
150         float minelem;
151
152         if (src[0])
153         {
154                 dst[0] = 0;
155                 if (src[1])
156                 {
157                         dst[1] = 0;
158                         if (src[2])
159                         {
160                                 dst[2] = 0;
161                                 pos = 0;
162                                 minelem = fabs(src[0]);
163                                 if (fabs(src[1]) < minelem)
164                                 {
165                                         pos = 1;
166                                         minelem = fabs(src[1]);
167                                 }
168                                 if (fabs(src[2]) < minelem)
169                                         pos = 2;
170
171                                 dst[pos] = 1;
172                                 dst[0] -= src[pos] * src[0];
173                                 dst[1] -= src[pos] * src[1];
174                                 dst[2] -= src[pos] * src[2];
175
176                                 // normalize the result
177                                 VectorNormalize(dst);
178                         }
179                         else
180                                 dst[2] = 1;
181                 }
182                 else
183                 {
184                         dst[1] = 1;
185                         dst[2] = 0;
186                 }
187         }
188         else
189         {
190                 dst[0] = 1;
191                 dst[1] = 0;
192                 dst[2] = 0;
193         }
194 }
195 #endif
196
197
198 // LordHavoc: like AngleVectors, but taking a forward vector instead of angles, useful!
199 void VectorVectors(const vec3_t forward, vec3_t right, vec3_t up)
200 {
201         // NOTE: this is consistent to AngleVectors applied to AnglesFromVectors
202         if (forward[0] == 0 && forward[1] == 0)
203         {
204                 if(forward[2] > 0)
205                 {
206                         VectorSet(right, 0, -1, 0);
207                         VectorSet(up, -1, 0, 0);
208                 }
209                 else
210                 {
211                         VectorSet(right, 0, -1, 0);
212                         VectorSet(up, 1, 0, 0);
213                 }
214         }
215         else
216         {
217                 right[0] = forward[1];
218                 right[1] = -forward[0];
219                 right[2] = 0;
220                 VectorNormalize(right);
221
222                 up[0] = (-forward[2]*forward[0]);
223                 up[1] = (-forward[2]*forward[1]);
224                 up[2] = (forward[0]*forward[0] + forward[1]*forward[1]);
225                 VectorNormalize(up);
226         }
227 }
228
229 void VectorVectorsDouble(const double *forward, double *right, double *up)
230 {
231         if (forward[0] == 0 && forward[1] == 0)
232         {
233                 if(forward[2] > 0)
234                 {
235                         VectorSet(right, 0, -1, 0);
236                         VectorSet(up, -1, 0, 0);
237                 }
238                 else
239                 {
240                         VectorSet(right, 0, -1, 0);
241                         VectorSet(up, 1, 0, 0);
242                 }
243         }
244         else
245         {
246                 right[0] = forward[1];
247                 right[1] = -forward[0];
248                 right[2] = 0;
249                 VectorNormalize(right);
250
251                 up[0] = (-forward[2]*forward[0]);
252                 up[1] = (-forward[2]*forward[1]);
253                 up[2] = (forward[0]*forward[0] + forward[1]*forward[1]);
254                 VectorNormalize(up);
255         }
256 }
257
258 void RotatePointAroundVector( vec3_t dst, const vec3_t dir, const vec3_t point, float degrees )
259 {
260         float t0, t1;
261         float angle, c, s;
262         vec3_t vr, vu, vf;
263
264         angle = DEG2RAD(degrees);
265         c = cos(angle);
266         s = sin(angle);
267         VectorCopy(dir, vf);
268         VectorVectors(vf, vr, vu);
269
270         t0 = vr[0] *  c + vu[0] * -s;
271         t1 = vr[0] *  s + vu[0] *  c;
272         dst[0] = (t0 * vr[0] + t1 * vu[0] + vf[0] * vf[0]) * point[0]
273                + (t0 * vr[1] + t1 * vu[1] + vf[0] * vf[1]) * point[1]
274                + (t0 * vr[2] + t1 * vu[2] + vf[0] * vf[2]) * point[2];
275
276         t0 = vr[1] *  c + vu[1] * -s;
277         t1 = vr[1] *  s + vu[1] *  c;
278         dst[1] = (t0 * vr[0] + t1 * vu[0] + vf[1] * vf[0]) * point[0]
279                + (t0 * vr[1] + t1 * vu[1] + vf[1] * vf[1]) * point[1]
280                + (t0 * vr[2] + t1 * vu[2] + vf[1] * vf[2]) * point[2];
281
282         t0 = vr[2] *  c + vu[2] * -s;
283         t1 = vr[2] *  s + vu[2] *  c;
284         dst[2] = (t0 * vr[0] + t1 * vu[0] + vf[2] * vf[0]) * point[0]
285                + (t0 * vr[1] + t1 * vu[1] + vf[2] * vf[1]) * point[1]
286                + (t0 * vr[2] + t1 * vu[2] + vf[2] * vf[2]) * point[2];
287 }
288
289 /*-----------------------------------------------------------------*/
290
291 // returns the smallest integer greater than or equal to "value", or 0 if "value" is too big
292 unsigned int CeilPowerOf2(unsigned int value)
293 {
294         unsigned int ceilvalue;
295
296         if (value > (1U << (sizeof(int) * 8 - 1)))
297                 return 0;
298
299         ceilvalue = 1;
300         while (ceilvalue < value)
301                 ceilvalue <<= 1;
302
303         return ceilvalue;
304 }
305
306
307 /*-----------------------------------------------------------------*/
308
309
310 void PlaneClassify(mplane_t *p)
311 {
312         // for optimized plane comparisons
313         if (p->normal[0] == 1)
314                 p->type = 0;
315         else if (p->normal[1] == 1)
316                 p->type = 1;
317         else if (p->normal[2] == 1)
318                 p->type = 2;
319         else
320                 p->type = 3;
321         // for BoxOnPlaneSide
322         p->signbits = 0;
323         if (p->normal[0] < 0) // 1
324                 p->signbits |= 1;
325         if (p->normal[1] < 0) // 2
326                 p->signbits |= 2;
327         if (p->normal[2] < 0) // 4
328                 p->signbits |= 4;
329 }
330
331 int BoxOnPlaneSide(const vec3_t emins, const vec3_t emaxs, const mplane_t *p)
332 {
333         if (p->type < 3)
334                 return ((emaxs[p->type] >= p->dist) | ((emins[p->type] < p->dist) << 1));
335         switch(p->signbits)
336         {
337         default:
338         case 0: return (((p->normal[0] * emaxs[0] + p->normal[1] * emaxs[1] + p->normal[2] * emaxs[2]) >= p->dist) | (((p->normal[0] * emins[0] + p->normal[1] * emins[1] + p->normal[2] * emins[2]) < p->dist) << 1));
339         case 1: return (((p->normal[0] * emins[0] + p->normal[1] * emaxs[1] + p->normal[2] * emaxs[2]) >= p->dist) | (((p->normal[0] * emaxs[0] + p->normal[1] * emins[1] + p->normal[2] * emins[2]) < p->dist) << 1));
340         case 2: return (((p->normal[0] * emaxs[0] + p->normal[1] * emins[1] + p->normal[2] * emaxs[2]) >= p->dist) | (((p->normal[0] * emins[0] + p->normal[1] * emaxs[1] + p->normal[2] * emins[2]) < p->dist) << 1));
341         case 3: return (((p->normal[0] * emins[0] + p->normal[1] * emins[1] + p->normal[2] * emaxs[2]) >= p->dist) | (((p->normal[0] * emaxs[0] + p->normal[1] * emaxs[1] + p->normal[2] * emins[2]) < p->dist) << 1));
342         case 4: return (((p->normal[0] * emaxs[0] + p->normal[1] * emaxs[1] + p->normal[2] * emins[2]) >= p->dist) | (((p->normal[0] * emins[0] + p->normal[1] * emins[1] + p->normal[2] * emaxs[2]) < p->dist) << 1));
343         case 5: return (((p->normal[0] * emins[0] + p->normal[1] * emaxs[1] + p->normal[2] * emins[2]) >= p->dist) | (((p->normal[0] * emaxs[0] + p->normal[1] * emins[1] + p->normal[2] * emaxs[2]) < p->dist) << 1));
344         case 6: return (((p->normal[0] * emaxs[0] + p->normal[1] * emins[1] + p->normal[2] * emins[2]) >= p->dist) | (((p->normal[0] * emins[0] + p->normal[1] * emaxs[1] + p->normal[2] * emaxs[2]) < p->dist) << 1));
345         case 7: return (((p->normal[0] * emins[0] + p->normal[1] * emins[1] + p->normal[2] * emins[2]) >= p->dist) | (((p->normal[0] * emaxs[0] + p->normal[1] * emaxs[1] + p->normal[2] * emaxs[2]) < p->dist) << 1));
346         }
347 }
348
349 #if 0
350 int BoxOnPlaneSide_Separate(const vec3_t emins, const vec3_t emaxs, const vec3_t normal, const vec_t dist)
351 {
352         switch((normal[0] < 0) | ((normal[1] < 0) << 1) | ((normal[2] < 0) << 2))
353         {
354         default:
355         case 0: return (((normal[0] * emaxs[0] + normal[1] * emaxs[1] + normal[2] * emaxs[2]) >= dist) | (((normal[0] * emins[0] + normal[1] * emins[1] + normal[2] * emins[2]) < dist) << 1));
356         case 1: return (((normal[0] * emins[0] + normal[1] * emaxs[1] + normal[2] * emaxs[2]) >= dist) | (((normal[0] * emaxs[0] + normal[1] * emins[1] + normal[2] * emins[2]) < dist) << 1));
357         case 2: return (((normal[0] * emaxs[0] + normal[1] * emins[1] + normal[2] * emaxs[2]) >= dist) | (((normal[0] * emins[0] + normal[1] * emaxs[1] + normal[2] * emins[2]) < dist) << 1));
358         case 3: return (((normal[0] * emins[0] + normal[1] * emins[1] + normal[2] * emaxs[2]) >= dist) | (((normal[0] * emaxs[0] + normal[1] * emaxs[1] + normal[2] * emins[2]) < dist) << 1));
359         case 4: return (((normal[0] * emaxs[0] + normal[1] * emaxs[1] + normal[2] * emins[2]) >= dist) | (((normal[0] * emins[0] + normal[1] * emins[1] + normal[2] * emaxs[2]) < dist) << 1));
360         case 5: return (((normal[0] * emins[0] + normal[1] * emaxs[1] + normal[2] * emins[2]) >= dist) | (((normal[0] * emaxs[0] + normal[1] * emins[1] + normal[2] * emaxs[2]) < dist) << 1));
361         case 6: return (((normal[0] * emaxs[0] + normal[1] * emins[1] + normal[2] * emins[2]) >= dist) | (((normal[0] * emins[0] + normal[1] * emaxs[1] + normal[2] * emaxs[2]) < dist) << 1));
362         case 7: return (((normal[0] * emins[0] + normal[1] * emins[1] + normal[2] * emins[2]) >= dist) | (((normal[0] * emaxs[0] + normal[1] * emaxs[1] + normal[2] * emaxs[2]) < dist) << 1));
363         }
364 }
365 #endif
366
367 void BoxPlaneCorners(const vec3_t emins, const vec3_t emaxs, const mplane_t *p, vec3_t outnear, vec3_t outfar)
368 {
369         if (p->type < 3)
370         {
371                 outnear[0] = outnear[1] = outnear[2] = outfar[0] = outfar[1] = outfar[2] = 0;
372                 outnear[p->type] = emins[p->type];
373                 outfar[p->type] = emaxs[p->type];
374                 return;
375         }
376         switch(p->signbits)
377         {
378         default:
379         case 0: outnear[0] = emaxs[0];outnear[1] = emaxs[1];outnear[2] = emaxs[2];outfar[0] = emins[0];outfar[1] = emins[1];outfar[2] = emins[2];break;
380         case 1: outnear[0] = emins[0];outnear[1] = emaxs[1];outnear[2] = emaxs[2];outfar[0] = emaxs[0];outfar[1] = emins[1];outfar[2] = emins[2];break;
381         case 2: outnear[0] = emaxs[0];outnear[1] = emins[1];outnear[2] = emaxs[2];outfar[0] = emins[0];outfar[1] = emaxs[1];outfar[2] = emins[2];break;
382         case 3: outnear[0] = emins[0];outnear[1] = emins[1];outnear[2] = emaxs[2];outfar[0] = emaxs[0];outfar[1] = emaxs[1];outfar[2] = emins[2];break;
383         case 4: outnear[0] = emaxs[0];outnear[1] = emaxs[1];outnear[2] = emins[2];outfar[0] = emins[0];outfar[1] = emins[1];outfar[2] = emaxs[2];break;
384         case 5: outnear[0] = emins[0];outnear[1] = emaxs[1];outnear[2] = emins[2];outfar[0] = emaxs[0];outfar[1] = emins[1];outfar[2] = emaxs[2];break;
385         case 6: outnear[0] = emaxs[0];outnear[1] = emins[1];outnear[2] = emins[2];outfar[0] = emins[0];outfar[1] = emaxs[1];outfar[2] = emaxs[2];break;
386         case 7: outnear[0] = emins[0];outnear[1] = emins[1];outnear[2] = emins[2];outfar[0] = emaxs[0];outfar[1] = emaxs[1];outfar[2] = emaxs[2];break;
387         }
388 }
389
390 void BoxPlaneCorners_Separate(const vec3_t emins, const vec3_t emaxs, const vec3_t normal, vec3_t outnear, vec3_t outfar)
391 {
392         switch((normal[0] < 0) | ((normal[1] < 0) << 1) | ((normal[2] < 0) << 2))
393         {
394         default:
395         case 0: outnear[0] = emaxs[0];outnear[1] = emaxs[1];outnear[2] = emaxs[2];outfar[0] = emins[0];outfar[1] = emins[1];outfar[2] = emins[2];break;
396         case 1: outnear[0] = emins[0];outnear[1] = emaxs[1];outnear[2] = emaxs[2];outfar[0] = emaxs[0];outfar[1] = emins[1];outfar[2] = emins[2];break;
397         case 2: outnear[0] = emaxs[0];outnear[1] = emins[1];outnear[2] = emaxs[2];outfar[0] = emins[0];outfar[1] = emaxs[1];outfar[2] = emins[2];break;
398         case 3: outnear[0] = emins[0];outnear[1] = emins[1];outnear[2] = emaxs[2];outfar[0] = emaxs[0];outfar[1] = emaxs[1];outfar[2] = emins[2];break;
399         case 4: outnear[0] = emaxs[0];outnear[1] = emaxs[1];outnear[2] = emins[2];outfar[0] = emins[0];outfar[1] = emins[1];outfar[2] = emaxs[2];break;
400         case 5: outnear[0] = emins[0];outnear[1] = emaxs[1];outnear[2] = emins[2];outfar[0] = emaxs[0];outfar[1] = emins[1];outfar[2] = emaxs[2];break;
401         case 6: outnear[0] = emaxs[0];outnear[1] = emins[1];outnear[2] = emins[2];outfar[0] = emins[0];outfar[1] = emaxs[1];outfar[2] = emaxs[2];break;
402         case 7: outnear[0] = emins[0];outnear[1] = emins[1];outnear[2] = emins[2];outfar[0] = emaxs[0];outfar[1] = emaxs[1];outfar[2] = emaxs[2];break;
403         }
404 }
405
406 void BoxPlaneCornerDistances(const vec3_t emins, const vec3_t emaxs, const mplane_t *p, vec_t *outneardist, vec_t *outfardist)
407 {
408         if (p->type < 3)
409         {
410                 *outneardist = emins[p->type] - p->dist;
411                 *outfardist = emaxs[p->type] - p->dist;
412                 return;
413         }
414         switch(p->signbits)
415         {
416         default:
417         case 0: *outneardist = p->normal[0] * emaxs[0] + p->normal[1] * emaxs[1] + p->normal[2] * emaxs[2] - p->dist;*outfardist = p->normal[0] * emins[0] + p->normal[1] * emins[1] + p->normal[2] * emins[2] - p->dist;break;
418         case 1: *outneardist = p->normal[0] * emins[0] + p->normal[1] * emaxs[1] + p->normal[2] * emaxs[2] - p->dist;*outfardist = p->normal[0] * emaxs[0] + p->normal[1] * emins[1] + p->normal[2] * emins[2] - p->dist;break;
419         case 2: *outneardist = p->normal[0] * emaxs[0] + p->normal[1] * emins[1] + p->normal[2] * emaxs[2] - p->dist;*outfardist = p->normal[0] * emins[0] + p->normal[1] * emaxs[1] + p->normal[2] * emins[2] - p->dist;break;
420         case 3: *outneardist = p->normal[0] * emins[0] + p->normal[1] * emins[1] + p->normal[2] * emaxs[2] - p->dist;*outfardist = p->normal[0] * emaxs[0] + p->normal[1] * emaxs[1] + p->normal[2] * emins[2] - p->dist;break;
421         case 4: *outneardist = p->normal[0] * emaxs[0] + p->normal[1] * emaxs[1] + p->normal[2] * emins[2] - p->dist;*outfardist = p->normal[0] * emins[0] + p->normal[1] * emins[1] + p->normal[2] * emaxs[2] - p->dist;break;
422         case 5: *outneardist = p->normal[0] * emins[0] + p->normal[1] * emaxs[1] + p->normal[2] * emins[2] - p->dist;*outfardist = p->normal[0] * emaxs[0] + p->normal[1] * emins[1] + p->normal[2] * emaxs[2] - p->dist;break;
423         case 6: *outneardist = p->normal[0] * emaxs[0] + p->normal[1] * emins[1] + p->normal[2] * emins[2] - p->dist;*outfardist = p->normal[0] * emins[0] + p->normal[1] * emaxs[1] + p->normal[2] * emaxs[2] - p->dist;break;
424         case 7: *outneardist = p->normal[0] * emins[0] + p->normal[1] * emins[1] + p->normal[2] * emins[2] - p->dist;*outfardist = p->normal[0] * emaxs[0] + p->normal[1] * emaxs[1] + p->normal[2] * emaxs[2] - p->dist;break;
425         }
426 }
427
428 void BoxPlaneCornerDistances_Separate(const vec3_t emins, const vec3_t emaxs, const vec3_t normal, vec_t *outneardist, vec_t *outfardist)
429 {
430         switch((normal[0] < 0) | ((normal[1] < 0) << 1) | ((normal[2] < 0) << 2))
431         {
432         default:
433         case 0: *outneardist = normal[0] * emaxs[0] + normal[1] * emaxs[1] + normal[2] * emaxs[2];*outfardist = normal[0] * emins[0] + normal[1] * emins[1] + normal[2] * emins[2];break;
434         case 1: *outneardist = normal[0] * emins[0] + normal[1] * emaxs[1] + normal[2] * emaxs[2];*outfardist = normal[0] * emaxs[0] + normal[1] * emins[1] + normal[2] * emins[2];break;
435         case 2: *outneardist = normal[0] * emaxs[0] + normal[1] * emins[1] + normal[2] * emaxs[2];*outfardist = normal[0] * emins[0] + normal[1] * emaxs[1] + normal[2] * emins[2];break;
436         case 3: *outneardist = normal[0] * emins[0] + normal[1] * emins[1] + normal[2] * emaxs[2];*outfardist = normal[0] * emaxs[0] + normal[1] * emaxs[1] + normal[2] * emins[2];break;
437         case 4: *outneardist = normal[0] * emaxs[0] + normal[1] * emaxs[1] + normal[2] * emins[2];*outfardist = normal[0] * emins[0] + normal[1] * emins[1] + normal[2] * emaxs[2];break;
438         case 5: *outneardist = normal[0] * emins[0] + normal[1] * emaxs[1] + normal[2] * emins[2];*outfardist = normal[0] * emaxs[0] + normal[1] * emins[1] + normal[2] * emaxs[2];break;
439         case 6: *outneardist = normal[0] * emaxs[0] + normal[1] * emins[1] + normal[2] * emins[2];*outfardist = normal[0] * emins[0] + normal[1] * emaxs[1] + normal[2] * emaxs[2];break;
440         case 7: *outneardist = normal[0] * emins[0] + normal[1] * emins[1] + normal[2] * emins[2];*outfardist = normal[0] * emaxs[0] + normal[1] * emaxs[1] + normal[2] * emaxs[2];break;
441         }
442 }
443
444 void AngleVectors (const vec3_t angles, vec3_t forward, vec3_t right, vec3_t up)
445 {
446         double angle, sr, sp, sy, cr, cp, cy;
447
448         angle = angles[YAW] * (M_PI*2 / 360);
449         sy = sin(angle);
450         cy = cos(angle);
451         angle = angles[PITCH] * (M_PI*2 / 360);
452         sp = sin(angle);
453         cp = cos(angle);
454         if (forward)
455         {
456                 forward[0] = cp*cy;
457                 forward[1] = cp*sy;
458                 forward[2] = -sp;
459         }
460         if (right || up)
461         {
462                 if (angles[ROLL])
463                 {
464                         angle = angles[ROLL] * (M_PI*2 / 360);
465                         sr = sin(angle);
466                         cr = cos(angle);
467                         if (right)
468                         {
469                                 right[0] = -1*(sr*sp*cy+cr*-sy);
470                                 right[1] = -1*(sr*sp*sy+cr*cy);
471                                 right[2] = -1*(sr*cp);
472                         }
473                         if (up)
474                         {
475                                 up[0] = (cr*sp*cy+-sr*-sy);
476                                 up[1] = (cr*sp*sy+-sr*cy);
477                                 up[2] = cr*cp;
478                         }
479                 }
480                 else
481                 {
482                         if (right)
483                         {
484                                 right[0] = sy;
485                                 right[1] = -cy;
486                                 right[2] = 0;
487                         }
488                         if (up)
489                         {
490                                 up[0] = (sp*cy);
491                                 up[1] = (sp*sy);
492                                 up[2] = cp;
493                         }
494                 }
495         }
496 }
497
498 void AngleVectorsFLU (const vec3_t angles, vec3_t forward, vec3_t left, vec3_t up)
499 {
500         double angle, sr, sp, sy, cr, cp, cy;
501
502         angle = angles[YAW] * (M_PI*2 / 360);
503         sy = sin(angle);
504         cy = cos(angle);
505         angle = angles[PITCH] * (M_PI*2 / 360);
506         sp = sin(angle);
507         cp = cos(angle);
508         if (forward)
509         {
510                 forward[0] = cp*cy;
511                 forward[1] = cp*sy;
512                 forward[2] = -sp;
513         }
514         if (left || up)
515         {
516                 if (angles[ROLL])
517                 {
518                         angle = angles[ROLL] * (M_PI*2 / 360);
519                         sr = sin(angle);
520                         cr = cos(angle);
521                         if (left)
522                         {
523                                 left[0] = sr*sp*cy+cr*-sy;
524                                 left[1] = sr*sp*sy+cr*cy;
525                                 left[2] = sr*cp;
526                         }
527                         if (up)
528                         {
529                                 up[0] = cr*sp*cy+-sr*-sy;
530                                 up[1] = cr*sp*sy+-sr*cy;
531                                 up[2] = cr*cp;
532                         }
533                 }
534                 else
535                 {
536                         if (left)
537                         {
538                                 left[0] = -sy;
539                                 left[1] = cy;
540                                 left[2] = 0;
541                         }
542                         if (up)
543                         {
544                                 up[0] = sp*cy;
545                                 up[1] = sp*sy;
546                                 up[2] = cp;
547                         }
548                 }
549         }
550 }
551
552 void AngleVectorsDuke3DFLU (const vec3_t angles, vec3_t forward, vec3_t left, vec3_t up, double maxShearAngle)
553 {
554         double angle, sr, sy, cr, cy;
555         double sxx, sxz, szx, szz;
556         double cosMaxShearAngle = cos(maxShearAngle * (M_PI*2 / 360));
557         double tanMaxShearAngle = tan(maxShearAngle * (M_PI*2 / 360));
558
559         angle = angles[YAW] * (M_PI*2 / 360);
560         sy = sin(angle);
561         cy = cos(angle);
562         angle = angles[PITCH] * (M_PI*2 / 360);
563
564         // We will calculate a shear matrix pitch = [[sxx sxz][szx szz]].
565
566         if (fabs(cos(angle)) > cosMaxShearAngle)
567         {
568                 // Pure shear. Keep the original sign of the coefficients.
569                 sxx = 1;
570                 sxz = 0;
571                 szx = -tan(angle);
572                 szz = 1;
573                 // Covering angle per screen coordinate:
574                 // d/dt arctan((sxz + t*szz) / (sxx + t*szx)) @ t=0
575                 // d_angle = det(S) / (sxx*sxx + szx*szx)
576                 //         = 1 / (1 + tan^2 angle)
577                 //         = cos^2 angle.
578         }
579         else
580         {
581                 // A mix of shear and rotation. Implementation-wise, we're
582                 // looking at a capsule, and making the screen surface
583                 // tangential to it... and if we get here, we're looking at the
584                 // two half-spheres of the capsule (and the cylinder part is
585                 // handled above).
586                 double x, y, h, t, d, f;
587                 h = tanMaxShearAngle;
588                 x = cos(angle);
589                 y = sin(angle);
590                 t = h * fabs(y) + sqrt(1 - (h * x) * (h * x));
591                 sxx =  x * t;
592                 sxz =  y * t - h * (y > 0 ? 1.0 : -1.0);
593                 szx = -y * t;
594                 szz =  x * t;
595                 // BUT: keep the amount of a sphere we see in pitch direction
596                 // invariant.
597                 // Covering angle per screen coordinate:
598                 // d_angle = det(S) / (sxx*sxx + szx*szx)
599                 d = (sxx * szz - sxz * szx) / (sxx * sxx + szx * szx);
600                 f = cosMaxShearAngle * cosMaxShearAngle / d;
601                 sxz *= f;
602                 szz *= f;
603         }
604
605         if (forward)
606         {
607                 forward[0] = sxx*cy;
608                 forward[1] = sxx*sy;
609                 forward[2] = szx;
610         }
611         if (left || up)
612         {
613                 if (angles[ROLL])
614                 {
615                         angle = angles[ROLL] * (M_PI*2 / 360);
616                         sr = sin(angle);
617                         cr = cos(angle);
618                         if (left)
619                         {
620                                 left[0] = sr*sxz*cy+cr*-sy;
621                                 left[1] = sr*sxz*sy+cr*cy;
622                                 left[2] = sr*szz;
623                         }
624                         if (up)
625                         {
626                                 up[0] = cr*sxz*cy+-sr*-sy;
627                                 up[1] = cr*sxz*sy+-sr*cy;
628                                 up[2] = cr*szz;
629                         }
630                 }
631                 else
632                 {
633                         if (left)
634                         {
635                                 left[0] = -sy;
636                                 left[1] = cy;
637                                 left[2] = 0;
638                         }
639                         if (up)
640                         {
641                                 up[0] = sxz*cy;
642                                 up[1] = sxz*sy;
643                                 up[2] = szz;
644                         }
645                 }
646         }
647 }
648
649 // LordHavoc: calculates pitch/yaw/roll angles from forward and up vectors
650 void AnglesFromVectors (vec3_t angles, const vec3_t forward, const vec3_t up, qboolean flippitch)
651 {
652         if (forward[0] == 0 && forward[1] == 0)
653         {
654                 if(forward[2] > 0)
655                 {
656                         angles[PITCH] = -M_PI * 0.5;
657                         angles[YAW] = up ? atan2(-up[1], -up[0]) : 0;
658                 }
659                 else
660                 {
661                         angles[PITCH] = M_PI * 0.5;
662                         angles[YAW] = up ? atan2(up[1], up[0]) : 0;
663                 }
664                 angles[ROLL] = 0;
665         }
666         else
667         {
668                 angles[YAW] = atan2(forward[1], forward[0]);
669                 angles[PITCH] = -atan2(forward[2], sqrt(forward[0]*forward[0] + forward[1]*forward[1]));
670                 // note: we know that angles[PITCH] is in ]-pi/2..pi/2[ due to atan2(anything, positive)
671                 if (up)
672                 {
673                         vec_t cp = cos(angles[PITCH]), sp = sin(angles[PITCH]);
674                         // note: we know cp > 0, due to the range angles[pitch] is in
675                         vec_t cy = cos(angles[YAW]), sy = sin(angles[YAW]);
676                         vec3_t tleft, tup;
677                         tleft[0] = -sy;
678                         tleft[1] = cy;
679                         tleft[2] = 0;
680                         tup[0] = sp*cy;
681                         tup[1] = sp*sy;
682                         tup[2] = cp;
683                         angles[ROLL] = -atan2(DotProduct(up, tleft), DotProduct(up, tup));
684                         // for up == '0 0 1', this is
685                         // angles[ROLL] = -atan2(0, cp);
686                         // which is 0
687                 }
688                 else
689                         angles[ROLL] = 0;
690
691                 // so no up vector is equivalent to '1 0 0'!
692         }
693
694         // now convert radians to degrees, and make all values positive
695         VectorScale(angles, 180.0 / M_PI, angles);
696         if (flippitch)
697                 angles[PITCH] *= -1;
698         if (angles[PITCH] < 0) angles[PITCH] += 360;
699         if (angles[YAW] < 0) angles[YAW] += 360;
700         if (angles[ROLL] < 0) angles[ROLL] += 360;
701
702 #if 0
703 {
704         // debugging code
705         vec3_t tforward, tleft, tup, nforward, nup;
706         VectorCopy(forward, nforward);
707         VectorNormalize(nforward);
708         if (up)
709         {
710                 VectorCopy(up, nup);
711                 VectorNormalize(nup);
712                 AngleVectors(angles, tforward, tleft, tup);
713                 if (VectorDistance(tforward, nforward) > 0.01 || VectorDistance(tup, nup) > 0.01)
714                 {
715                         Con_Printf("vectoangles('%f %f %f', '%f %f %f') = %f %f %f\n", nforward[0], nforward[1], nforward[2], nup[0], nup[1], nup[2], angles[0], angles[1], angles[2]);
716                         Con_Printf("^3But that is '%f %f %f', '%f %f %f'\n", tforward[0], tforward[1], tforward[2], tup[0], tup[1], tup[2]);
717                 }
718         }
719         else
720         {
721                 AngleVectors(angles, tforward, tleft, tup);
722                 if (VectorDistance(tforward, nforward) > 0.01)
723                 {
724                         Con_Printf("vectoangles('%f %f %f') = %f %f %f\n", nforward[0], nforward[1], nforward[2], angles[0], angles[1], angles[2]);
725                         Con_Printf("^3But that is '%f %f %f'\n", tforward[0], tforward[1], tforward[2]);
726                 }
727         }
728 }
729 #endif
730 }
731
732 #if 0
733 void AngleMatrix (const vec3_t angles, const vec3_t translate, vec_t matrix[][4])
734 {
735         double angle, sr, sp, sy, cr, cp, cy;
736
737         angle = angles[YAW] * (M_PI*2 / 360);
738         sy = sin(angle);
739         cy = cos(angle);
740         angle = angles[PITCH] * (M_PI*2 / 360);
741         sp = sin(angle);
742         cp = cos(angle);
743         angle = angles[ROLL] * (M_PI*2 / 360);
744         sr = sin(angle);
745         cr = cos(angle);
746         matrix[0][0] = cp*cy;
747         matrix[0][1] = sr*sp*cy+cr*-sy;
748         matrix[0][2] = cr*sp*cy+-sr*-sy;
749         matrix[0][3] = translate[0];
750         matrix[1][0] = cp*sy;
751         matrix[1][1] = sr*sp*sy+cr*cy;
752         matrix[1][2] = cr*sp*sy+-sr*cy;
753         matrix[1][3] = translate[1];
754         matrix[2][0] = -sp;
755         matrix[2][1] = sr*cp;
756         matrix[2][2] = cr*cp;
757         matrix[2][3] = translate[2];
758 }
759 #endif
760
761
762 // LordHavoc: renamed this to Length, and made the normal one a #define
763 float VectorNormalizeLength (vec3_t v)
764 {
765         float length, ilength;
766
767         length = v[0]*v[0] + v[1]*v[1] + v[2]*v[2];
768         length = sqrt (length);
769
770         if (length)
771         {
772                 ilength = 1/length;
773                 v[0] *= ilength;
774                 v[1] *= ilength;
775                 v[2] *= ilength;
776         }
777
778         return length;
779
780 }
781
782
783 /*
784 ================
785 R_ConcatRotations
786 ================
787 */
788 void R_ConcatRotations (const float in1[3*3], const float in2[3*3], float out[3*3])
789 {
790         out[0*3+0] = in1[0*3+0] * in2[0*3+0] + in1[0*3+1] * in2[1*3+0] + in1[0*3+2] * in2[2*3+0];
791         out[0*3+1] = in1[0*3+0] * in2[0*3+1] + in1[0*3+1] * in2[1*3+1] + in1[0*3+2] * in2[2*3+1];
792         out[0*3+2] = in1[0*3+0] * in2[0*3+2] + in1[0*3+1] * in2[1*3+2] + in1[0*3+2] * in2[2*3+2];
793         out[1*3+0] = in1[1*3+0] * in2[0*3+0] + in1[1*3+1] * in2[1*3+0] + in1[1*3+2] * in2[2*3+0];
794         out[1*3+1] = in1[1*3+0] * in2[0*3+1] + in1[1*3+1] * in2[1*3+1] + in1[1*3+2] * in2[2*3+1];
795         out[1*3+2] = in1[1*3+0] * in2[0*3+2] + in1[1*3+1] * in2[1*3+2] + in1[1*3+2] * in2[2*3+2];
796         out[2*3+0] = in1[2*3+0] * in2[0*3+0] + in1[2*3+1] * in2[1*3+0] + in1[2*3+2] * in2[2*3+0];
797         out[2*3+1] = in1[2*3+0] * in2[0*3+1] + in1[2*3+1] * in2[1*3+1] + in1[2*3+2] * in2[2*3+1];
798         out[2*3+2] = in1[2*3+0] * in2[0*3+2] + in1[2*3+1] * in2[1*3+2] + in1[2*3+2] * in2[2*3+2];
799 }
800
801
802 /*
803 ================
804 R_ConcatTransforms
805 ================
806 */
807 void R_ConcatTransforms (const float in1[3*4], const float in2[3*4], float out[3*4])
808 {
809         out[0*4+0] = in1[0*4+0] * in2[0*4+0] + in1[0*4+1] * in2[1*4+0] + in1[0*4+2] * in2[2*4+0];
810         out[0*4+1] = in1[0*4+0] * in2[0*4+1] + in1[0*4+1] * in2[1*4+1] + in1[0*4+2] * in2[2*4+1];
811         out[0*4+2] = in1[0*4+0] * in2[0*4+2] + in1[0*4+1] * in2[1*4+2] + in1[0*4+2] * in2[2*4+2];
812         out[0*4+3] = in1[0*4+0] * in2[0*4+3] + in1[0*4+1] * in2[1*4+3] + in1[0*4+2] * in2[2*4+3] + in1[0*4+3];
813         out[1*4+0] = in1[1*4+0] * in2[0*4+0] + in1[1*4+1] * in2[1*4+0] + in1[1*4+2] * in2[2*4+0];
814         out[1*4+1] = in1[1*4+0] * in2[0*4+1] + in1[1*4+1] * in2[1*4+1] + in1[1*4+2] * in2[2*4+1];
815         out[1*4+2] = in1[1*4+0] * in2[0*4+2] + in1[1*4+1] * in2[1*4+2] + in1[1*4+2] * in2[2*4+2];
816         out[1*4+3] = in1[1*4+0] * in2[0*4+3] + in1[1*4+1] * in2[1*4+3] + in1[1*4+2] * in2[2*4+3] + in1[1*4+3];
817         out[2*4+0] = in1[2*4+0] * in2[0*4+0] + in1[2*4+1] * in2[1*4+0] + in1[2*4+2] * in2[2*4+0];
818         out[2*4+1] = in1[2*4+0] * in2[0*4+1] + in1[2*4+1] * in2[1*4+1] + in1[2*4+2] * in2[2*4+1];
819         out[2*4+2] = in1[2*4+0] * in2[0*4+2] + in1[2*4+1] * in2[1*4+2] + in1[2*4+2] * in2[2*4+2];
820         out[2*4+3] = in1[2*4+0] * in2[0*4+3] + in1[2*4+1] * in2[1*4+3] + in1[2*4+2] * in2[2*4+3] + in1[2*4+3];
821 }
822
823 float RadiusFromBounds (const vec3_t mins, const vec3_t maxs)
824 {
825         vec3_t m1, m2;
826         VectorMultiply(mins, mins, m1);
827         VectorMultiply(maxs, maxs, m2);
828         return sqrt(max(m1[0], m2[0]) + max(m1[1], m2[1]) + max(m1[2], m2[2]));
829 }
830
831 float RadiusFromBoundsAndOrigin (const vec3_t mins, const vec3_t maxs, const vec3_t origin)
832 {
833         vec3_t m1, m2;
834         VectorSubtract(mins, origin, m1);VectorMultiply(m1, m1, m1);
835         VectorSubtract(maxs, origin, m2);VectorMultiply(m2, m2, m2);
836         return sqrt(max(m1[0], m2[0]) + max(m1[1], m2[1]) + max(m1[2], m2[2]));
837 }
838
839 static void Math_RandomSeed_UnitTests(void);
840 void Mathlib_Init(void)
841 {
842         int a;
843
844         // LordHavoc: setup 1.0f / N table for quick recipricols of integers
845         ixtable[0] = 0;
846         for (a = 1;a < 4096;a++)
847                 ixtable[a] = 1.0f / a;
848
849         Math_RandomSeed_UnitTests();
850 }
851
852 #include "matrixlib.h"
853
854 void Matrix4x4_Print(const matrix4x4_t *in)
855 {
856         Con_Printf("%f %f %f %f\n%f %f %f %f\n%f %f %f %f\n%f %f %f %f\n"
857         , in->m[0][0], in->m[0][1], in->m[0][2], in->m[0][3]
858         , in->m[1][0], in->m[1][1], in->m[1][2], in->m[1][3]
859         , in->m[2][0], in->m[2][1], in->m[2][2], in->m[2][3]
860         , in->m[3][0], in->m[3][1], in->m[3][2], in->m[3][3]);
861 }
862
863 int Math_atov(const char *s, prvm_vec3_t out)
864 {
865         int i;
866         VectorClear(out);
867         if (*s == '\'')
868                 s++;
869         for (i = 0;i < 3;i++)
870         {
871                 while (*s == ' ' || *s == '\t')
872                         s++;
873                 out[i] = atof (s);
874                 if (out[i] == 0 && *s != '-' && *s != '+' && (*s < '0' || *s > '9'))
875                         break; // not a number
876                 while (*s && *s != ' ' && *s !='\t' && *s != '\'')
877                         s++;
878                 if (*s == '\'')
879                         break;
880         }
881         return i;
882 }
883
884 void BoxFromPoints(vec3_t mins, vec3_t maxs, int numpoints, vec_t *point3f)
885 {
886         int i;
887         VectorCopy(point3f, mins);
888         VectorCopy(point3f, maxs);
889         for (i = 1, point3f += 3;i < numpoints;i++, point3f += 3)
890         {
891                 mins[0] = min(mins[0], point3f[0]);maxs[0] = max(maxs[0], point3f[0]);
892                 mins[1] = min(mins[1], point3f[1]);maxs[1] = max(maxs[1], point3f[1]);
893                 mins[2] = min(mins[2], point3f[2]);maxs[2] = max(maxs[2], point3f[2]);
894         }
895 }
896
897 // LordHavoc: this has to be done right or you get severe precision breakdown
898 int LoopingFrameNumberFromDouble(double t, int loopframes)
899 {
900         if (loopframes)
901                 return (int)(t - floor(t/loopframes)*loopframes);
902         else
903                 return (int)t;
904 }
905
906 static unsigned int mul_Lecuyer[4] = { 0x12e15e35, 0xb500f16e, 0x2e714eb2, 0xb37916a5 };
907
908 static void mul128(const unsigned int a[], const unsigned int b[], unsigned int dest[4])
909 {
910 #if 0 //defined(__GNUC__) && defined(__x86_64__)
911         unsigned __int128 ia = ((__int128)a[0] << 96) | ((__int128)a[1] << 64) | ((__int128)a[2] << 32) | (a[3]);
912         unsigned __int128 ib = ((__int128)b[0] << 96) | ((__int128)b[1] << 64) | ((__int128)b[2] << 32) | (b[3]);
913         unsigned __int128 id = ia * ib;
914         dest[0] = (id >> 96) & 0xffffffff;
915         dest[1] = (id >> 64) & 0xffffffff;
916         dest[2] = (id >> 32) & 0xffffffff;
917         dest[3] = (id) & 0xffffffff;
918 #else
919         unsigned long long t[4];
920
921         // this multiply chain is relatively straightforward - a[] is repeatedly
922         // added with shifts based on b[] and the results stored into uint64,
923         // but due to C limitations (no access to carry flag) we have to resolve
924         // carries in a really lame way which wastes a fair number of ops
925         // (repeatedly iterating MSB to LSB, rather than LSB to MSB with carry),
926         // an alternative would be to use 16bit multiplies and resolve carries
927         // only at the end, but that would be twice as many multiplies...
928         //
929         // note: >> 32 is a function call in win32 MSVS2015 debug builds.
930         t[0] = (unsigned long long)a[0] * b[3];
931         t[1] = (unsigned long long)a[1] * b[3];
932         t[2] = (unsigned long long)a[2] * b[3];
933         t[3] = (unsigned long long)a[3] * b[3];
934         t[0] += t[1] >> 32;
935         t[1] &= 0xffffffff;
936         t[1] += t[2] >> 32;
937         t[2] &= 0xffffffff;
938         t[2] += t[3] >> 32;
939
940         t[0] += t[1] >> 32;
941         t[1] &= 0xffffffff;
942         t[1] += t[2] >> 32;
943         t[2] &= 0xffffffff;
944
945         t[0] += t[1] >> 32;
946         t[1] &= 0xffffffff;
947
948         t[0] += (unsigned long long)a[1] * b[2];
949         t[1] += (unsigned long long)a[2] * b[2];
950         t[2] += (unsigned long long)a[3] * b[2];
951         t[0] += t[1] >> 32;
952         t[1] &= 0xffffffff;
953         t[1] += t[2] >> 32;
954
955         t[0] += t[1] >> 32;
956         t[1] &= 0xffffffff;
957
958         t[0] += (unsigned long long)a[2] * b[1];
959         t[1] += (unsigned long long)a[3] * b[1];
960         t[0] += t[1] >> 32;
961
962         t[0] += (unsigned long long)a[3] * b[0];
963
964         dest[0] = t[0] & 0xffffffff;
965         dest[1] = t[1] & 0xffffffff;
966         dest[2] = t[2] & 0xffffffff;
967         dest[3] = t[3] & 0xffffffff;
968 #endif
969 }
970
971 static void testmul128(unsigned int a0, unsigned int a1, unsigned int a2, unsigned int a3, unsigned int b0, unsigned int b1, unsigned int b2, unsigned int b3, unsigned int x0, unsigned int x1, unsigned int x2, unsigned int x3)
972 {
973         unsigned int a[4];
974         unsigned int b[4];
975         unsigned int expected[4];
976         unsigned int result[4];
977         a[0] = a0;
978         a[1] = a1;
979         a[2] = a2;
980         a[3] = a3;
981         b[0] = b0;
982         b[1] = b1;
983         b[2] = b2;
984         b[3] = b3;
985         expected[0] = x0;
986         expected[1] = x1;
987         expected[2] = x2;
988         expected[3] = x3;
989         mul128(a, b, result);
990         if (result[0] != expected[0]
991          || result[1] != expected[1]
992          || result[2] != expected[2]
993          || result[3] != expected[3])
994                 Con_Printf("testmul128(\na = %08x %08x %08x %08x,\nb = %08x %08x %08x %08x,\nx = %08x %08x %08x %08x) instead computed\nc = %08x %08x %08x %08x\n", a[0], a[1], a[2], a[3], b[0], b[1], b[2], b[3], expected[0], expected[1], expected[2], expected[3], result[0], result[1], result[2], result[3]);
995 }
996
997 void Math_RandomSeed_UnitTests(void)
998 {
999         testmul128(
1000                 0x00000000, 0x00000000, 0x00000000, 0x00000001,
1001                 0x00000000, 0x00000000, 0x00000000, 0x00000001,
1002                 0x00000000, 0x00000000, 0x00000000, 0x00000001);
1003         testmul128(
1004                 0x00000000, 0x00000000, 0x00000000, 0x00000001,
1005                 0x00000000, 0x00000000, 0x00000001, 0x00000000,
1006                 0x00000000, 0x00000000, 0x00000001, 0x00000000);
1007         testmul128(
1008                 0x00000000, 0x00000000, 0x00000001, 0x00000000,
1009                 0x00000000, 0x00000000, 0x00000000, 0x00000001,
1010                 0x00000000, 0x00000000, 0x00000001, 0x00000000);
1011         testmul128(
1012                 0x00000000, 0x00000000, 0x00000000, 0x00000001,
1013                 0x00000001, 0x00000001, 0x00000001, 0x00000001,
1014                 0x00000001, 0x00000001, 0x00000001, 0x00000001);
1015         testmul128(
1016                 0x00000000, 0x00000000, 0x00000000, 0x00000002,
1017                 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff,
1018                 0xffffffff, 0xffffffff, 0xffffffff, 0xfffffffe);
1019         testmul128(
1020                 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff,
1021                 0x00000000, 0x00000000, 0x00000000, 0x00000002,
1022                 0xffffffff, 0xffffffff, 0xffffffff, 0xfffffffe);
1023         testmul128(
1024                 0x00000000, 0x00000000, 0xffffffff, 0xffffffff,
1025                 0x00000000, 0x00000000, 0x00000002, 0x00000000,
1026                 0x00000001, 0xffffffff, 0xfffffffe, 0x00000000);
1027         testmul128(
1028                 0x00000000, 0x00000000, 0x00000002, 0x00000000,
1029                 0x00000000, 0x00000000, 0xffffffff, 0xffffffff,
1030                 0x00000001, 0xffffffff, 0xfffffffe, 0x00000000);
1031 }
1032
1033 void Math_RandomSeed_Reset(randomseed_t *r)
1034 {
1035         r->s[0] = 1;
1036         r->s[1] = 0;
1037         r->s[2] = 0;
1038         r->s[3] = 0;
1039 }
1040
1041 void Math_RandomSeed_FromInts(randomseed_t *r, unsigned int s0, unsigned int s1, unsigned int s2, unsigned int s3)
1042 {
1043         r->s[0] = s0;
1044         r->s[1] = s1;
1045         r->s[2] = s2;
1046         r->s[3] = s3 | 1; // the Lehmer RNG requires that the seed be odd
1047 }
1048
1049 unsigned long long Math_rand64(randomseed_t *r)
1050 {
1051         unsigned int o[4];
1052         mul128(r->s, mul_Lecuyer, o);
1053         r->s[0] = o[0];
1054         r->s[1] = o[1];
1055         r->s[2] = o[2];
1056         r->s[3] = o[3];
1057         return ((unsigned long long)o[3] << 32) + o[2];
1058 }
1059
1060 float Math_randomf(randomseed_t *r)
1061 {
1062         unsigned long long n = Math_rand64(r);
1063         return n * (0.25f / 0x80000000 / 0x80000000);
1064 }
1065
1066 float Math_crandomf(randomseed_t *r)
1067 {
1068         // do this with a signed number and double the result, so we make use of all parts of the cow
1069         long long n = (long long)Math_rand64(r);
1070         return n * (0.5f / 0x80000000 / 0x80000000);
1071 }
1072
1073 float Math_randomrangef(randomseed_t *r, float minf, float maxf)
1074 {
1075         return Math_randomf(r) * (maxf - minf) + minf;
1076 }
1077
1078 int Math_randomrangei(randomseed_t *r, int mini, int maxi)
1079 {
1080         unsigned long long n = Math_rand64(r);
1081         return (int)(((n >> 33) * (maxi - mini) + mini) >> 31);
1082 }