libs/mathlib/mathlib.c

   1 /*
   2 Copyright (C) 1999-2006 Id Software, Inc. and contributors.
   3 For a list of contributors, see the accompanying CONTRIBUTORS file.
   4
   5 This file is part of GtkRadiant.
   6
   7 GtkRadiant is free software; you can redistribute it and/or modify
   8 it under the terms of the GNU General Public License as published by
   9 the Free Software Foundation; either version 2 of the License, or
  10 (at your option) any later version.
  11
  12 GtkRadiant is distributed in the hope that it will be useful,
  13 but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 GNU General Public License for more details.
  16
  17 You should have received a copy of the GNU General Public License
  18 along with GtkRadiant; if not, write to the Free Software
  19 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
  20 */
  21
  22 // mathlib.c -- math primitives
  23 #include "mathlib.h"
  24 // we use memcpy and memset
  25 #include <memory.h>
  26
  27 const vec3_t vec3_origin = {0.0f,0.0f,0.0f};
  28
  29 const vec3_t g_vec3_axis_x = { 1, 0, 0, };
  30 const vec3_t g_vec3_axis_y = { 0, 1, 0, };
  31 const vec3_t g_vec3_axis_z = { 0, 0, 1, };
  32
  33 /*
  34 ================
  35 VectorIsOnAxis
  36 ================
  37 */
  38 qboolean VectorIsOnAxis(vec3_t v)
  39 {
  40         int     i, zeroComponentCount;
  41
  42         zeroComponentCount = 0;
  43         for (i = 0; i < 3; i++)
  44         {
  45                 if (v[i] == 0.0)
  46                 {
  47                         zeroComponentCount++;
  48                 }
  49         }
  50
  51         if (zeroComponentCount > 1)
  52         {
  53                 // The zero vector will be on axis.
  54                 return qtrue;
  55         }
  56
  57         return qfalse;
  58 }
  59
  60 /*
  61 ================
  62 VectorIsOnAxialPlane
  63 ================
  64 */
  65 qboolean VectorIsOnAxialPlane(vec3_t v)
  66 {
  67         int     i;
  68
  69         for (i = 0; i < 3; i++)
  70         {
  71                 if (v[i] == 0.0)
  72                 {
  73                         // The zero vector will be on axial plane.
  74                         return qtrue;
  75                 }
  76         }
  77
  78         return qfalse;
  79 }
  80
  81 /*
  82 ================
  83 VectorIsOnAxis
  84 ================
  85 */
  86 qboolean VectorIsOnAxis(vec3_t v)
  87 {
  88         int     i, zeroComponentCount;
  89
  90         zeroComponentCount = 0;
  91         for (i = 0; i < 3; i++)
  92         {
  93                 if (v[i] == 0.0)
  94                 {
  95                         zeroComponentCount++;
  96                 }
  97         }
  98
  99         if (zeroComponentCount > 1)
 100         {
 101                 // The zero vector will be on axis.
 102                 return qtrue;
 103         }
 104
 105         return qfalse;
 106 }
 107
 108 /*
 109 ================
 110 VectorIsOnAxialPlane
 111 ================
 112 */
 113 qboolean VectorIsOnAxialPlane(vec3_t v)
 114 {
 115         int     i;
 116
 117         for (i = 0; i < 3; i++)
 118         {
 119                 if (v[i] == 0.0)
 120                 {
 121                         // The zero vector will be on axial plane.
 122                         return qtrue;
 123                 }
 124         }
 125
 126         return qfalse;
 127 }
 128
 129 /*
 130 ================
 131 MakeNormalVectors
 132
 133 Given a normalized forward vector, create two
 134 other perpendicular vectors
 135 ================
 136 */
 137 void MakeNormalVectors (vec3_t forward, vec3_t right, vec3_t up)
 138 {
 139         float           d;
 140
 141         // this rotate and negate guarantees a vector
 142         // not colinear with the original
 143         right[1] = -forward[0];
 144         right[2] = forward[1];
 145         right[0] = forward[2];
 146
 147         d = DotProduct (right, forward);
 148         VectorMA (right, -d, forward, right);
 149         VectorNormalize (right, right);
 150         CrossProduct (right, forward, up);
 151 }
 152
 153 vec_t VectorLength(const vec3_t v)
 154 {
 155         int             i;
 156         float   length;
 157
 158         length = 0.0f;
 159         for (i=0 ; i< 3 ; i++)
 160                 length += v[i]*v[i];
 161         length = (float)sqrt (length);
 162
 163         return length;
 164 }
 165
 166 qboolean VectorCompare (const vec3_t v1, const vec3_t v2)
 167 {
 168         int             i;
 169
 170         for (i=0 ; i<3 ; i++)
 171                 if (fabs(v1[i]-v2[i]) > EQUAL_EPSILON)
 172                         return qfalse;
 173
 174         return qtrue;
 175 }
 176
 177 void VectorMA( const vec3_t va, vec_t scale, const vec3_t vb, vec3_t vc )
 178 {
 179         vc[0] = va[0] + scale*vb[0];
 180         vc[1] = va[1] + scale*vb[1];
 181         vc[2] = va[2] + scale*vb[2];
 182 }
 183
 184 void _CrossProduct (vec3_t v1, vec3_t v2, vec3_t cross)
 185 {
 186         cross[0] = v1[1]*v2[2] - v1[2]*v2[1];
 187         cross[1] = v1[2]*v2[0] - v1[0]*v2[2];
 188         cross[2] = v1[0]*v2[1] - v1[1]*v2[0];
 189 }
 190
 191 vec_t _DotProduct (vec3_t v1, vec3_t v2)
 192 {
 193         return v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2];
 194 }
 195
 196 void _VectorSubtract (vec3_t va, vec3_t vb, vec3_t out)
 197 {
 198         out[0] = va[0]-vb[0];
 199         out[1] = va[1]-vb[1];
 200         out[2] = va[2]-vb[2];
 201 }
 202
 203 void _VectorAdd (vec3_t va, vec3_t vb, vec3_t out)
 204 {
 205         out[0] = va[0]+vb[0];
 206         out[1] = va[1]+vb[1];
 207         out[2] = va[2]+vb[2];
 208 }
 209
 210 void _VectorCopy (vec3_t in, vec3_t out)
 211 {
 212         out[0] = in[0];
 213         out[1] = in[1];
 214         out[2] = in[2];
 215 }
 216
 217 vec_t VectorNormalize( const vec3_t in, vec3_t out ) {
 218
 219 #if MATHLIB_VECTOR_NORMALIZE_PRECISION_FIX
 220
 221         // The sqrt() function takes double as an input and returns double as an
 222         // output according the the man pages on Debian and on FreeBSD.  Therefore,
 223         // I don't see a reason why using a double outright (instead of using the
 224         // vec_accu_t alias for example) could possibly be frowned upon.
 225
 226         double  x, y, z, length;
 227
 228         x = (double) in[0];
 229         y = (double) in[1];
 230         z = (double) in[2];
 231
 232         length = sqrt((x * x) + (y * y) + (z * z));
 233         if (length == 0)
 234         {
 235                 VectorClear (out);
 236                 return 0;
 237         }
 238
 239         out[0] = (vec_t) (x / length);
 240         out[1] = (vec_t) (y / length);
 241         out[2] = (vec_t) (z / length);
 242
 243         return (vec_t) length;
 244
 245 #else
 246
 247         vec_t   length, ilength;
 248
 249         length = (vec_t)sqrt (in[0]*in[0] + in[1]*in[1] + in[2]*in[2]);
 250         if (length == 0)
 251         {
 252                 VectorClear (out);
 253                 return 0;
 254         }
 255
 256         ilength = 1.0f/length;
 257         out[0] = in[0]*ilength;
 258         out[1] = in[1]*ilength;
 259         out[2] = in[2]*ilength;
 260
 261         return length;
 262
 263 #endif
 264
 265 }
 266
 267 vec_t ColorNormalize( const vec3_t in, vec3_t out ) {
 268         float   max, scale;
 269
 270         max = in[0];
 271         if (in[1] > max)
 272                 max = in[1];
 273         if (in[2] > max)
 274                 max = in[2];
 275
 276         if (max == 0) {
 277                 out[0] = out[1] = out[2] = 1.0;
 278                 return 0;
 279         }
 280
 281         scale = 1.0f / max;
 282
 283         VectorScale (in, scale, out);
 284
 285         return max;
 286 }
 287
 288 void VectorInverse (vec3_t v)
 289 {
 290         v[0] = -v[0];
 291         v[1] = -v[1];
 292         v[2] = -v[2];
 293 }
 294
 295 /*
 296 void VectorScale (vec3_t v, vec_t scale, vec3_t out)
 297 {
 298         out[0] = v[0] * scale;
 299         out[1] = v[1] * scale;
 300         out[2] = v[2] * scale;
 301 }
 302 */
 303
 304 void VectorRotate (vec3_t vIn, vec3_t vRotation, vec3_t out)
 305 {
 306   vec3_t vWork, va;
 307   int nIndex[3][2];
 308   int i;
 309
 310   VectorCopy(vIn, va);
 311   VectorCopy(va, vWork);
 312   nIndex[0][0] = 1; nIndex[0][1] = 2;
 313   nIndex[1][0] = 2; nIndex[1][1] = 0;
 314   nIndex[2][0] = 0; nIndex[2][1] = 1;
 315
 316   for (i = 0; i < 3; i++)
 317   {
 318     if (vRotation[i] != 0)
 319     {
 320       float dAngle = vRotation[i] * Q_PI / 180.0f;
 321             float c = (vec_t)cos(dAngle);
 322       float s = (vec_t)sin(dAngle);
 323       vWork[nIndex[i][0]] = va[nIndex[i][0]] * c - va[nIndex[i][1]] * s;
 324       vWork[nIndex[i][1]] = va[nIndex[i][0]] * s + va[nIndex[i][1]] * c;
 325     }
 326     VectorCopy(vWork, va);
 327   }
 328   VectorCopy(vWork, out);
 329 }
 330
 331 void VectorRotateOrigin (vec3_t vIn, vec3_t vRotation, vec3_t vOrigin, vec3_t out)
 332 {
 333   vec3_t vTemp, vTemp2;
 334
 335   VectorSubtract(vIn, vOrigin, vTemp);
 336   VectorRotate(vTemp, vRotation, vTemp2);
 337   VectorAdd(vTemp2, vOrigin, out);
 338 }
 339
 340 void VectorPolar(vec3_t v, float radius, float theta, float phi)
 341 {
 342         v[0]=(float)(radius * cos(theta) * cos(phi));
 343         v[1]=(float)(radius * sin(theta) * cos(phi));
 344         v[2]=(float)(radius * sin(phi));
 345 }
 346
 347 void VectorSnap(vec3_t v)
 348 {
 349   int i;
 350   for (i = 0; i < 3; i++)
 351   {
 352     v[i] = (vec_t)FLOAT_TO_INTEGER(v[i]);
 353   }
 354 }
 355
 356 void VectorISnap(vec3_t point, int snap)
 357 {
 358   int i;
 359         for (i = 0 ;i < 3 ; i++)
 360         {
 361                 point[i] = (vec_t)FLOAT_SNAP(point[i], snap);
 362         }
 363 }
 364
 365 void VectorFSnap(vec3_t point, float snap)
 366 {
 367   int i;
 368         for (i = 0 ;i < 3 ; i++)
 369         {
 370                 point[i] = (vec_t)FLOAT_SNAP(point[i], snap);
 371         }
 372 }
 373
 374 void _Vector5Add (vec5_t va, vec5_t vb, vec5_t out)
 375 {
 376         out[0] = va[0]+vb[0];
 377         out[1] = va[1]+vb[1];
 378         out[2] = va[2]+vb[2];
 379         out[3] = va[3]+vb[3];
 380         out[4] = va[4]+vb[4];
 381 }
 382
 383 void _Vector5Scale (vec5_t v, vec_t scale, vec5_t out)
 384 {
 385         out[0] = v[0] * scale;
 386         out[1] = v[1] * scale;
 387         out[2] = v[2] * scale;
 388         out[3] = v[3] * scale;
 389         out[4] = v[4] * scale;
 390 }
 391
 392 void _Vector53Copy (vec5_t in, vec3_t out)
 393 {
 394         out[0] = in[0];
 395         out[1] = in[1];
 396         out[2] = in[2];
 397 }
 398
 399 // NOTE: added these from Ritual's Q3Radiant
 400 #define INVALID_BOUNDS 99999
 401 void ClearBounds (vec3_t mins, vec3_t maxs)
 402 {
 403         mins[0] = mins[1] = mins[2] = +INVALID_BOUNDS;
 404         maxs[0] = maxs[1] = maxs[2] = -INVALID_BOUNDS;
 405 }
 406
 407 void AddPointToBounds (vec3_t v, vec3_t mins, vec3_t maxs)
 408 {
 409         int             i;
 410         vec_t   val;
 411
 412         if(mins[0] == +INVALID_BOUNDS)
 413         if(maxs[0] == -INVALID_BOUNDS)
 414         {
 415                 VectorCopy(v, mins);
 416                 VectorCopy(v, maxs);
 417         }
 418
 419         for (i=0 ; i<3 ; i++)
 420         {
 421                 val = v[i];
 422                 if (val < mins[i])
 423                         mins[i] = val;
 424                 if (val > maxs[i])
 425                         maxs[i] = val;
 426         }
 427 }
 428
 429 void AngleVectors (vec3_t angles, vec3_t forward, vec3_t right, vec3_t up)
 430 {
 431         float           angle;
 432         static float            sr, sp, sy, cr, cp, cy;
 433         // static to help MS compiler fp bugs
 434
 435         angle = angles[YAW] * (Q_PI*2.0f / 360.0f);
 436         sy = (vec_t)sin(angle);
 437         cy = (vec_t)cos(angle);
 438         angle = angles[PITCH] * (Q_PI*2.0f / 360.0f);
 439         sp = (vec_t)sin(angle);
 440         cp = (vec_t)cos(angle);
 441         angle = angles[ROLL] * (Q_PI*2.0f / 360.0f);
 442         sr = (vec_t)sin(angle);
 443         cr = (vec_t)cos(angle);
 444
 445         if (forward)
 446         {
 447                 forward[0] = cp*cy;
 448                 forward[1] = cp*sy;
 449                 forward[2] = -sp;
 450         }
 451         if (right)
 452         {
 453                 right[0] = -sr*sp*cy+cr*sy;
 454                 right[1] = -sr*sp*sy-cr*cy;
 455                 right[2] = -sr*cp;
 456         }
 457         if (up)
 458         {
 459                 up[0] = cr*sp*cy+sr*sy;
 460                 up[1] = cr*sp*sy-sr*cy;
 461                 up[2] = cr*cp;
 462         }
 463 }
 464
 465 void VectorToAngles( vec3_t vec, vec3_t angles )
 466 {
 467         float forward;
 468         float yaw, pitch;
 469
 470         if ( ( vec[ 0 ] == 0 ) && ( vec[ 1 ] == 0 ) )
 471         {
 472                 yaw = 0;
 473                 if ( vec[ 2 ] > 0 )
 474                 {
 475                         pitch = 90;
 476                 }
 477                 else
 478                 {
 479                         pitch = 270;
 480                 }
 481         }
 482         else
 483         {
 484                 yaw = (vec_t)atan2( vec[ 1 ], vec[ 0 ] ) * 180 / Q_PI;
 485                 if ( yaw < 0 )
 486                 {
 487                         yaw += 360;
 488                 }
 489
 490                 forward = ( float )sqrt( vec[ 0 ] * vec[ 0 ] + vec[ 1 ] * vec[ 1 ] );
 491                 pitch = (vec_t)atan2( vec[ 2 ], forward ) * 180 / Q_PI;
 492                 if ( pitch < 0 )
 493                 {
 494                         pitch += 360;
 495                 }
 496         }
 497
 498         angles[ 0 ] = pitch;
 499         angles[ 1 ] = yaw;
 500         angles[ 2 ] = 0;
 501 }
 502
 503 /*
 504 =====================
 505 PlaneFromPoints
 506
 507 Returns false if the triangle is degenrate.
 508 The normal will point out of the clock for clockwise ordered points
 509 =====================
 510 */
 511 qboolean PlaneFromPoints( vec4_t plane, const vec3_t a, const vec3_t b, const vec3_t c ) {
 512         vec3_t  d1, d2;
 513
 514         VectorSubtract( b, a, d1 );
 515         VectorSubtract( c, a, d2 );
 516         CrossProduct( d2, d1, plane );
 517         if ( VectorNormalize( plane, plane ) == 0 ) {
 518                 return qfalse;
 519         }
 520
 521         plane[3] = DotProduct( a, plane );
 522         return qtrue;
 523 }
 524
 525 /*
 526 ** NormalToLatLong
 527 **
 528 ** We use two byte encoded normals in some space critical applications.
 529 ** Lat = 0 at (1,0,0) to 360 (-1,0,0), encoded in 8-bit sine table format
 530 ** Lng = 0 at (0,0,1) to 180 (0,0,-1), encoded in 8-bit sine table format
 531 **
 532 */
 533 void NormalToLatLong( const vec3_t normal, byte bytes[2] ) {
 534         // check for singularities
 535         if ( normal[0] == 0 && normal[1] == 0 ) {
 536                 if ( normal[2] > 0 ) {
 537                         bytes[0] = 0;
 538                         bytes[1] = 0;           // lat = 0, long = 0
 539                 } else {
 540                         bytes[0] = 128;
 541                         bytes[1] = 0;           // lat = 0, long = 128
 542                 }
 543         } else {
 544                 int     a, b;
 545
 546                 a = (int)( RAD2DEG( atan2( normal[1], normal[0] ) ) * (255.0f / 360.0f ) );
 547                 a &= 0xff;
 548
 549                 b = (int)( RAD2DEG( acos( normal[2] ) ) * ( 255.0f / 360.0f ) );
 550                 b &= 0xff;
 551
 552                 bytes[0] = b;   // longitude
 553                 bytes[1] = a;   // lattitude
 554         }
 555 }
 556
 557 /*
 558 =================
 559 PlaneTypeForNormal
 560 =================
 561 */
 562 int     PlaneTypeForNormal (vec3_t normal) {
 563         if (normal[0] == 1.0 || normal[0] == -1.0)
 564                 return PLANE_X;
 565         if (normal[1] == 1.0 || normal[1] == -1.0)
 566                 return PLANE_Y;
 567         if (normal[2] == 1.0 || normal[2] == -1.0)
 568                 return PLANE_Z;
 569
 570         return PLANE_NON_AXIAL;
 571 }
 572
 573 /*
 574 ================
 575 MatrixMultiply
 576 ================
 577 */
 578 void MatrixMultiply(float in1[3][3], float in2[3][3], float out[3][3]) {
 579         out[0][0] = in1[0][0] * in2[0][0] + in1[0][1] * in2[1][0] +
 580                                 in1[0][2] * in2[2][0];
 581         out[0][1] = in1[0][0] * in2[0][1] + in1[0][1] * in2[1][1] +
 582                                 in1[0][2] * in2[2][1];
 583         out[0][2] = in1[0][0] * in2[0][2] + in1[0][1] * in2[1][2] +
 584                                 in1[0][2] * in2[2][2];
 585         out[1][0] = in1[1][0] * in2[0][0] + in1[1][1] * in2[1][0] +
 586                                 in1[1][2] * in2[2][0];
 587         out[1][1] = in1[1][0] * in2[0][1] + in1[1][1] * in2[1][1] +
 588                                 in1[1][2] * in2[2][1];
 589         out[1][2] = in1[1][0] * in2[0][2] + in1[1][1] * in2[1][2] +
 590                                 in1[1][2] * in2[2][2];
 591         out[2][0] = in1[2][0] * in2[0][0] + in1[2][1] * in2[1][0] +
 592                                 in1[2][2] * in2[2][0];
 593         out[2][1] = in1[2][0] * in2[0][1] + in1[2][1] * in2[1][1] +
 594                                 in1[2][2] * in2[2][1];
 595         out[2][2] = in1[2][0] * in2[0][2] + in1[2][1] * in2[1][2] +
 596                                 in1[2][2] * in2[2][2];
 597 }
 598
 599 void ProjectPointOnPlane( vec3_t dst, const vec3_t p, const vec3_t normal )
 600 {
 601         float d;
 602         vec3_t n;
 603         float inv_denom;
 604
 605         inv_denom = 1.0F / DotProduct( normal, normal );
 606
 607         d = DotProduct( normal, p ) * inv_denom;
 608
 609         n[0] = normal[0] * inv_denom;
 610         n[1] = normal[1] * inv_denom;
 611         n[2] = normal[2] * inv_denom;
 612
 613         dst[0] = p[0] - d * n[0];
 614         dst[1] = p[1] - d * n[1];
 615         dst[2] = p[2] - d * n[2];
 616 }
 617
 618 /*
 619 ** assumes "src" is normalized
 620 */
 621 void PerpendicularVector( vec3_t dst, const vec3_t src )
 622 {
 623         int     pos;
 624         int i;
 625         vec_t minelem = 1.0F;
 626         vec3_t tempvec;
 627
 628         /*
 629         ** find the smallest magnitude axially aligned vector
 630         */
 631         for ( pos = 0, i = 0; i < 3; i++ )
 632         {
 633                 if ( fabs( src[i] ) < minelem )
 634                 {
 635                         pos = i;
 636                         minelem = (vec_t)fabs( src[i] );
 637                 }
 638         }
 639         tempvec[0] = tempvec[1] = tempvec[2] = 0.0F;
 640         tempvec[pos] = 1.0F;
 641
 642         /*
 643         ** project the point onto the plane defined by src
 644         */
 645         ProjectPointOnPlane( dst, tempvec, src );
 646
 647         /*
 648         ** normalize the result
 649         */
 650         VectorNormalize( dst, dst );
 651 }
 652
 653 /*
 654 ===============
 655 RotatePointAroundVector
 656
 657 This is not implemented very well...
 658 ===============
 659 */
 660 void RotatePointAroundVector( vec3_t dst, const vec3_t dir, const vec3_t point,
 661                                                          float degrees ) {
 662         float   m[3][3];
 663         float   im[3][3];
 664         float   zrot[3][3];
 665         float   tmpmat[3][3];
 666         float   rot[3][3];
 667         int     i;
 668         vec3_t vr, vup, vf;
 669         float   rad;
 670
 671         vf[0] = dir[0];
 672         vf[1] = dir[1];
 673         vf[2] = dir[2];
 674
 675         PerpendicularVector( vr, dir );
 676         CrossProduct( vr, vf, vup );
 677
 678         m[0][0] = vr[0];
 679         m[1][0] = vr[1];
 680         m[2][0] = vr[2];
 681
 682         m[0][1] = vup[0];
 683         m[1][1] = vup[1];
 684         m[2][1] = vup[2];
 685
 686         m[0][2] = vf[0];
 687         m[1][2] = vf[1];
 688         m[2][2] = vf[2];
 689
 690         memcpy( im, m, sizeof( im ) );
 691
 692         im[0][1] = m[1][0];
 693         im[0][2] = m[2][0];
 694         im[1][0] = m[0][1];
 695         im[1][2] = m[2][1];
 696         im[2][0] = m[0][2];
 697         im[2][1] = m[1][2];
 698
 699         memset( zrot, 0, sizeof( zrot ) );
 700         zrot[0][0] = zrot[1][1] = zrot[2][2] = 1.0F;
 701
 702         rad = (float)DEG2RAD( degrees );
 703         zrot[0][0] = (vec_t)cos( rad );
 704         zrot[0][1] = (vec_t)sin( rad );
 705         zrot[1][0] = (vec_t)-sin( rad );
 706         zrot[1][1] = (vec_t)cos( rad );
 707
 708         MatrixMultiply( m, zrot, tmpmat );
 709         MatrixMultiply( tmpmat, im, rot );
 710
 711         for ( i = 0; i < 3; i++ ) {
 712                 dst[i] = rot[i][0] * point[0] + rot[i][1] * point[1] + rot[i][2] * point[2];
 713         }
 714 }
 715
 716
 717 ////////////////////////////////////////////////////////////////////////////////
 718 // Below is double-precision math stuff.  This was initially needed by the new
 719 // "base winding" code in q3map2 brush processing in order to fix the famous
 720 // "disappearing triangles" issue.  These definitions can be used wherever extra
 721 // precision is needed.
 722 ////////////////////////////////////////////////////////////////////////////////
 723
 724 /*
 725 =================
 726 VectorLengthAccu
 727 =================
 728 */
 729 vec_accu_t VectorLengthAccu(const vec3_accu_t v)
 730 {
 731         return (vec_accu_t) sqrt((v[0] * v[0]) + (v[1] * v[1]) + (v[2] * v[2]));
 732 }
 733
 734 /*
 735 =================
 736 DotProductAccu
 737 =================
 738 */
 739 vec_accu_t DotProductAccu(const vec3_accu_t a, const vec3_accu_t b)
 740 {
 741         return (a[0] * b[0]) + (a[1] * b[1]) + (a[2] * b[2]);
 742 }
 743
 744 /*
 745 =================
 746 VectorSubtractAccu
 747 =================
 748 */
 749 void VectorSubtractAccu(const vec3_accu_t a, const vec3_accu_t b, vec3_accu_t out)
 750 {
 751         out[0] = a[0] - b[0];
 752         out[1] = a[1] - b[1];
 753         out[2] = a[2] - b[2];
 754 }
 755
 756 /*
 757 =================
 758 VectorAddAccu
 759 =================
 760 */
 761 void VectorAddAccu(const vec3_accu_t a, const vec3_accu_t b, vec3_accu_t out)
 762 {
 763         out[0] = a[0] + b[0];
 764         out[1] = a[1] + b[1];
 765         out[2] = a[2] + b[2];
 766 }
 767
 768 /*
 769 =================
 770 VectorCopyAccu
 771 =================
 772 */
 773 void VectorCopyAccu(const vec3_accu_t in, vec3_accu_t out)
 774 {
 775         out[0] = in[0];
 776         out[1] = in[1];
 777         out[2] = in[2];
 778 }
 779
 780 /*
 781 =================
 782 VectorScaleAccu
 783 =================
 784 */
 785 void VectorScaleAccu(const vec3_accu_t in, vec_accu_t scaleFactor, vec3_accu_t out)
 786 {
 787         out[0] = in[0] * scaleFactor;
 788         out[1] = in[1] * scaleFactor;
 789         out[2] = in[2] * scaleFactor;
 790 }
 791
 792 /*
 793 =================
 794 CrossProductAccu
 795 =================
 796 */
 797 void CrossProductAccu(const vec3_accu_t a, const vec3_accu_t b, vec3_accu_t out)
 798 {
 799         out[0] = (a[1] * b[2]) - (a[2] * b[1]);
 800         out[1] = (a[2] * b[0]) - (a[0] * b[2]);
 801         out[2] = (a[0] * b[1]) - (a[1] * b[0]);
 802 }
 803
 804 /*
 805 =================
 806 Q_rintAccu
 807 =================
 808 */
 809 vec_accu_t Q_rintAccu(vec_accu_t val)
 810 {
 811         return (vec_accu_t) floor(val + 0.5);
 812 }
 813
 814 /*
 815 =================
 816 VectorCopyAccuToRegular
 817 =================
 818 */
 819 void VectorCopyAccuToRegular(const vec3_accu_t in, vec3_t out)
 820 {
 821         out[0] = (vec_t) in[0];
 822         out[1] = (vec_t) in[1];
 823         out[2] = (vec_t) in[2];
 824 }
 825
 826 /*
 827 =================
 828 VectorCopyRegularToAccu
 829 =================
 830 */
 831 void VectorCopyRegularToAccu(const vec3_t in, vec3_accu_t out)
 832 {
 833         out[0] = (vec_accu_t) in[0];
 834         out[1] = (vec_accu_t) in[1];
 835         out[2] = (vec_accu_t) in[2];
 836 }
 837
 838 /*
 839 =================
 840 VectorNormalizeAccu
 841 =================
 842 */
 843 vec_accu_t VectorNormalizeAccu(const vec3_accu_t in, vec3_accu_t out)
 844 {
 845         // The sqrt() function takes double as an input and returns double as an
 846         // output according the the man pages on Debian and on FreeBSD.  Therefore,
 847         // I don't see a reason why using a double outright (instead of using the
 848         // vec_accu_t alias for example) could possibly be frowned upon.
 849
 850         vec_accu_t      length;
 851
 852         length = (vec_accu_t) sqrt((in[0] * in[0]) + (in[1] * in[1]) + (in[2] * in[2]));
 853         if (length == 0)
 854         {
 855                 VectorClear(out);
 856                 return 0;
 857         }
 858
 859         out[0] = in[0] / length;
 860         out[1] = in[1] / length;
 861         out[2] = in[2] / length;
 862
 863         return length;
 864 }
 865
 866
 867 ////////////////////////////////////////////////////////////////////////////////
 868 // Below is double-precision math stuff.  This was initially needed by the new
 869 // "base winding" code in q3map2 brush processing in order to fix the famous
 870 // "disappearing triangles" issue.  These definitions can be used wherever extra
 871 // precision is needed.
 872 ////////////////////////////////////////////////////////////////////////////////
 873
 874 /*
 875 =================
 876 VectorLengthAccu
 877 =================
 878 */
 879 vec_accu_t VectorLengthAccu(const vec3_accu_t v)
 880 {
 881         return (vec_accu_t) sqrt((v[0] * v[0]) + (v[1] * v[1]) + (v[2] * v[2]));
 882 }
 883
 884 /*
 885 =================
 886 DotProductAccu
 887 =================
 888 */
 889 vec_accu_t DotProductAccu(const vec3_accu_t a, const vec3_accu_t b)
 890 {
 891         return (a[0] * b[0]) + (a[1] * b[1]) + (a[2] * b[2]);
 892 }
 893
 894 /*
 895 =================
 896 VectorSubtractAccu
 897 =================
 898 */
 899 void VectorSubtractAccu(const vec3_accu_t a, const vec3_accu_t b, vec3_accu_t out)
 900 {
 901         out[0] = a[0] - b[0];
 902         out[1] = a[1] - b[1];
 903         out[2] = a[2] - b[2];
 904 }
 905
 906 /*
 907 =================
 908 VectorAddAccu
 909 =================
 910 */
 911 void VectorAddAccu(const vec3_accu_t a, const vec3_accu_t b, vec3_accu_t out)
 912 {
 913         out[0] = a[0] + b[0];
 914         out[1] = a[1] + b[1];
 915         out[2] = a[2] + b[2];
 916 }
 917
 918 /*
 919 =================
 920 VectorCopyAccu
 921 =================
 922 */
 923 void VectorCopyAccu(const vec3_accu_t in, vec3_accu_t out)
 924 {
 925         out[0] = in[0];
 926         out[1] = in[1];
 927         out[2] = in[2];
 928 }
 929
 930 /*
 931 =================
 932 VectorScaleAccu
 933 =================
 934 */
 935 void VectorScaleAccu(const vec3_accu_t in, vec_accu_t scaleFactor, vec3_accu_t out)
 936 {
 937         out[0] = in[0] * scaleFactor;
 938         out[1] = in[1] * scaleFactor;
 939         out[2] = in[2] * scaleFactor;
 940 }
 941
 942 /*
 943 =================
 944 CrossProductAccu
 945 =================
 946 */
 947 void CrossProductAccu(const vec3_accu_t a, const vec3_accu_t b, vec3_accu_t out)
 948 {
 949         out[0] = (a[1] * b[2]) - (a[2] * b[1]);
 950         out[1] = (a[2] * b[0]) - (a[0] * b[2]);
 951         out[2] = (a[0] * b[1]) - (a[1] * b[0]);
 952 }
 953
 954 /*
 955 =================
 956 Q_rintAccu
 957 =================
 958 */
 959 vec_accu_t Q_rintAccu(vec_accu_t val)
 960 {
 961         return (vec_accu_t) floor(val + 0.5);
 962 }
 963
 964 /*
 965 =================
 966 VectorCopyAccuToRegular
 967 =================
 968 */
 969 void VectorCopyAccuToRegular(const vec3_accu_t in, vec3_t out)
 970 {
 971         out[0] = (vec_t) in[0];
 972         out[1] = (vec_t) in[1];
 973         out[2] = (vec_t) in[2];
 974 }
 975
 976 /*
 977 =================
 978 VectorCopyRegularToAccu
 979 =================
 980 */
 981 void VectorCopyRegularToAccu(const vec3_t in, vec3_accu_t out)
 982 {
 983         out[0] = (vec_accu_t) in[0];
 984         out[1] = (vec_accu_t) in[1];
 985         out[2] = (vec_accu_t) in[2];
 986 }
 987
 988 /*
 989 =================
 990 VectorNormalizeAccu
 991 =================
 992 */
 993 vec_accu_t VectorNormalizeAccu(const vec3_accu_t in, vec3_accu_t out)
 994 {
 995         // The sqrt() function takes double as an input and returns double as an
 996         // output according the the man pages on Debian and on FreeBSD.  Therefore,
 997         // I don't see a reason why using a double outright (instead of using the
 998         // vec_accu_t alias for example) could possibly be frowned upon.
 999
1000         vec_accu_t      length;
1001
1002         length = (vec_accu_t) sqrt((in[0] * in[0]) + (in[1] * in[1]) + (in[2] * in[2]));
1003         if (length == 0)
1004         {
1005                 VectorClear(out);
1006                 return 0;
1007         }
1008
1009         out[0] = in[0] / length;
1010         out[1] = in[1] / length;
1011         out[2] = in[2] / length;
1012
1013         return length;
1014 }