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