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         int i, zeroComponentCount;
  40
  41         zeroComponentCount = 0;
  42         for ( i = 0; i < 3; i++ )
  43         {
  44                 if ( v[i] == 0.0 ) {
  45                         zeroComponentCount++;
  46                 }
  47         }
  48
  49         if ( zeroComponentCount > 1 ) {
  50                 // The zero vector will be on axis.
  51                 return qtrue;
  52         }
  53
  54         return qfalse;
  55 }
  56
  57 /*
  58    ================
  59    VectorIsOnAxialPlane
  60    ================
  61  */
  62 qboolean VectorIsOnAxialPlane( vec3_t v ){
  63         int i;
  64
  65         for ( i = 0; i < 3; i++ )
  66         {
  67                 if ( v[i] == 0.0 ) {
  68                         // The zero vector will be on axial plane.
  69                         return qtrue;
  70                 }
  71         }
  72
  73         return qfalse;
  74 }
  75
  76 /*
  77    ================
  78    MakeNormalVectors
  79
  80    Given a normalized forward vector, create two
  81    other perpendicular vectors
  82    ================
  83  */
  84 void MakeNormalVectors( vec3_t forward, vec3_t right, vec3_t up ){
  85         float d;
  86
  87         // this rotate and negate guarantees a vector
  88         // not colinear with the original
  89         right[1] = -forward[0];
  90         right[2] = forward[1];
  91         right[0] = forward[2];
  92
  93         d = DotProduct( right, forward );
  94         VectorMA( right, -d, forward, right );
  95         VectorNormalize( right, right );
  96         CrossProduct( right, forward, up );
  97 }
  98
  99 vec_t VectorLength( const vec3_t v ){
 100         int i;
 101         float length;
 102
 103         length = 0.0f;
 104         for ( i = 0 ; i < 3 ; i++ )
 105                 length += v[i] * v[i];
 106         length = (float)sqrt( length );
 107
 108         return length;
 109 }
 110
 111 qboolean VectorCompare( const vec3_t v1, const vec3_t v2 ){
 112         int i;
 113
 114         for ( i = 0 ; i < 3 ; i++ )
 115                 if ( fabs( v1[i] - v2[i] ) > EQUAL_EPSILON ) {
 116                         return qfalse;
 117                 }
 118
 119         return qtrue;
 120 }
 121
 122 void VectorMA( const vec3_t va, vec_t scale, const vec3_t vb, vec3_t vc ){
 123         vc[0] = va[0] + scale * vb[0];
 124         vc[1] = va[1] + scale * vb[1];
 125         vc[2] = va[2] + scale * vb[2];
 126 }
 127
 128 void _CrossProduct( vec3_t v1, vec3_t v2, vec3_t cross ){
 129         cross[0] = v1[1] * v2[2] - v1[2] * v2[1];
 130         cross[1] = v1[2] * v2[0] - v1[0] * v2[2];
 131         cross[2] = v1[0] * v2[1] - v1[1] * v2[0];
 132 }
 133
 134 vec_t _DotProduct( vec3_t v1, vec3_t v2 ){
 135         return v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2];
 136 }
 137
 138 void _VectorSubtract( vec3_t va, vec3_t vb, vec3_t out ){
 139         out[0] = va[0] - vb[0];
 140         out[1] = va[1] - vb[1];
 141         out[2] = va[2] - vb[2];
 142 }
 143
 144 void _VectorAdd( vec3_t va, vec3_t vb, vec3_t out ){
 145         out[0] = va[0] + vb[0];
 146         out[1] = va[1] + vb[1];
 147         out[2] = va[2] + vb[2];
 148 }
 149
 150 void _VectorCopy( vec3_t in, vec3_t out ){
 151         out[0] = in[0];
 152         out[1] = in[1];
 153         out[2] = in[2];
 154 }
 155
 156 vec_t VectorNormalize( const vec3_t in, vec3_t out ) {
 157
 158 #if MATHLIB_VECTOR_NORMALIZE_PRECISION_FIX
 159
 160         // The sqrt() function takes double as an input and returns double as an
 161         // output according the the man pages on Debian and on FreeBSD.  Therefore,
 162         // I don't see a reason why using a double outright (instead of using the
 163         // vec_accu_t alias for example) could possibly be frowned upon.
 164
 165         double x, y, z, length;
 166
 167         x = (double) in[0];
 168         y = (double) in[1];
 169         z = (double) in[2];
 170
 171         length = sqrt( ( x * x ) + ( y * y ) + ( z * z ) );
 172         if ( length == 0 ) {
 173                 VectorClear( out );
 174                 return 0;
 175         }
 176
 177         out[0] = (vec_t) ( x / length );
 178         out[1] = (vec_t) ( y / length );
 179         out[2] = (vec_t) ( z / length );
 180
 181         return (vec_t) length;
 182
 183 #else
 184
 185         vec_t length, ilength;
 186
 187         length = (vec_t)sqrt( in[0] * in[0] + in[1] * in[1] + in[2] * in[2] );
 188         if ( length == 0 ) {
 189                 VectorClear( out );
 190                 return 0;
 191         }
 192
 193         ilength = 1.0f / length;
 194         out[0] = in[0] * ilength;
 195         out[1] = in[1] * ilength;
 196         out[2] = in[2] * ilength;
 197
 198         return length;
 199
 200 #endif
 201
 202 }
 203
 204 vec_t ColorNormalize( const vec3_t in, vec3_t out ) {
 205         float max, scale;
 206
 207         max = in[0];
 208         if ( in[1] > max ) {
 209                 max = in[1];
 210         }
 211         if ( in[2] > max ) {
 212                 max = in[2];
 213         }
 214
 215         if ( max == 0 ) {
 216                 out[0] = out[1] = out[2] = 1.0;
 217                 return 0;
 218         }
 219
 220         scale = 1.0f / max;
 221
 222         VectorScale( in, scale, out );
 223
 224         return max;
 225 }
 226
 227 void VectorInverse( vec3_t v ){
 228         v[0] = -v[0];
 229         v[1] = -v[1];
 230         v[2] = -v[2];
 231 }
 232
 233 /*
 234    void VectorScale (vec3_t v, vec_t scale, vec3_t out)
 235    {
 236     out[0] = v[0] * scale;
 237     out[1] = v[1] * scale;
 238     out[2] = v[2] * scale;
 239    }
 240  */
 241
 242 void VectorRotate( vec3_t vIn, vec3_t vRotation, vec3_t out ){
 243         vec3_t vWork, va;
 244         int nIndex[3][2];
 245         int i;
 246
 247         VectorCopy( vIn, va );
 248         VectorCopy( va, vWork );
 249         nIndex[0][0] = 1; nIndex[0][1] = 2;
 250         nIndex[1][0] = 2; nIndex[1][1] = 0;
 251         nIndex[2][0] = 0; nIndex[2][1] = 1;
 252
 253         for ( i = 0; i < 3; i++ )
 254         {
 255                 if ( vRotation[i] != 0 ) {
 256                         float dAngle = vRotation[i] * Q_PI / 180.0f;
 257                         float c = (vec_t)cos( dAngle );
 258                         float s = (vec_t)sin( dAngle );
 259                         vWork[nIndex[i][0]] = va[nIndex[i][0]] * c - va[nIndex[i][1]] * s;
 260                         vWork[nIndex[i][1]] = va[nIndex[i][0]] * s + va[nIndex[i][1]] * c;
 261                 }
 262                 VectorCopy( vWork, va );
 263         }
 264         VectorCopy( vWork, out );
 265 }
 266
 267 void VectorRotateOrigin( vec3_t vIn, vec3_t vRotation, vec3_t vOrigin, vec3_t out ){
 268         vec3_t vTemp, vTemp2;
 269
 270         VectorSubtract( vIn, vOrigin, vTemp );
 271         VectorRotate( vTemp, vRotation, vTemp2 );
 272         VectorAdd( vTemp2, vOrigin, out );
 273 }
 274
 275 void VectorPolar( vec3_t v, float radius, float theta, float phi ){
 276         v[0] = (float)( radius * cos( theta ) * cos( phi ) );
 277         v[1] = (float)( radius * sin( theta ) * cos( phi ) );
 278         v[2] = (float)( radius * sin( phi ) );
 279 }
 280
 281 void VectorSnap( vec3_t v ){
 282         int i;
 283         for ( i = 0; i < 3; i++ )
 284         {
 285                 v[i] = (vec_t)FLOAT_TO_INTEGER( v[i] );
 286         }
 287 }
 288
 289 void VectorISnap( vec3_t point, int snap ){
 290         int i;
 291         for ( i = 0 ; i < 3 ; i++ )
 292         {
 293                 point[i] = (vec_t)FLOAT_SNAP( point[i], snap );
 294         }
 295 }
 296
 297 void VectorFSnap( vec3_t point, float snap ){
 298         int i;
 299         for ( i = 0 ; i < 3 ; i++ )
 300         {
 301                 point[i] = (vec_t)FLOAT_SNAP( point[i], snap );
 302         }
 303 }
 304
 305 void _Vector5Add( vec5_t va, vec5_t vb, vec5_t out ){
 306         out[0] = va[0] + vb[0];
 307         out[1] = va[1] + vb[1];
 308         out[2] = va[2] + vb[2];
 309         out[3] = va[3] + vb[3];
 310         out[4] = va[4] + vb[4];
 311 }
 312
 313 void _Vector5Scale( vec5_t v, vec_t scale, vec5_t out ){
 314         out[0] = v[0] * scale;
 315         out[1] = v[1] * scale;
 316         out[2] = v[2] * scale;
 317         out[3] = v[3] * scale;
 318         out[4] = v[4] * scale;
 319 }
 320
 321 void _Vector53Copy( vec5_t in, vec3_t out ){
 322         out[0] = in[0];
 323         out[1] = in[1];
 324         out[2] = in[2];
 325 }
 326
 327 // NOTE: added these from Ritual's Q3Radiant
 328 #define INVALID_BOUNDS 99999
 329 void ClearBounds( vec3_t mins, vec3_t maxs ){
 330         mins[0] = mins[1] = mins[2] = +INVALID_BOUNDS;
 331         maxs[0] = maxs[1] = maxs[2] = -INVALID_BOUNDS;
 332 }
 333
 334 void AddPointToBounds( vec3_t v, vec3_t mins, vec3_t maxs ){
 335         int i;
 336         vec_t val;
 337
 338         if ( mins[0] == +INVALID_BOUNDS ) {
 339                 if ( maxs[0] == -INVALID_BOUNDS ) {
 340                         VectorCopy( v, mins );
 341                         VectorCopy( v, maxs );
 342                 }
 343         }
 344
 345         for ( i = 0 ; i < 3 ; i++ )
 346         {
 347                 val = v[i];
 348                 if ( val < mins[i] ) {
 349                         mins[i] = val;
 350                 }
 351                 if ( val > maxs[i] ) {
 352                         maxs[i] = val;
 353                 }
 354         }
 355 }
 356
 357 void AngleVectors( vec3_t angles, vec3_t forward, vec3_t right, vec3_t up ){
 358         float angle;
 359         static float sr, sp, sy, cr, cp, cy;
 360         // static to help MS compiler fp bugs
 361
 362         angle = angles[YAW] * ( Q_PI * 2.0f / 360.0f );
 363         sy = (vec_t)sin( angle );
 364         cy = (vec_t)cos( angle );
 365         angle = angles[PITCH] * ( Q_PI * 2.0f / 360.0f );
 366         sp = (vec_t)sin( angle );
 367         cp = (vec_t)cos( angle );
 368         angle = angles[ROLL] * ( Q_PI * 2.0f / 360.0f );
 369         sr = (vec_t)sin( angle );
 370         cr = (vec_t)cos( angle );
 371
 372         if ( forward ) {
 373                 forward[0] = cp * cy;
 374                 forward[1] = cp * sy;
 375                 forward[2] = -sp;
 376         }
 377         if ( right ) {
 378                 right[0] = -sr * sp * cy + cr * sy;
 379                 right[1] = -sr * sp * sy - cr * cy;
 380                 right[2] = -sr * cp;
 381         }
 382         if ( up ) {
 383                 up[0] = cr * sp * cy + sr * sy;
 384                 up[1] = cr * sp * sy - sr * cy;
 385                 up[2] = cr * cp;
 386         }
 387 }
 388
 389 void VectorToAngles( vec3_t vec, vec3_t angles ){
 390         float forward;
 391         float yaw, pitch;
 392
 393         if ( ( vec[ 0 ] == 0 ) && ( vec[ 1 ] == 0 ) ) {
 394                 yaw = 0;
 395                 if ( vec[ 2 ] > 0 ) {
 396                         pitch = 90;
 397                 }
 398                 else
 399                 {
 400                         pitch = 270;
 401                 }
 402         }
 403         else
 404         {
 405                 yaw = (vec_t)atan2( vec[ 1 ], vec[ 0 ] ) * 180 / Q_PI;
 406                 if ( yaw < 0 ) {
 407                         yaw += 360;
 408                 }
 409
 410                 forward = ( float )sqrt( vec[ 0 ] * vec[ 0 ] + vec[ 1 ] * vec[ 1 ] );
 411                 pitch = (vec_t)atan2( vec[ 2 ], forward ) * 180 / Q_PI;
 412                 if ( pitch < 0 ) {
 413                         pitch += 360;
 414                 }
 415         }
 416
 417         angles[ 0 ] = pitch;
 418         angles[ 1 ] = yaw;
 419         angles[ 2 ] = 0;
 420 }
 421
 422 /*
 423    =====================
 424    PlaneFromPoints
 425
 426    Returns false if the triangle is degenrate.
 427    The normal will point out of the clock for clockwise ordered points
 428    =====================
 429  */
 430 qboolean PlaneFromPoints( vec4_t plane, const vec3_t a, const vec3_t b, const vec3_t c ) {
 431         vec3_t d1, d2;
 432
 433         VectorSubtract( b, a, d1 );
 434         VectorSubtract( c, a, d2 );
 435         CrossProduct( d2, d1, plane );
 436         if ( VectorNormalize( plane, plane ) == 0 ) {
 437                 return qfalse;
 438         }
 439
 440         plane[3] = DotProduct( a, plane );
 441         return qtrue;
 442 }
 443
 444 /*
 445 ** NormalToLatLong
 446 **
 447 ** We use two byte encoded normals in some space critical applications.
 448 ** Lat = 0 at (1,0,0) to 360 (-1,0,0), encoded in 8-bit sine table format
 449 ** Lng = 0 at (0,0,1) to 180 (0,0,-1), encoded in 8-bit sine table format
 450 **
 451 */
 452 void NormalToLatLong( const vec3_t normal, byte bytes[2] ) {
 453         // check for singularities
 454         if ( normal[0] == 0 && normal[1] == 0 ) {
 455                 if ( normal[2] > 0 ) {
 456                         bytes[0] = 0;
 457                         bytes[1] = 0;       // lat = 0, long = 0
 458                 }
 459                 else {
 460                         bytes[0] = 128;
 461                         bytes[1] = 0;       // lat = 0, long = 128
 462                 }
 463         }
 464         else {
 465                 int a, b;
 466
 467                 a = (int)( RAD2DEG( atan2( normal[1], normal[0] ) ) * ( 255.0f / 360.0f ) );
 468                 a &= 0xff;
 469
 470                 b = (int)( RAD2DEG( acos( normal[2] ) ) * ( 255.0f / 360.0f ) );
 471                 b &= 0xff;
 472
 473                 bytes[0] = b;   // longitude
 474                 bytes[1] = a;   // lattitude
 475         }
 476 }
 477
 478 /*
 479    =================
 480    PlaneTypeForNormal
 481    =================
 482  */
 483 int PlaneTypeForNormal( vec3_t normal ) {
 484         if ( normal[0] == 1.0 || normal[0] == -1.0 ) {
 485                 return PLANE_X;
 486         }
 487         if ( normal[1] == 1.0 || normal[1] == -1.0 ) {
 488                 return PLANE_Y;
 489         }
 490         if ( normal[2] == 1.0 || normal[2] == -1.0 ) {
 491                 return PLANE_Z;
 492         }
 493
 494         return PLANE_NON_AXIAL;
 495 }
 496
 497 /*
 498    ================
 499    MatrixMultiply
 500    ================
 501  */
 502 void MatrixMultiply( float in1[3][3], float in2[3][3], float out[3][3] ) {
 503         out[0][0] = in1[0][0] * in2[0][0] + in1[0][1] * in2[1][0] +
 504                                 in1[0][2] * in2[2][0];
 505         out[0][1] = in1[0][0] * in2[0][1] + in1[0][1] * in2[1][1] +
 506                                 in1[0][2] * in2[2][1];
 507         out[0][2] = in1[0][0] * in2[0][2] + in1[0][1] * in2[1][2] +
 508                                 in1[0][2] * in2[2][2];
 509         out[1][0] = in1[1][0] * in2[0][0] + in1[1][1] * in2[1][0] +
 510                                 in1[1][2] * in2[2][0];
 511         out[1][1] = in1[1][0] * in2[0][1] + in1[1][1] * in2[1][1] +
 512                                 in1[1][2] * in2[2][1];
 513         out[1][2] = in1[1][0] * in2[0][2] + in1[1][1] * in2[1][2] +
 514                                 in1[1][2] * in2[2][2];
 515         out[2][0] = in1[2][0] * in2[0][0] + in1[2][1] * in2[1][0] +
 516                                 in1[2][2] * in2[2][0];
 517         out[2][1] = in1[2][0] * in2[0][1] + in1[2][1] * in2[1][1] +
 518                                 in1[2][2] * in2[2][1];
 519         out[2][2] = in1[2][0] * in2[0][2] + in1[2][1] * in2[1][2] +
 520                                 in1[2][2] * in2[2][2];
 521 }
 522
 523 void ProjectPointOnPlane( vec3_t dst, const vec3_t p, const vec3_t normal ){
 524         float d;
 525         vec3_t n;
 526         float inv_denom;
 527
 528         inv_denom = 1.0F / DotProduct( normal, normal );
 529
 530         d = DotProduct( normal, p ) * inv_denom;
 531
 532         n[0] = normal[0] * inv_denom;
 533         n[1] = normal[1] * inv_denom;
 534         n[2] = normal[2] * inv_denom;
 535
 536         dst[0] = p[0] - d * n[0];
 537         dst[1] = p[1] - d * n[1];
 538         dst[2] = p[2] - d * n[2];
 539 }
 540
 541 /*
 542 ** assumes "src" is normalized
 543 */
 544 void PerpendicularVector( vec3_t dst, const vec3_t src ){
 545         int pos;
 546         int i;
 547         vec_t minelem = 1.0F;
 548         vec3_t tempvec;
 549
 550         /*
 551         ** find the smallest magnitude axially aligned vector
 552         */
 553         for ( pos = 0, i = 0; i < 3; i++ )
 554         {
 555                 if ( fabs( src[i] ) < minelem ) {
 556                         pos = i;
 557                         minelem = (vec_t)fabs( src[i] );
 558                 }
 559         }
 560         tempvec[0] = tempvec[1] = tempvec[2] = 0.0F;
 561         tempvec[pos] = 1.0F;
 562
 563         /*
 564         ** project the point onto the plane defined by src
 565         */
 566         ProjectPointOnPlane( dst, tempvec, src );
 567
 568         /*
 569         ** normalize the result
 570         */
 571         VectorNormalize( dst, dst );
 572 }
 573
 574 /*
 575    ===============
 576    RotatePointAroundVector
 577
 578    This is not implemented very well...
 579    ===============
 580  */
 581 void RotatePointAroundVector( vec3_t dst, const vec3_t dir, const vec3_t point,
 582                                                           float degrees ) {
 583         float m[3][3];
 584         float im[3][3];
 585         float zrot[3][3];
 586         float tmpmat[3][3];
 587         float rot[3][3];
 588         int i;
 589         vec3_t vr, vup, vf;
 590         float rad;
 591
 592         vf[0] = dir[0];
 593         vf[1] = dir[1];
 594         vf[2] = dir[2];
 595
 596         PerpendicularVector( vr, dir );
 597         CrossProduct( vr, vf, vup );
 598
 599         m[0][0] = vr[0];
 600         m[1][0] = vr[1];
 601         m[2][0] = vr[2];
 602
 603         m[0][1] = vup[0];
 604         m[1][1] = vup[1];
 605         m[2][1] = vup[2];
 606
 607         m[0][2] = vf[0];
 608         m[1][2] = vf[1];
 609         m[2][2] = vf[2];
 610
 611         memcpy( im, m, sizeof( im ) );
 612
 613         im[0][1] = m[1][0];
 614         im[0][2] = m[2][0];
 615         im[1][0] = m[0][1];
 616         im[1][2] = m[2][1];
 617         im[2][0] = m[0][2];
 618         im[2][1] = m[1][2];
 619
 620         memset( zrot, 0, sizeof( zrot ) );
 621         zrot[0][0] = zrot[1][1] = zrot[2][2] = 1.0F;
 622
 623         rad = (float)DEG2RAD( degrees );
 624         zrot[0][0] = (vec_t)cos( rad );
 625         zrot[0][1] = (vec_t)sin( rad );
 626         zrot[1][0] = (vec_t)-sin( rad );
 627         zrot[1][1] = (vec_t)cos( rad );
 628
 629         MatrixMultiply( m, zrot, tmpmat );
 630         MatrixMultiply( tmpmat, im, rot );
 631
 632         for ( i = 0; i < 3; i++ ) {
 633                 dst[i] = rot[i][0] * point[0] + rot[i][1] * point[1] + rot[i][2] * point[2];
 634         }
 635 }
 636
 637
 638 ////////////////////////////////////////////////////////////////////////////////
 639 // Below is double-precision math stuff.  This was initially needed by the new
 640 // "base winding" code in q3map2 brush processing in order to fix the famous
 641 // "disappearing triangles" issue.  These definitions can be used wherever extra
 642 // precision is needed.
 643 ////////////////////////////////////////////////////////////////////////////////
 644
 645 /*
 646    =================
 647    VectorLengthAccu
 648    =================
 649  */
 650 vec_accu_t VectorLengthAccu( const vec3_accu_t v ){
 651         return (vec_accu_t) sqrt( ( v[0] * v[0] ) + ( v[1] * v[1] ) + ( v[2] * v[2] ) );
 652 }
 653
 654 /*
 655    =================
 656    DotProductAccu
 657    =================
 658  */
 659 vec_accu_t DotProductAccu( const vec3_accu_t a, const vec3_accu_t b ){
 660         return ( a[0] * b[0] ) + ( a[1] * b[1] ) + ( a[2] * b[2] );
 661 }
 662
 663 /*
 664    =================
 665    VectorSubtractAccu
 666    =================
 667  */
 668 void VectorSubtractAccu( const vec3_accu_t a, const vec3_accu_t b, vec3_accu_t out ){
 669         out[0] = a[0] - b[0];
 670         out[1] = a[1] - b[1];
 671         out[2] = a[2] - b[2];
 672 }
 673
 674 /*
 675    =================
 676    VectorAddAccu
 677    =================
 678  */
 679 void VectorAddAccu( const vec3_accu_t a, const vec3_accu_t b, vec3_accu_t out ){
 680         out[0] = a[0] + b[0];
 681         out[1] = a[1] + b[1];
 682         out[2] = a[2] + b[2];
 683 }
 684
 685 /*
 686    =================
 687    VectorCopyAccu
 688    =================
 689  */
 690 void VectorCopyAccu( const vec3_accu_t in, vec3_accu_t out ){
 691         out[0] = in[0];
 692         out[1] = in[1];
 693         out[2] = in[2];
 694 }
 695
 696 /*
 697    =================
 698    VectorScaleAccu
 699    =================
 700  */
 701 void VectorScaleAccu( const vec3_accu_t in, vec_accu_t scaleFactor, vec3_accu_t out ){
 702         out[0] = in[0] * scaleFactor;
 703         out[1] = in[1] * scaleFactor;
 704         out[2] = in[2] * scaleFactor;
 705 }
 706
 707 /*
 708    =================
 709    CrossProductAccu
 710    =================
 711  */
 712 void CrossProductAccu( const vec3_accu_t a, const vec3_accu_t b, vec3_accu_t out ){
 713         out[0] = ( a[1] * b[2] ) - ( a[2] * b[1] );
 714         out[1] = ( a[2] * b[0] ) - ( a[0] * b[2] );
 715         out[2] = ( a[0] * b[1] ) - ( a[1] * b[0] );
 716 }
 717
 718 /*
 719    =================
 720    Q_rintAccu
 721    =================
 722  */
 723 vec_accu_t Q_rintAccu( vec_accu_t val ){
 724         return (vec_accu_t) floor( val + 0.5 );
 725 }
 726
 727 /*
 728    =================
 729    VectorCopyAccuToRegular
 730    =================
 731  */
 732 void VectorCopyAccuToRegular( const vec3_accu_t in, vec3_t out ){
 733         out[0] = (vec_t) in[0];
 734         out[1] = (vec_t) in[1];
 735         out[2] = (vec_t) in[2];
 736 }
 737
 738 /*
 739    =================
 740    VectorCopyRegularToAccu
 741    =================
 742  */
 743 void VectorCopyRegularToAccu( const vec3_t in, vec3_accu_t out ){
 744         out[0] = (vec_accu_t) in[0];
 745         out[1] = (vec_accu_t) in[1];
 746         out[2] = (vec_accu_t) in[2];
 747 }
 748
 749 /*
 750    =================
 751    VectorNormalizeAccu
 752    =================
 753  */
 754 vec_accu_t VectorNormalizeAccu( const vec3_accu_t in, vec3_accu_t out ){
 755         // The sqrt() function takes double as an input and returns double as an
 756         // output according the the man pages on Debian and on FreeBSD.  Therefore,
 757         // I don't see a reason why using a double outright (instead of using the
 758         // vec_accu_t alias for example) could possibly be frowned upon.
 759
 760         vec_accu_t length;
 761
 762         length = (vec_accu_t) sqrt( ( in[0] * in[0] ) + ( in[1] * in[1] ) + ( in[2] * in[2] ) );
 763         if ( length == 0 ) {
 764                 VectorClear( out );
 765                 return 0;
 766         }
 767
 768         out[0] = in[0] / length;
 769         out[1] = in[1] / length;
 770         out[2] = in[2] / length;
 771
 772         return length;
 773 }