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 <math.h>
  23 #include "quakedef.h"
  24
  25 void Sys_Error (char *error, ...);
  26
  27 vec3_t vec3_origin = {0,0,0};
  28 int nanmask = 255<<23;
  29 float ixtable[4096];
  30
  31 /*-----------------------------------------------------------------*/
  32
  33 float m_bytenormals[NUMVERTEXNORMALS][3] =
  34 {
  35 {-0.525731, 0.000000, 0.850651}, {-0.442863, 0.238856, 0.864188},
  36 {-0.295242, 0.000000, 0.955423}, {-0.309017, 0.500000, 0.809017},
  37 {-0.162460, 0.262866, 0.951056}, {0.000000, 0.000000, 1.000000},
  38 {0.000000, 0.850651, 0.525731}, {-0.147621, 0.716567, 0.681718},
  39 {0.147621, 0.716567, 0.681718}, {0.000000, 0.525731, 0.850651},
  40 {0.309017, 0.500000, 0.809017}, {0.525731, 0.000000, 0.850651},
  41 {0.295242, 0.000000, 0.955423}, {0.442863, 0.238856, 0.864188},
  42 {0.162460, 0.262866, 0.951056}, {-0.681718, 0.147621, 0.716567},
  43 {-0.809017, 0.309017, 0.500000}, {-0.587785, 0.425325, 0.688191},
  44 {-0.850651, 0.525731, 0.000000}, {-0.864188, 0.442863, 0.238856},
  45 {-0.716567, 0.681718, 0.147621}, {-0.688191, 0.587785, 0.425325},
  46 {-0.500000, 0.809017, 0.309017}, {-0.238856, 0.864188, 0.442863},
  47 {-0.425325, 0.688191, 0.587785}, {-0.716567, 0.681718, -0.147621},
  48 {-0.500000, 0.809017, -0.309017}, {-0.525731, 0.850651, 0.000000},
  49 {0.000000, 0.850651, -0.525731}, {-0.238856, 0.864188, -0.442863},
  50 {0.000000, 0.955423, -0.295242}, {-0.262866, 0.951056, -0.162460},
  51 {0.000000, 1.000000, 0.000000}, {0.000000, 0.955423, 0.295242},
  52 {-0.262866, 0.951056, 0.162460}, {0.238856, 0.864188, 0.442863},
  53 {0.262866, 0.951056, 0.162460}, {0.500000, 0.809017, 0.309017},
  54 {0.238856, 0.864188, -0.442863}, {0.262866, 0.951056, -0.162460},
  55 {0.500000, 0.809017, -0.309017}, {0.850651, 0.525731, 0.000000},
  56 {0.716567, 0.681718, 0.147621}, {0.716567, 0.681718, -0.147621},
  57 {0.525731, 0.850651, 0.000000}, {0.425325, 0.688191, 0.587785},
  58 {0.864188, 0.442863, 0.238856}, {0.688191, 0.587785, 0.425325},
  59 {0.809017, 0.309017, 0.500000}, {0.681718, 0.147621, 0.716567},
  60 {0.587785, 0.425325, 0.688191}, {0.955423, 0.295242, 0.000000},
  61 {1.000000, 0.000000, 0.000000}, {0.951056, 0.162460, 0.262866},
  62 {0.850651, -0.525731, 0.000000}, {0.955423, -0.295242, 0.000000},
  63 {0.864188, -0.442863, 0.238856}, {0.951056, -0.162460, 0.262866},
  64 {0.809017, -0.309017, 0.500000}, {0.681718, -0.147621, 0.716567},
  65 {0.850651, 0.000000, 0.525731}, {0.864188, 0.442863, -0.238856},
  66 {0.809017, 0.309017, -0.500000}, {0.951056, 0.162460, -0.262866},
  67 {0.525731, 0.000000, -0.850651}, {0.681718, 0.147621, -0.716567},
  68 {0.681718, -0.147621, -0.716567}, {0.850651, 0.000000, -0.525731},
  69 {0.809017, -0.309017, -0.500000}, {0.864188, -0.442863, -0.238856},
  70 {0.951056, -0.162460, -0.262866}, {0.147621, 0.716567, -0.681718},
  71 {0.309017, 0.500000, -0.809017}, {0.425325, 0.688191, -0.587785},
  72 {0.442863, 0.238856, -0.864188}, {0.587785, 0.425325, -0.688191},
  73 {0.688191, 0.587785, -0.425325}, {-0.147621, 0.716567, -0.681718},
  74 {-0.309017, 0.500000, -0.809017}, {0.000000, 0.525731, -0.850651},
  75 {-0.525731, 0.000000, -0.850651}, {-0.442863, 0.238856, -0.864188},
  76 {-0.295242, 0.000000, -0.955423}, {-0.162460, 0.262866, -0.951056},
  77 {0.000000, 0.000000, -1.000000}, {0.295242, 0.000000, -0.955423},
  78 {0.162460, 0.262866, -0.951056}, {-0.442863, -0.238856, -0.864188},
  79 {-0.309017, -0.500000, -0.809017}, {-0.162460, -0.262866, -0.951056},
  80 {0.000000, -0.850651, -0.525731}, {-0.147621, -0.716567, -0.681718},
  81 {0.147621, -0.716567, -0.681718}, {0.000000, -0.525731, -0.850651},
  82 {0.309017, -0.500000, -0.809017}, {0.442863, -0.238856, -0.864188},
  83 {0.162460, -0.262866, -0.951056}, {0.238856, -0.864188, -0.442863},
  84 {0.500000, -0.809017, -0.309017}, {0.425325, -0.688191, -0.587785},
  85 {0.716567, -0.681718, -0.147621}, {0.688191, -0.587785, -0.425325},
  86 {0.587785, -0.425325, -0.688191}, {0.000000, -0.955423, -0.295242},
  87 {0.000000, -1.000000, 0.000000}, {0.262866, -0.951056, -0.162460},
  88 {0.000000, -0.850651, 0.525731}, {0.000000, -0.955423, 0.295242},
  89 {0.238856, -0.864188, 0.442863}, {0.262866, -0.951056, 0.162460},
  90 {0.500000, -0.809017, 0.309017}, {0.716567, -0.681718, 0.147621},
  91 {0.525731, -0.850651, 0.000000}, {-0.238856, -0.864188, -0.442863},
  92 {-0.500000, -0.809017, -0.309017}, {-0.262866, -0.951056, -0.162460},
  93 {-0.850651, -0.525731, 0.000000}, {-0.716567, -0.681718, -0.147621},
  94 {-0.716567, -0.681718, 0.147621}, {-0.525731, -0.850651, 0.000000},
  95 {-0.500000, -0.809017, 0.309017}, {-0.238856, -0.864188, 0.442863},
  96 {-0.262866, -0.951056, 0.162460}, {-0.864188, -0.442863, 0.238856},
  97 {-0.809017, -0.309017, 0.500000}, {-0.688191, -0.587785, 0.425325},
  98 {-0.681718, -0.147621, 0.716567}, {-0.442863, -0.238856, 0.864188},
  99 {-0.587785, -0.425325, 0.688191}, {-0.309017, -0.500000, 0.809017},
 100 {-0.147621, -0.716567, 0.681718}, {-0.425325, -0.688191, 0.587785},
 101 {-0.162460, -0.262866, 0.951056}, {0.442863, -0.238856, 0.864188},
 102 {0.162460, -0.262866, 0.951056}, {0.309017, -0.500000, 0.809017},
 103 {0.147621, -0.716567, 0.681718}, {0.000000, -0.525731, 0.850651},
 104 {0.425325, -0.688191, 0.587785}, {0.587785, -0.425325, 0.688191},
 105 {0.688191, -0.587785, 0.425325}, {-0.955423, 0.295242, 0.000000},
 106 {-0.951056, 0.162460, 0.262866}, {-1.000000, 0.000000, 0.000000},
 107 {-0.850651, 0.000000, 0.525731}, {-0.955423, -0.295242, 0.000000},
 108 {-0.951056, -0.162460, 0.262866}, {-0.864188, 0.442863, -0.238856},
 109 {-0.951056, 0.162460, -0.262866}, {-0.809017, 0.309017, -0.500000},
 110 {-0.864188, -0.442863, -0.238856}, {-0.951056, -0.162460, -0.262866},
 111 {-0.809017, -0.309017, -0.500000}, {-0.681718, 0.147621, -0.716567},
 112 {-0.681718, -0.147621, -0.716567}, {-0.850651, 0.000000, -0.525731},
 113 {-0.688191, 0.587785, -0.425325}, {-0.587785, 0.425325, -0.688191},
 114 {-0.425325, 0.688191, -0.587785}, {-0.425325, -0.688191, -0.587785},
 115 {-0.587785, -0.425325, -0.688191}, {-0.688191, -0.587785, -0.425325},
 116 };
 117
 118 byte NormalToByte(vec3_t n)
 119 {
 120         int i, best;
 121         float bestdistance, distance;
 122
 123         best = 0;
 124         bestdistance = DotProduct (n, m_bytenormals[0]);
 125         for (i = 1;i < NUMVERTEXNORMALS;i++)
 126         {
 127                 distance = DotProduct (n, m_bytenormals[i]);
 128                 if (distance > bestdistance)
 129                 {
 130                         bestdistance = distance;
 131                         best = i;
 132                 }
 133         }
 134         return best;
 135 }
 136
 137 // note: uses byte partly to force unsigned for the validity check
 138 void ByteToNormal(byte num, vec3_t n)
 139 {
 140         if (num < NUMVERTEXNORMALS)
 141                 VectorCopy(m_bytenormals[num], n)
 142         else
 143                 VectorClear(n) // FIXME: complain?
 144 }
 145
 146 float Q_RSqrt(float number)
 147 {
 148         float y;
 149
 150         if (number == 0.0f)
 151                 return 0.0f;
 152
 153         *((long *)&y) = 0x5f3759df - ((* (long *) &number) >> 1);
 154         return y * (1.5f - (number * 0.5f * y * y));
 155 }
 156
 157 void _VectorNormalizeFast(vec3_t v)
 158 {
 159         float y, number;
 160
 161         number = DotProduct(v, v);
 162
 163         if (number != 0.0)
 164         {
 165                 *((long *)&y) = 0x5f3759df - ((* (long *) &number) >> 1);
 166                 y = y * (1.5f - (number * 0.5f * y * y));
 167
 168                 VectorScale(v, y, v);
 169         }
 170 }
 171
 172 #if 0
 173 // LordHavoc: no longer used at all
 174 void ProjectPointOnPlane( vec3_t dst, const vec3_t p, const vec3_t normal )
 175 {
 176 #if 0
 177         // LordHavoc: the old way...
 178         float d;
 179         vec3_t n;
 180         float inv_denom;
 181
 182         inv_denom = 1.0F / DotProduct( normal, normal );
 183
 184         d = DotProduct( normal, p ) * inv_denom;
 185
 186         n[0] = normal[0] * inv_denom;
 187         n[1] = normal[1] * inv_denom;
 188         n[2] = normal[2] * inv_denom;
 189
 190         dst[0] = p[0] - d * n[0];
 191         dst[1] = p[1] - d * n[1];
 192         dst[2] = p[2] - d * n[2];
 193 #else
 194         // LordHavoc: optimized to death and beyond
 195         float d;
 196
 197         // LordHavoc: the normal is a unit vector by definition,
 198         //            therefore inv_denom was pointless.
 199         d = DotProduct(normal, p);
 200         dst[0] = p[0] - d * normal[0];
 201         dst[1] = p[1] - d * normal[1];
 202         dst[2] = p[2] - d * normal[2];
 203 #endif
 204 }
 205 #endif
 206
 207 // assumes "src" is normalized
 208 void PerpendicularVector( vec3_t dst, const vec3_t src )
 209 {
 210 #if 0
 211         // LordHavoc: the old way...
 212         int     pos;
 213         int i;
 214         float minelem, d;
 215         vec3_t tempvec;
 216
 217         // find the smallest magnitude axially aligned vector
 218         minelem = 1.0F;
 219         for ( pos = 0, i = 0; i < 3; i++ )
 220         {
 221                 if ( fabs( src[i] ) < minelem )
 222                 {
 223                         pos = i;
 224                         minelem = fabs( src[i] );
 225                 }
 226         }
 227         VectorClear(tempvec);
 228         tempvec[pos] = 1.0F;
 229
 230         // project the point onto the plane defined by src
 231         ProjectPointOnPlane( dst, tempvec, src );
 232
 233         // normalize the result
 234         VectorNormalize(dst);
 235 #else
 236         // LordHavoc: optimized to death and beyond
 237         int     pos;
 238         float minelem;
 239
 240         if (src[0])
 241         {
 242                 dst[0] = 0;
 243                 if (src[1])
 244                 {
 245                         dst[1] = 0;
 246                         if (src[2])
 247                         {
 248                                 dst[2] = 0;
 249                                 pos = 0;
 250                                 minelem = fabs(src[0]);
 251                                 if (fabs(src[1]) < minelem)
 252                                 {
 253                                         pos = 1;
 254                                         minelem = fabs(src[1]);
 255                                 }
 256                                 if (fabs(src[2]) < minelem)
 257                                         pos = 2;
 258
 259                                 dst[pos] = 1;
 260                                 dst[0] -= src[pos] * src[0];
 261                                 dst[1] -= src[pos] * src[1];
 262                                 dst[2] -= src[pos] * src[2];
 263
 264                                 // normalize the result
 265                                 VectorNormalize(dst);
 266                         }
 267                         else
 268                                 dst[2] = 1;
 269                 }
 270                 else
 271                 {
 272                         dst[1] = 1;
 273                         dst[2] = 0;
 274                 }
 275         }
 276         else
 277         {
 278                 dst[0] = 1;
 279                 dst[1] = 0;
 280                 dst[2] = 0;
 281         }
 282 #endif
 283 }
 284
 285
 286 #ifdef _WIN32
 287 #pragma optimize( "", off )
 288 #endif
 289
 290
 291 // LordHavoc: like AngleVectors, but taking a forward vector instead of angles, useful!
 292 void VectorVectors(const vec3_t forward, vec3_t right, vec3_t up)
 293 {
 294         float d;
 295
 296         right[0] = forward[2];
 297         right[1] = -forward[0];
 298         right[2] = forward[1];
 299
 300         d = DotProduct(forward, right);
 301         right[0] -= d * forward[0];
 302         right[1] -= d * forward[1];
 303         right[2] -= d * forward[2];
 304         VectorNormalize(right);
 305         CrossProduct(right, forward, up);
 306 }
 307
 308 void RotatePointAroundVector( vec3_t dst, const vec3_t dir, const vec3_t point, float degrees )
 309 {
 310 #if 0
 311         // LordHavoc: the old way, cryptic brute force...
 312         float   m[3][3];
 313         float   im[3][3];
 314         float   zrot[3][3];
 315         float   tmpmat[3][3];
 316         float   rot[3][3];
 317         int     i;
 318         vec3_t vr, vup, vf;
 319
 320         vf[0] = dir[0];
 321         vf[1] = dir[1];
 322         vf[2] = dir[2];
 323
 324         PerpendicularVector( vr, dir );
 325         CrossProduct( vr, vf, vup );
 326
 327         m[0][0] = vr[0];
 328         m[1][0] = vr[1];
 329         m[2][0] = vr[2];
 330
 331         m[0][1] = vup[0];
 332         m[1][1] = vup[1];
 333         m[2][1] = vup[2];
 334
 335         m[0][2] = vf[0];
 336         m[1][2] = vf[1];
 337         m[2][2] = vf[2];
 338
 339         memcpy( im, m, sizeof( im ) );
 340
 341         im[0][1] = m[1][0];
 342         im[0][2] = m[2][0];
 343         im[1][0] = m[0][1];
 344         im[1][2] = m[2][1];
 345         im[2][0] = m[0][2];
 346         im[2][1] = m[1][2];
 347
 348         memset( zrot, 0, sizeof( zrot ) );
 349         zrot[0][0] = zrot[1][1] = zrot[2][2] = 1.0F;
 350
 351         zrot[0][0] = cos( DEG2RAD( degrees ) );
 352         zrot[0][1] = sin( DEG2RAD( degrees ) );
 353         zrot[1][0] = -sin( DEG2RAD( degrees ) );
 354         zrot[1][1] = cos( DEG2RAD( degrees ) );
 355
 356         R_ConcatRotations( m, zrot, tmpmat );
 357         R_ConcatRotations( tmpmat, im, rot );
 358
 359         for ( i = 0; i < 3; i++ )
 360                 dst[i] = rot[i][0] * point[0] + rot[i][1] * point[1] + rot[i][2] * point[2];
 361 #elif 0
 362         // LordHavoc: on the path to unintelligible code...
 363 //      float   m[3][3];
 364 //      float   im[3][3];
 365 //      float   zrot[3][3];
 366         float   tmpmat[3][3];
 367 //      float   rot[3][3];
 368         float   angle, c, s;
 369 //      int     i;
 370         vec3_t vr, vu, vf;
 371
 372         angle = DEG2RAD(degrees);
 373
 374         c = cos(angle);
 375         s = sin(angle);
 376
 377         vf[0] = dir[0];
 378         vf[1] = dir[1];
 379         vf[2] = dir[2];
 380
 381         PerpendicularVector(vr, dir);
 382         CrossProduct(vr, vf, vu);
 383
 384 //      m   [0][0] = vr[0];m   [0][1] = vu[0];m   [0][2] = vf[0];
 385 //      m   [1][0] = vr[1];m   [1][1] = vu[1];m   [1][2] = vf[1];
 386 //      m   [2][0] = vr[2];m   [2][1] = vu[2];m   [2][2] = vf[2];
 387 //      im  [0][0] = vr[0];im  [0][1] = vr[1];im  [0][2] = vr[2];
 388 //      im  [1][0] = vu[0];im  [1][1] = vu[1];im  [1][2] = vu[2];
 389 //      im  [2][0] = vf[0];im  [2][1] = vf[1];im  [2][2] = vf[2];
 390 //      zrot[0][0] =     c;zrot[0][1] =     s;zrot[0][2] =     0;
 391 //      zrot[1][0] =    -s;zrot[1][1] =     c;zrot[1][2] =     0;
 392 //      zrot[2][0] =     0;zrot[2][1] =     0;zrot[2][2] =     1;
 393
 394 //      tmpmat[0][0] = m[0][0] * zrot[0][0] + m[0][1] * zrot[1][0] + m[0][2] * zrot[2][0];
 395 //      tmpmat[0][1] = m[0][0] * zrot[0][1] + m[0][1] * zrot[1][1] + m[0][2] * zrot[2][1];
 396 //      tmpmat[0][2] = m[0][0] * zrot[0][2] + m[0][1] * zrot[1][2] + m[0][2] * zrot[2][2];
 397 //      tmpmat[1][0] = m[1][0] * zrot[0][0] + m[1][1] * zrot[1][0] + m[1][2] * zrot[2][0];
 398 //      tmpmat[1][1] = m[1][0] * zrot[0][1] + m[1][1] * zrot[1][1] + m[1][2] * zrot[2][1];
 399 //      tmpmat[1][2] = m[1][0] * zrot[0][2] + m[1][1] * zrot[1][2] + m[1][2] * zrot[2][2];
 400 //      tmpmat[2][0] = m[2][0] * zrot[0][0] + m[2][1] * zrot[1][0] + m[2][2] * zrot[2][0];
 401 //      tmpmat[2][1] = m[2][0] * zrot[0][1] + m[2][1] * zrot[1][1] + m[2][2] * zrot[2][1];
 402 //      tmpmat[2][2] = m[2][0] * zrot[0][2] + m[2][1] * zrot[1][2] + m[2][2] * zrot[2][2];
 403
 404         tmpmat[0][0] = vr[0] *  c + vu[0] * -s;
 405         tmpmat[0][1] = vr[0] *  s + vu[0] *  c;
 406 //      tmpmat[0][2] =                           vf[0];
 407         tmpmat[1][0] = vr[1] *  c + vu[1] * -s;
 408         tmpmat[1][1] = vr[1] *  s + vu[1] *  c;
 409 //      tmpmat[1][2] =                           vf[1];
 410         tmpmat[2][0] = vr[2] *  c + vu[2] * -s;
 411         tmpmat[2][1] = vr[2] *  s + vu[2] *  c;
 412 //      tmpmat[2][2] =                           vf[2];
 413
 414 //      rot[0][0] = tmpmat[0][0] * vr[0] + tmpmat[0][1] * vu[0] + tmpmat[0][2] * vf[0];
 415 //      rot[0][1] = tmpmat[0][0] * vr[1] + tmpmat[0][1] * vu[1] + tmpmat[0][2] * vf[1];
 416 //      rot[0][2] = tmpmat[0][0] * vr[2] + tmpmat[0][1] * vu[2] + tmpmat[0][2] * vf[2];
 417 //      rot[1][0] = tmpmat[1][0] * vr[0] + tmpmat[1][1] * vu[0] + tmpmat[1][2] * vf[0];
 418 //      rot[1][1] = tmpmat[1][0] * vr[1] + tmpmat[1][1] * vu[1] + tmpmat[1][2] * vf[1];
 419 //      rot[1][2] = tmpmat[1][0] * vr[2] + tmpmat[1][1] * vu[2] + tmpmat[1][2] * vf[2];
 420 //      rot[2][0] = tmpmat[2][0] * vr[0] + tmpmat[2][1] * vu[0] + tmpmat[2][2] * vf[0];
 421 //      rot[2][1] = tmpmat[2][0] * vr[1] + tmpmat[2][1] * vu[1] + tmpmat[2][2] * vf[1];
 422 //      rot[2][2] = tmpmat[2][0] * vr[2] + tmpmat[2][1] * vu[2] + tmpmat[2][2] * vf[2];
 423
 424 //      rot[0][0] = tmpmat[0][0] * vr[0] + tmpmat[0][1] * vu[0] + vf[0] * vf[0];
 425 //      rot[0][1] = tmpmat[0][0] * vr[1] + tmpmat[0][1] * vu[1] + vf[0] * vf[1];
 426 //      rot[0][2] = tmpmat[0][0] * vr[2] + tmpmat[0][1] * vu[2] + vf[0] * vf[2];
 427 //      rot[1][0] = tmpmat[1][0] * vr[0] + tmpmat[1][1] * vu[0] + vf[1] * vf[0];
 428 //      rot[1][1] = tmpmat[1][0] * vr[1] + tmpmat[1][1] * vu[1] + vf[1] * vf[1];
 429 //      rot[1][2] = tmpmat[1][0] * vr[2] + tmpmat[1][1] * vu[2] + vf[1] * vf[2];
 430 //      rot[2][0] = tmpmat[2][0] * vr[0] + tmpmat[2][1] * vu[0] + vf[2] * vf[0];
 431 //      rot[2][1] = tmpmat[2][0] * vr[1] + tmpmat[2][1] * vu[1] + vf[2] * vf[1];
 432 //      rot[2][2] = tmpmat[2][0] * vr[2] + tmpmat[2][1] * vu[2] + vf[2] * vf[2];
 433
 434 //      dst[0] = rot[0][0] * point[0] + rot[0][1] * point[1] + rot[0][2] * point[2];
 435 //      dst[1] = rot[1][0] * point[0] + rot[1][1] * point[1] + rot[1][2] * point[2];
 436 //      dst[2] = rot[2][0] * point[0] + rot[2][1] * point[1] + rot[2][2] * point[2];
 437
 438         dst[0] = (tmpmat[0][0] * vr[0] + tmpmat[0][1] * vu[0] + vf[0] * vf[0]) * point[0]
 439                + (tmpmat[0][0] * vr[1] + tmpmat[0][1] * vu[1] + vf[0] * vf[1]) * point[1]
 440                + (tmpmat[0][0] * vr[2] + tmpmat[0][1] * vu[2] + vf[0] * vf[2]) * point[2];
 441         dst[1] = (tmpmat[1][0] * vr[0] + tmpmat[1][1] * vu[0] + vf[1] * vf[0]) * point[0]
 442                + (tmpmat[1][0] * vr[1] + tmpmat[1][1] * vu[1] + vf[1] * vf[1]) * point[1]
 443                + (tmpmat[1][0] * vr[2] + tmpmat[1][1] * vu[2] + vf[1] * vf[2]) * point[2];
 444         dst[2] = (tmpmat[2][0] * vr[0] + tmpmat[2][1] * vu[0] + vf[2] * vf[0]) * point[0]
 445                + (tmpmat[2][0] * vr[1] + tmpmat[2][1] * vu[1] + vf[2] * vf[1]) * point[1]
 446                + (tmpmat[2][0] * vr[2] + tmpmat[2][1] * vu[2] + vf[2] * vf[2]) * point[2];
 447 #else
 448         // LordHavoc: optimized to death and beyond, cryptic in an entirely new way
 449         float   t0, t1;
 450         float   angle, c, s;
 451         vec3_t  vr, vu, vf;
 452
 453         angle = DEG2RAD(degrees);
 454
 455         c = cos(angle);
 456         s = sin(angle);
 457
 458         vf[0] = dir[0];
 459         vf[1] = dir[1];
 460         vf[2] = dir[2];
 461
 462 //      PerpendicularVector(vr, dir);
 463 //      CrossProduct(vr, vf, vu);
 464         VectorVectors(vf, vr, vu);
 465
 466         t0 = vr[0] *  c + vu[0] * -s;
 467         t1 = vr[0] *  s + vu[0] *  c;
 468         dst[0] = (t0 * vr[0] + t1 * vu[0] + vf[0] * vf[0]) * point[0]
 469                + (t0 * vr[1] + t1 * vu[1] + vf[0] * vf[1]) * point[1]
 470                + (t0 * vr[2] + t1 * vu[2] + vf[0] * vf[2]) * point[2];
 471
 472         t0 = vr[1] *  c + vu[1] * -s;
 473         t1 = vr[1] *  s + vu[1] *  c;
 474         dst[1] = (t0 * vr[0] + t1 * vu[0] + vf[1] * vf[0]) * point[0]
 475                + (t0 * vr[1] + t1 * vu[1] + vf[1] * vf[1]) * point[1]
 476                + (t0 * vr[2] + t1 * vu[2] + vf[1] * vf[2]) * point[2];
 477
 478         t0 = vr[2] *  c + vu[2] * -s;
 479         t1 = vr[2] *  s + vu[2] *  c;
 480         dst[2] = (t0 * vr[0] + t1 * vu[0] + vf[2] * vf[0]) * point[0]
 481                + (t0 * vr[1] + t1 * vu[1] + vf[2] * vf[1]) * point[1]
 482                + (t0 * vr[2] + t1 * vu[2] + vf[2] * vf[2]) * point[2];
 483 #endif
 484 }
 485
 486 #ifdef _WIN32
 487 #pragma optimize( "", on )
 488 #endif
 489
 490 /*-----------------------------------------------------------------*/
 491
 492
 493 // LordHavoc note 1:
 494 // BoxOnPlaneSide did a switch on a 'signbits' value and had optimized
 495 // assembly in an attempt to accelerate it further, very inefficient
 496 // considering that signbits of the frustum planes only changed each
 497 // frame, and the world planes changed only at load time.
 498 // So, to optimize it further I took the obvious route of storing a function
 499 // pointer in the plane struct itself, and shrunk each of the individual
 500 // cases to a single return statement.
 501 // LordHavoc note 2:
 502 // realized axial cases would be a nice speedup for world geometry, although
 503 // never useful for the frustum planes.
 504 int BoxOnPlaneSideX (vec3_t emins, vec3_t emaxs, mplane_t *p) {return p->dist <= emins[0] ? 1 : (p->dist >= emaxs[0] ? 2 : 3);}
 505 int BoxOnPlaneSideY (vec3_t emins, vec3_t emaxs, mplane_t *p) {return p->dist <= emins[1] ? 1 : (p->dist >= emaxs[1] ? 2 : 3);}
 506 int BoxOnPlaneSideZ (vec3_t emins, vec3_t emaxs, mplane_t *p) {return p->dist <= emins[2] ? 1 : (p->dist >= emaxs[2] ? 2 : 3);}
 507 int BoxOnPlaneSide0 (vec3_t emins, vec3_t emaxs, mplane_t *p) {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));}
 508 int BoxOnPlaneSide1 (vec3_t emins, vec3_t emaxs, mplane_t *p) {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));}
 509 int BoxOnPlaneSide2 (vec3_t emins, vec3_t emaxs, mplane_t *p) {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));}
 510 int BoxOnPlaneSide3 (vec3_t emins, vec3_t emaxs, mplane_t *p) {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));}
 511 int BoxOnPlaneSide4 (vec3_t emins, vec3_t emaxs, mplane_t *p) {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));}
 512 int BoxOnPlaneSide5 (vec3_t emins, vec3_t emaxs, mplane_t *p) {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));}
 513 int BoxOnPlaneSide6 (vec3_t emins, vec3_t emaxs, mplane_t *p) {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));}
 514 int BoxOnPlaneSide7 (vec3_t emins, vec3_t emaxs, mplane_t *p) {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));}
 515
 516 void BoxOnPlaneSideClassify(mplane_t *p)
 517 {
 518         switch(p->type)
 519         {
 520         case 0: // x axis
 521                 p->BoxOnPlaneSideFunc = BoxOnPlaneSideX;
 522                 break;
 523         case 1: // y axis
 524                 p->BoxOnPlaneSideFunc = BoxOnPlaneSideY;
 525                 break;
 526         case 2: // z axis
 527                 p->BoxOnPlaneSideFunc = BoxOnPlaneSideZ;
 528                 break;
 529         default:
 530                 if (p->normal[2] < 0) // 4
 531                 {
 532                         if (p->normal[1] < 0) // 2
 533                         {
 534                                 if (p->normal[0] < 0) // 1
 535                                         p->BoxOnPlaneSideFunc = BoxOnPlaneSide7;
 536                                 else
 537                                         p->BoxOnPlaneSideFunc = BoxOnPlaneSide6;
 538                         }
 539                         else
 540                         {
 541                                 if (p->normal[0] < 0) // 1
 542                                         p->BoxOnPlaneSideFunc = BoxOnPlaneSide5;
 543                                 else
 544                                         p->BoxOnPlaneSideFunc = BoxOnPlaneSide4;
 545                         }
 546                 }
 547                 else
 548                 {
 549                         if (p->normal[1] < 0) // 2
 550                         {
 551                                 if (p->normal[0] < 0) // 1
 552                                         p->BoxOnPlaneSideFunc = BoxOnPlaneSide3;
 553                                 else
 554                                         p->BoxOnPlaneSideFunc = BoxOnPlaneSide2;
 555                         }
 556                         else
 557                         {
 558                                 if (p->normal[0] < 0) // 1
 559                                         p->BoxOnPlaneSideFunc = BoxOnPlaneSide1;
 560                                 else
 561                                         p->BoxOnPlaneSideFunc = BoxOnPlaneSide0;
 562                         }
 563                 }
 564                 break;
 565         }
 566 }
 567
 568 void PlaneClassify(mplane_t *p)
 569 {
 570         if (p->normal[0] == 1)
 571                 p->type = 0;
 572         else if (p->normal[1] == 1)
 573                 p->type = 1;
 574         else if (p->normal[2] == 1)
 575                 p->type = 2;
 576         else
 577                 p->type = 3;
 578         BoxOnPlaneSideClassify(p);
 579 }
 580
 581 void AngleVectors (vec3_t angles, vec3_t forward, vec3_t right, vec3_t up)
 582 {
 583         float           angle;
 584         float           sr, sp, sy, cr, cp, cy;
 585
 586         angle = angles[YAW] * (M_PI*2 / 360);
 587         sy = sin(angle);
 588         cy = cos(angle);
 589         angle = angles[PITCH] * (M_PI*2 / 360);
 590         sp = sin(angle);
 591         cp = cos(angle);
 592         // LordHavoc: this is only to hush up gcc complaining about 'might be used uninitialized' variables
 593         // (they are NOT used uninitialized, but oh well)
 594         cr = 0;
 595         sr = 0;
 596         if (right || up)
 597         {
 598                 angle = angles[ROLL] * (M_PI*2 / 360);
 599                 sr = sin(angle);
 600                 cr = cos(angle);
 601         }
 602
 603         if (forward)
 604         {
 605                 forward[0] = cp*cy;
 606                 forward[1] = cp*sy;
 607                 forward[2] = -sp;
 608         }
 609         if (right)
 610         {
 611                 right[0] = (-1*sr*sp*cy+-1*cr*-sy);
 612                 right[1] = (-1*sr*sp*sy+-1*cr*cy);
 613                 right[2] = -1*sr*cp;
 614         }
 615         if (up)
 616         {
 617                 up[0] = (cr*sp*cy+-sr*-sy);
 618                 up[1] = (cr*sp*sy+-sr*cy);
 619                 up[2] = cr*cp;
 620         }
 621 }
 622
 623 int VectorCompare (vec3_t v1, vec3_t v2)
 624 {
 625         int             i;
 626
 627         for (i=0 ; i<3 ; i++)
 628                 if (v1[i] != v2[i])
 629                         return 0;
 630
 631         return 1;
 632 }
 633
 634 void VectorMASlow (vec3_t veca, float scale, vec3_t vecb, vec3_t vecc)
 635 {
 636         vecc[0] = veca[0] + scale*vecb[0];
 637         vecc[1] = veca[1] + scale*vecb[1];
 638         vecc[2] = veca[2] + scale*vecb[2];
 639 }
 640
 641
 642 vec_t _DotProduct (vec3_t v1, vec3_t v2)
 643 {
 644         return v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2];
 645 }
 646
 647 void _VectorSubtract (vec3_t veca, vec3_t vecb, vec3_t out)
 648 {
 649         out[0] = veca[0]-vecb[0];
 650         out[1] = veca[1]-vecb[1];
 651         out[2] = veca[2]-vecb[2];
 652 }
 653
 654 void _VectorAdd (vec3_t veca, vec3_t vecb, vec3_t out)
 655 {
 656         out[0] = veca[0]+vecb[0];
 657         out[1] = veca[1]+vecb[1];
 658         out[2] = veca[2]+vecb[2];
 659 }
 660
 661 void _VectorCopy (vec3_t in, vec3_t out)
 662 {
 663         out[0] = in[0];
 664         out[1] = in[1];
 665         out[2] = in[2];
 666 }
 667
 668 // LordHavoc: changed CrossProduct to a #define
 669 /*
 670 void CrossProduct (vec3_t v1, vec3_t v2, vec3_t cross)
 671 {
 672         cross[0] = v1[1]*v2[2] - v1[2]*v2[1];
 673         cross[1] = v1[2]*v2[0] - v1[0]*v2[2];
 674         cross[2] = v1[0]*v2[1] - v1[1]*v2[0];
 675 }
 676 */
 677
 678 double sqrt(double x);
 679
 680 vec_t Length(vec3_t v)
 681 {
 682         int             i;
 683         float   length;
 684
 685         length = 0;
 686         for (i=0 ; i< 3 ; i++)
 687                 length += v[i]*v[i];
 688         length = sqrt (length);         // FIXME
 689
 690         return length;
 691 }
 692
 693 /*
 694 // LordHavoc: fixme: do more research on gcc assembly so that qftol_minushalf and result will not be considered unused
 695 static double qftol_minushalf = -0.5;
 696
 697 int qftol(double v)
 698 {
 699         int result;
 700 #ifdef _MSC_VER
 701         __asm
 702         {
 703                 fld             v
 704                 fadd    qftol_minushalf
 705                 fistp   result
 706         }
 707 #else // gcc hopefully
 708         asm("fldl v\n\tfaddl qftol_minushalf\n\tfistpl result");
 709 #endif
 710         return result;
 711 }
 712 */
 713
 714 // LordHavoc: renamed these to Length, and made the normal ones #define
 715 float VectorNormalizeLength (vec3_t v)
 716 {
 717         float   length, ilength;
 718
 719         length = v[0]*v[0] + v[1]*v[1] + v[2]*v[2];
 720         length = sqrt (length);         // FIXME
 721
 722         if (length)
 723         {
 724                 ilength = 1/length;
 725                 v[0] *= ilength;
 726                 v[1] *= ilength;
 727                 v[2] *= ilength;
 728         }
 729
 730         return length;
 731
 732 }
 733
 734 float VectorNormalizeLength2 (vec3_t v, vec3_t dest) // LordHavoc: added to allow copying while doing the calculation...
 735 {
 736         float   length, ilength;
 737
 738         length = v[0]*v[0] + v[1]*v[1] + v[2]*v[2];
 739         length = sqrt (length);         // FIXME
 740
 741         if (length)
 742         {
 743                 ilength = 1/length;
 744                 dest[0] = v[0] * ilength;
 745                 dest[1] = v[1] * ilength;
 746                 dest[2] = v[2] * ilength;
 747         }
 748         else
 749                 dest[0] = dest[1] = dest[2] = 0;
 750
 751         return length;
 752
 753 }
 754
 755 void _VectorInverse (vec3_t v)
 756 {
 757         v[0] = -v[0];
 758         v[1] = -v[1];
 759         v[2] = -v[2];
 760 }
 761
 762 void _VectorScale (vec3_t in, vec_t scale, vec3_t out)
 763 {
 764         out[0] = in[0]*scale;
 765         out[1] = in[1]*scale;
 766         out[2] = in[2]*scale;
 767 }
 768
 769
 770 int Q_log2(int val)
 771 {
 772         int answer=0;
 773         while (val>>=1)
 774                 answer++;
 775         return answer;
 776 }
 777
 778
 779 /*
 780 ================
 781 R_ConcatRotations
 782 ================
 783 */
 784 void R_ConcatRotations (float in1[3][3], float in2[3][3], float out[3][3])
 785 {
 786         out[0][0] = in1[0][0] * in2[0][0] + in1[0][1] * in2[1][0] + in1[0][2] * in2[2][0];
 787         out[0][1] = in1[0][0] * in2[0][1] + in1[0][1] * in2[1][1] + in1[0][2] * in2[2][1];
 788         out[0][2] = in1[0][0] * in2[0][2] + in1[0][1] * in2[1][2] + in1[0][2] * in2[2][2];
 789         out[1][0] = in1[1][0] * in2[0][0] + in1[1][1] * in2[1][0] + in1[1][2] * in2[2][0];
 790         out[1][1] = in1[1][0] * in2[0][1] + in1[1][1] * in2[1][1] + in1[1][2] * in2[2][1];
 791         out[1][2] = in1[1][0] * in2[0][2] + in1[1][1] * in2[1][2] + in1[1][2] * in2[2][2];
 792         out[2][0] = in1[2][0] * in2[0][0] + in1[2][1] * in2[1][0] + in1[2][2] * in2[2][0];
 793         out[2][1] = in1[2][0] * in2[0][1] + in1[2][1] * in2[1][1] + in1[2][2] * in2[2][1];
 794         out[2][2] = in1[2][0] * in2[0][2] + in1[2][1] * in2[1][2] + in1[2][2] * in2[2][2];
 795 }
 796
 797
 798 /*
 799 ================
 800 R_ConcatTransforms
 801 ================
 802 */
 803 void R_ConcatTransforms (float in1[3][4], float in2[3][4], float out[3][4])
 804 {
 805         out[0][0] = in1[0][0] * in2[0][0] + in1[0][1] * in2[1][0] + in1[0][2] * in2[2][0];
 806         out[0][1] = in1[0][0] * in2[0][1] + in1[0][1] * in2[1][1] + in1[0][2] * in2[2][1];
 807         out[0][2] = in1[0][0] * in2[0][2] + in1[0][1] * in2[1][2] + in1[0][2] * in2[2][2];
 808         out[0][3] = in1[0][0] * in2[0][3] + in1[0][1] * in2[1][3] + in1[0][2] * in2[2][3] + in1[0][3];
 809         out[1][0] = in1[1][0] * in2[0][0] + in1[1][1] * in2[1][0] + in1[1][2] * in2[2][0];
 810         out[1][1] = in1[1][0] * in2[0][1] + in1[1][1] * in2[1][1] + in1[1][2] * in2[2][1];
 811         out[1][2] = in1[1][0] * in2[0][2] + in1[1][1] * in2[1][2] + in1[1][2] * in2[2][2];
 812         out[1][3] = in1[1][0] * in2[0][3] + in1[1][1] * in2[1][3] + in1[1][2] * in2[2][3] + in1[1][3];
 813         out[2][0] = in1[2][0] * in2[0][0] + in1[2][1] * in2[1][0] + in1[2][2] * in2[2][0];
 814         out[2][1] = in1[2][0] * in2[0][1] + in1[2][1] * in2[1][1] + in1[2][2] * in2[2][1];
 815         out[2][2] = in1[2][0] * in2[0][2] + in1[2][1] * in2[1][2] + in1[2][2] * in2[2][2];
 816         out[2][3] = in1[2][0] * in2[0][3] + in1[2][1] * in2[1][3] + in1[2][2] * in2[2][3] + in1[2][3];
 817 }
 818
 819
 820 /*
 821 ===================
 822 FloorDivMod
 823
 824 Returns mathematically correct (floor-based) quotient and remainder for
 825 numer and denom, both of which should contain no fractional part. The
 826 quotient must fit in 32 bits.
 827 ====================
 828 */
 829
 830 void FloorDivMod (double numer, double denom, int *quotient,
 831                 int *rem)
 832 {
 833         int             q, r;
 834         double  x;
 835
 836 #ifndef PARANOID
 837         if (denom <= 0.0)
 838                 Sys_Error ("FloorDivMod: bad denominator %d\n", denom);
 839
 840 //      if ((floor(numer) != numer) || (floor(denom) != denom))
 841 //              Sys_Error ("FloorDivMod: non-integer numer or denom %f %f\n",
 842 //                              numer, denom);
 843 #endif
 844
 845         if (numer >= 0.0)
 846         {
 847
 848                 x = floor(numer / denom);
 849                 q = (int)x;
 850                 r = (int)floor(numer - (x * denom));
 851         }
 852         else
 853         {
 854         //
 855         // perform operations with positive values, and fix mod to make floor-based
 856         //
 857                 x = floor(-numer / denom);
 858                 q = -(int)x;
 859                 r = (int)floor(-numer - (x * denom));
 860                 if (r != 0)
 861                 {
 862                         q--;
 863                         r = (int)denom - r;
 864                 }
 865         }
 866
 867         *quotient = q;
 868         *rem = r;
 869 }
 870
 871
 872 /*
 873 ===================
 874 GreatestCommonDivisor
 875 ====================
 876 */
 877 int GreatestCommonDivisor (int i1, int i2)
 878 {
 879         if (i1 > i2)
 880         {
 881                 if (i2 == 0)
 882                         return (i1);
 883                 return GreatestCommonDivisor (i2, i1 % i2);
 884         }
 885         else
 886         {
 887                 if (i1 == 0)
 888                         return (i2);
 889                 return GreatestCommonDivisor (i1, i2 % i1);
 890         }
 891 }
 892
 893
 894 void Mathlib_Init(void)
 895 {
 896         int a;
 897
 898         // LordHavoc: setup 1.0f / N table for quick recipricols of integers
 899         ixtable[0] = 0;
 900         for (a = 1;a < 4096;a++)
 901                 ixtable[a] = 1.0f / a;
 902 }