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         // LadyHavoc: 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 // LadyHavoc: 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 // LadyHavoc: calculates pitch/yaw/roll angles from forward and up vectors
 650 void AnglesFromVectors (vec3_t angles, const vec3_t forward, const vec3_t up, qbool 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 // LadyHavoc: renamed this to Length, and made the normal one a #define
 763 float VectorNormalizeLength (vec3_t v)
 764 {
 765         float length;
 766
 767         length = sqrt(DotProduct(v,v));
 768
 769         if (length)
 770                 VectorScale(v, 1 / length, v);
 771
 772         return length;
 773 }
 774
 775
 776 /*
 777 ================
 778 R_ConcatRotations
 779 ================
 780 */
 781 void R_ConcatRotations (const float in1[3*3], const float in2[3*3], float out[3*3])
 782 {
 783         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];
 784         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];
 785         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];
 786         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];
 787         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];
 788         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];
 789         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];
 790         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];
 791         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];
 792 }
 793
 794
 795 /*
 796 ================
 797 R_ConcatTransforms
 798 ================
 799 */
 800 void R_ConcatTransforms (const float in1[3*4], const float in2[3*4], float out[3*4])
 801 {
 802         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];
 803         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];
 804         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];
 805         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];
 806         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];
 807         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];
 808         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];
 809         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];
 810         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];
 811         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];
 812         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];
 813         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];
 814 }
 815
 816 float RadiusFromBounds (const vec3_t mins, const vec3_t maxs)
 817 {
 818         vec3_t m1, m2;
 819         VectorMultiply(mins, mins, m1);
 820         VectorMultiply(maxs, maxs, m2);
 821         return sqrt(max(m1[0], m2[0]) + max(m1[1], m2[1]) + max(m1[2], m2[2]));
 822 }
 823
 824 float RadiusFromBoundsAndOrigin (const vec3_t mins, const vec3_t maxs, const vec3_t origin)
 825 {
 826         vec3_t m1, m2;
 827         VectorSubtract(mins, origin, m1);VectorMultiply(m1, m1, m1);
 828         VectorSubtract(maxs, origin, m2);VectorMultiply(m2, m2, m2);
 829         return sqrt(max(m1[0], m2[0]) + max(m1[1], m2[1]) + max(m1[2], m2[2]));
 830 }
 831
 832 static void Math_RandomSeed_UnitTests(void);
 833 void Mathlib_Init(void)
 834 {
 835         int a;
 836
 837         // LadyHavoc: setup 1.0f / N table for quick recipricols of integers
 838         ixtable[0] = 0;
 839         for (a = 1;a < 4096;a++)
 840                 ixtable[a] = 1.0f / a;
 841
 842         Math_RandomSeed_UnitTests();
 843 }
 844
 845 #include "matrixlib.h"
 846
 847 void Matrix4x4_Print(const matrix4x4_t *in)
 848 {
 849         Con_Printf("%f %f %f %f\n%f %f %f %f\n%f %f %f %f\n%f %f %f %f\n"
 850         , in->m[0][0], in->m[0][1], in->m[0][2], in->m[0][3]
 851         , in->m[1][0], in->m[1][1], in->m[1][2], in->m[1][3]
 852         , in->m[2][0], in->m[2][1], in->m[2][2], in->m[2][3]
 853         , in->m[3][0], in->m[3][1], in->m[3][2], in->m[3][3]);
 854 }
 855
 856 int Math_atov(const char *s, prvm_vec3_t out)
 857 {
 858         int i;
 859         VectorClear(out);
 860         if (*s == '\'')
 861                 s++;
 862         for (i = 0;i < 3;i++)
 863         {
 864                 while (*s == ' ' || *s == '\t')
 865                         s++;
 866                 out[i] = atof (s);
 867                 if (out[i] == 0 && *s != '-' && *s != '+' && (*s < '0' || *s > '9'))
 868                         break; // not a number
 869                 while (*s && *s != ' ' && *s !='\t' && *s != '\'')
 870                         s++;
 871                 if (*s == '\'')
 872                         break;
 873         }
 874         return i;
 875 }
 876
 877 void BoxFromPoints(vec3_t mins, vec3_t maxs, int numpoints, vec_t *point3f)
 878 {
 879         int i;
 880         VectorCopy(point3f, mins);
 881         VectorCopy(point3f, maxs);
 882         for (i = 1, point3f += 3;i < numpoints;i++, point3f += 3)
 883         {
 884                 mins[0] = min(mins[0], point3f[0]);maxs[0] = max(maxs[0], point3f[0]);
 885                 mins[1] = min(mins[1], point3f[1]);maxs[1] = max(maxs[1], point3f[1]);
 886                 mins[2] = min(mins[2], point3f[2]);maxs[2] = max(maxs[2], point3f[2]);
 887         }
 888 }
 889
 890 // LadyHavoc: this has to be done right or you get severe precision breakdown
 891 int LoopingFrameNumberFromDouble(double t, int loopframes)
 892 {
 893         if (loopframes)
 894                 return (int)(t - floor(t/loopframes)*loopframes);
 895         else
 896                 return (int)t;
 897 }
 898
 899 static unsigned int mul_Lecuyer[4] = { 0x12e15e35, 0xb500f16e, 0x2e714eb2, 0xb37916a5 };
 900
 901 static void mul128(const unsigned int a[], const unsigned int b[], unsigned int dest[4])
 902 {
 903 #if 0 //defined(__GNUC__) && defined(__x86_64__)
 904         unsigned __int128 ia = ((__int128)a[0] << 96) | ((__int128)a[1] << 64) | ((__int128)a[2] << 32) | (a[3]);
 905         unsigned __int128 ib = ((__int128)b[0] << 96) | ((__int128)b[1] << 64) | ((__int128)b[2] << 32) | (b[3]);
 906         unsigned __int128 id = ia * ib;
 907         dest[0] = (id >> 96) & 0xffffffff;
 908         dest[1] = (id >> 64) & 0xffffffff;
 909         dest[2] = (id >> 32) & 0xffffffff;
 910         dest[3] = (id) & 0xffffffff;
 911 #else
 912         unsigned long long t[4];
 913
 914         // this multiply chain is relatively straightforward - a[] is repeatedly
 915         // added with shifts based on b[] and the results stored into uint64,
 916         // but due to C limitations (no access to carry flag) we have to resolve
 917         // carries in a really lame way which wastes a fair number of ops
 918         // (repeatedly iterating MSB to LSB, rather than LSB to MSB with carry),
 919         // an alternative would be to use 16bit multiplies and resolve carries
 920         // only at the end, but that would be twice as many multiplies...
 921         //
 922         // note: >> 32 is a function call in win32 MSVS2015 debug builds.
 923         t[0] = (unsigned long long)a[0] * b[3];
 924         t[1] = (unsigned long long)a[1] * b[3];
 925         t[2] = (unsigned long long)a[2] * b[3];
 926         t[3] = (unsigned long long)a[3] * b[3];
 927         t[0] += t[1] >> 32;
 928         t[1] &= 0xffffffff;
 929         t[1] += t[2] >> 32;
 930         t[2] &= 0xffffffff;
 931         t[2] += t[3] >> 32;
 932
 933         t[0] += t[1] >> 32;
 934         t[1] &= 0xffffffff;
 935         t[1] += t[2] >> 32;
 936         t[2] &= 0xffffffff;
 937
 938         t[0] += t[1] >> 32;
 939         t[1] &= 0xffffffff;
 940
 941         t[0] += (unsigned long long)a[1] * b[2];
 942         t[1] += (unsigned long long)a[2] * b[2];
 943         t[2] += (unsigned long long)a[3] * b[2];
 944         t[0] += t[1] >> 32;
 945         t[1] &= 0xffffffff;
 946         t[1] += t[2] >> 32;
 947
 948         t[0] += t[1] >> 32;
 949         t[1] &= 0xffffffff;
 950
 951         t[0] += (unsigned long long)a[2] * b[1];
 952         t[1] += (unsigned long long)a[3] * b[1];
 953         t[0] += t[1] >> 32;
 954
 955         t[0] += (unsigned long long)a[3] * b[0];
 956
 957         dest[0] = t[0] & 0xffffffff;
 958         dest[1] = t[1] & 0xffffffff;
 959         dest[2] = t[2] & 0xffffffff;
 960         dest[3] = t[3] & 0xffffffff;
 961 #endif
 962 }
 963
 964 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)
 965 {
 966         unsigned int a[4];
 967         unsigned int b[4];
 968         unsigned int expected[4];
 969         unsigned int result[4];
 970         a[0] = a0;
 971         a[1] = a1;
 972         a[2] = a2;
 973         a[3] = a3;
 974         b[0] = b0;
 975         b[1] = b1;
 976         b[2] = b2;
 977         b[3] = b3;
 978         expected[0] = x0;
 979         expected[1] = x1;
 980         expected[2] = x2;
 981         expected[3] = x3;
 982         mul128(a, b, result);
 983         if (result[0] != expected[0]
 984          || result[1] != expected[1]
 985          || result[2] != expected[2]
 986          || result[3] != expected[3])
 987                 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]);
 988 }
 989
 990 void Math_RandomSeed_UnitTests(void)
 991 {
 992         testmul128(
 993                 0x00000000, 0x00000000, 0x00000000, 0x00000001,
 994                 0x00000000, 0x00000000, 0x00000000, 0x00000001,
 995                 0x00000000, 0x00000000, 0x00000000, 0x00000001);
 996         testmul128(
 997                 0x00000000, 0x00000000, 0x00000000, 0x00000001,
 998                 0x00000000, 0x00000000, 0x00000001, 0x00000000,
 999                 0x00000000, 0x00000000, 0x00000001, 0x00000000);
1000         testmul128(
1001                 0x00000000, 0x00000000, 0x00000001, 0x00000000,
1002                 0x00000000, 0x00000000, 0x00000000, 0x00000001,
1003                 0x00000000, 0x00000000, 0x00000001, 0x00000000);
1004         testmul128(
1005                 0x00000000, 0x00000000, 0x00000000, 0x00000001,
1006                 0x00000001, 0x00000001, 0x00000001, 0x00000001,
1007                 0x00000001, 0x00000001, 0x00000001, 0x00000001);
1008         testmul128(
1009                 0x00000000, 0x00000000, 0x00000000, 0x00000002,
1010                 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff,
1011                 0xffffffff, 0xffffffff, 0xffffffff, 0xfffffffe);
1012         testmul128(
1013                 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff,
1014                 0x00000000, 0x00000000, 0x00000000, 0x00000002,
1015                 0xffffffff, 0xffffffff, 0xffffffff, 0xfffffffe);
1016         testmul128(
1017                 0x00000000, 0x00000000, 0xffffffff, 0xffffffff,
1018                 0x00000000, 0x00000000, 0x00000002, 0x00000000,
1019                 0x00000001, 0xffffffff, 0xfffffffe, 0x00000000);
1020         testmul128(
1021                 0x00000000, 0x00000000, 0x00000002, 0x00000000,
1022                 0x00000000, 0x00000000, 0xffffffff, 0xffffffff,
1023                 0x00000001, 0xffffffff, 0xfffffffe, 0x00000000);
1024 }
1025
1026 void Math_RandomSeed_Reset(randomseed_t *r)
1027 {
1028         r->s[0] = 1;
1029         r->s[1] = 0;
1030         r->s[2] = 0;
1031         r->s[3] = 0;
1032 }
1033
1034 void Math_RandomSeed_FromInts(randomseed_t *r, unsigned int s0, unsigned int s1, unsigned int s2, unsigned int s3)
1035 {
1036         r->s[0] = s0;
1037         r->s[1] = s1;
1038         r->s[2] = s2;
1039         r->s[3] = s3 | 1; // the Lehmer RNG requires that the seed be odd
1040 }
1041
1042 unsigned long long Math_rand64(randomseed_t *r)
1043 {
1044         unsigned int o[4];
1045         mul128(r->s, mul_Lecuyer, o);
1046         r->s[0] = o[0];
1047         r->s[1] = o[1];
1048         r->s[2] = o[2];
1049         r->s[3] = o[3];
1050         return ((unsigned long long)o[3] << 32) + o[2];
1051 }
1052
1053 float Math_randomf(randomseed_t *r)
1054 {
1055         unsigned long long n = Math_rand64(r);
1056         return n * (0.25f / 0x80000000 / 0x80000000);
1057 }
1058
1059 float Math_crandomf(randomseed_t *r)
1060 {
1061         // do this with a signed number and double the result, so we make use of all parts of the cow
1062         long long n = (long long)Math_rand64(r);
1063         return n * (0.5f / 0x80000000 / 0x80000000);
1064 }
1065
1066 float Math_randomrangef(randomseed_t *r, float minf, float maxf)
1067 {
1068         return Math_randomf(r) * (maxf - minf) + minf;
1069 }
1070
1071 int Math_randomrangei(randomseed_t *r, int mini, int maxi)
1072 {
1073         unsigned long long n = Math_rand64(r);
1074         return (int)(((n >> 33) * (maxi - mini) + mini) >> 31);
1075 }