diff --git a/fix/fix.cpp b/fix/fix.cpp
index 731c26d5..7b9da014 100644
--- a/fix/fix.cpp
+++ b/fix/fix.cpp
@@ -69,24 +69,24 @@
 #include "fix.h"
 
 // Tables for trig functions
-float sincos_table[321]; // 256 entries + 64 sin-only + 1 for interpolation
+scalar sincos_table[321]; // 256 entries + 64 sin-only + 1 for interpolation
 angle asin_table[257];   // 1 quadrants worth, +1 for interpolation
 angle acos_table[257];
 
 // Generate the data for the trig tables
 void InitMathTables() {
   int i;
-  float rad, s, c;
+  scalar rad, s, c;
 
   for (i = 0; i < 321; i++) {
-    rad = (float)((double)i / 256.0 * 2 * PI);
+    rad = (scalar)((double)i / 256.0 * 2 * PI);
     sincos_table[i] = std::sin(rad);
   }
 
   for (i = 0; i < 256; i++) {
 
-    s = std::asin((float)i / 256.0f);
-    c = std::acos((float)i / 256.0f);
+    s = std::asin((scalar)i / 256.0f);
+    c = std::acos((scalar)i / 256.0f);
 
     s = (s / (PI * 2));
     c = (c / (PI * 2));
@@ -100,36 +100,36 @@ void InitMathTables() {
 }
 
 // Returns the sine of the given angle.  Linearly interpolates between two entries in a 256-entry table
-float FixSin(angle a) {
+scalar FixSin(angle a) {
   int i, f;
-  float s0, s1;
+  scalar s0, s1;
 
   i = (a >> 8) & 0xff;
   f = a & 0xff;
 
   s0 = sincos_table[i];
   s1 = sincos_table[i + 1];
-  return (float)(s0 + ((s1 - s0) * (double)f / 256.0));
+  return (scalar)(s0 + ((s1 - s0) * (double)f / 256.0));
 }
 
 // Returns the cosine of the given angle.  Linearly interpolates between two entries in a 256-entry table
-float FixCos(angle a) {
+scalar FixCos(angle a) {
   int i, f;
-  float c0, c1;
+  scalar c0, c1;
 
   i = (a >> 8) & 0xff;
   f = a & 0xff;
 
   c0 = sincos_table[i + 64];
   c1 = sincos_table[i + 64 + 1];
-  return (float)(c0 + ((c1 - c0) * (double)f / 256.0));
+  return (scalar)(c0 + ((c1 - c0) * (double)f / 256.0));
 }
 
 // Get rid of the "no return value" warnings in the next three functions
 #pragma warning(disable : 4035)
 
 // compute inverse sine
-angle FixAsin(float v) {
+angle FixAsin(scalar v) {
   fix vv;
   int i, f, aa;
 
@@ -151,7 +151,7 @@ angle FixAsin(float v) {
 }
 
 // compute inverse cosine
-angle FixAcos(float v) {
+angle FixAcos(scalar v) {
   fix vv;
   int i, f, aa;
 
@@ -177,8 +177,8 @@ angle FixAcos(float v) {
 // equal the ratio of the actual cos & sin for the result angle, but the parms
 // need not be the actual cos & sin.
 // NOTE: this is different from the standard C atan2, since it is left-handed.
-angle FixAtan2(float cos, float sin) {
-  float q, m;
+angle FixAtan2(scalar cos, scalar sin) {
+  scalar q, m;
   angle t;
 
   // find smaller of two
diff --git a/fix/fix.h b/fix/fix.h
index 62229fcf..c0164c4a 100644
--- a/fix/fix.h
+++ b/fix/fix.h
@@ -63,12 +63,17 @@
 
 #include <cmath>
 #include <cstdint>
+#include <type_traits>
+
+typedef std::make_signed<size_t>::type ssize_t;
+// the basic floating-point type
+using scalar = float;
 
 // Angles are unsigned shorts
-typedef uint16_t angle;
+using angle = uint16_t;
 
 // The basic fixed-point type
-typedef int32_t fix;
+using fix = int32_t;
 
 #define PI 3.141592654f
 #define PIOVER2 1.570796327 // DAJ
@@ -84,10 +89,10 @@ typedef int32_t fix;
 void InitMathTables();
 
 // Returns the sine of the given angle.  Linearly interpolates between two entries in a 256-entry table
-float FixSin(angle a);
+scalar FixSin(angle a);
 
 // Returns the cosine of the given angle.  Linearly interpolates between two entries in a 256-entry table
-float FixCos(angle a);
+scalar FixCos(angle a);
 
 #define Round(x) ((int)(x + 0.5))
 
@@ -96,11 +101,11 @@ float FixCos(angle a);
 #define FloatToFix(num) ((fix)((num) * FLOAT_SCALER))
 #define IntToFix(num) ((num) << FIX_SHIFT)
 #define ShortToFix(num) (((int32_t)(num)) << FIX_SHIFT)
-#define FixToFloat(num) (((float)(num)) / FLOAT_SCALER)
+#define FixToFloat(num) (((scalar)(num)) / FLOAT_SCALER)
 #define FixToInt(num) ((num) >> FIX_SHIFT)
 
-angle FixAtan2(float cos, float sin);
-angle FixAsin(float v);
-angle FixAcos(float v);
+angle FixAtan2(scalar cos, scalar sin);
+angle FixAsin(scalar v);
+angle FixAcos(scalar v);
 
 #endif
diff --git a/lib/vecmat.h b/lib/vecmat.h
index c406b894..e2dffdc4 100644
--- a/lib/vecmat.h
+++ b/lib/vecmat.h
@@ -153,8 +153,6 @@
 
 #include <cmath>
 
-#include "fix.h"
-
 // All structs, defines and inline functions are located in vecmat_external.h
 // vecmat_external.h is where anything that can be used by DLLs should be.
 #include "vecmat_external.h"
diff --git a/lib/vecmat_external.h b/lib/vecmat_external.h
index 474cb638..8b620c87 100644
--- a/lib/vecmat_external.h
+++ b/lib/vecmat_external.h
@@ -40,48 +40,176 @@
  * $NoKeywords: $
  */
 
-#ifndef VECMAT_EXTERNAL_H
-#define VECMAT_EXTERNAL_H
+#pragma once
 
-#include <cstdint>
-
-// Angles are unsigned shorts
-typedef uint16_t angle; // make sure this matches up with fix.h
+#include "fix.h"
 
 struct angvec {
-  angle p, h, b;
+using T = angle;
+constexpr static const size_t N       = 3;
+constexpr static const size_t PITCH   = 0;
+constexpr static const size_t HEADING = 1;
+constexpr static const size_t BANK    = 2;
+union {
+  T phb[N];
+  struct { T p, h, b; };
+};
+constexpr static inline angvec id(ssize_t i = -1)
+{
+  return angvec{
+  ((i % N) == PITCH)   ? (T)1 : (T)0,
+  ((i % N) == HEADING) ? (T)1 : (T)0,
+  ((i % N) == BANK)    ? (T)1 : (T)0
+  };
+}
+constexpr static inline angvec ne()
+{
+  return angvec{ (T)0, (T)0, (T)0 };
+}
 };
-
-#define IDENTITY_MATRIX                                                                                                \
-  {                                                                                                                    \
-    {1.0, 0, 0}, {0, 1.0, 0}, { 0, 0, 1.0 }                                                                            \
-  }
 
 struct vector {
-  float x, y, z;
+using T = scalar;
+constexpr static const size_t N = 3;
+constexpr static const size_t X = 0;
+constexpr static const size_t Y = 1;
+constexpr static const size_t Z = 2;
+union {
+  T xyz[N];
+  struct { T x, y, z; };
+  struct { T r, g, b; };
+};
+constexpr static inline const vector id(ssize_t i = -1)
+{
+  return vector{
+  ((i % N) == X) ? (T)1 : (T)0,
+  ((i % N) == Y) ? (T)1 : (T)0,
+  ((i % N) == Z) ? (T)1 : (T)0
+  };
+}
+constexpr static inline const vector ne()
+{
+  return vector{ (T)0, (T)0, (T)0 };
+}
 };
 
-struct vector4 {
-  float x, y, z, kat_pad;
+using vector_array = vector;
+
+struct alignas(sizeof(scalar) * 4) aligned_vector {
+using T = scalar;
+constexpr static const size_t N = 3;
+constexpr static const size_t X = 0;
+constexpr static const size_t Y = 1;
+constexpr static const size_t Z = 2;
+union {
+  T xyz[N];
+  struct { T x, y, z; };
+  struct { T r, g, b; };
+};
+constexpr static inline const aligned_vector id(ssize_t i = -1)
+{
+  return aligned_vector{
+  ((i % N) == X) ? (T)1 : (T)0,
+  ((i % N) == Y) ? (T)1 : (T)0,
+  ((i % N) == Z) ? (T)1 : (T)0,
+  };
+}
+constexpr static inline const aligned_vector ne()
+{
+  return aligned_vector{ (scalar)0, (scalar)0, (scalar)0 };
+}
 };
 
-struct vector_array {
-  float xyz[3];
+using aligned_vector_array = aligned_vector;
+
+struct alignas(sizeof(scalar) * 4) vector4 {
+using T = scalar;
+constexpr static const size_t N = 4;
+constexpr static const size_t X = 0;
+constexpr static const size_t Y = 1;
+constexpr static const size_t Z = 2;
+constexpr static const size_t W = 3;
+union {
+  T xyzw[N];
+  struct { T x, y, z; union { T w; T kat_pad; }; };
+  struct { union { T r; T l; }; T g, b, a; };
+  struct { T u,v; union { T u2; T s; }; union { T v2; T t; }; };
 };
 
+constexpr static inline const vector4 id(ssize_t i = -1)
+{
+  return vector4{
+  ((i % N) == X) ? (T)1 : (T)0,
+  ((i % N) == Y) ? (T)1 : (T)0,
+  ((i % N) == Z) ? (T)1 : (T)0,
+  ((i % N) == W) ? (T)1 : (T)0
+  };
+}
+constexpr static inline const vector4 ne()
+{
+  return vector4{ (T)0, (T)0, (T)0, (T)0 };
+}
+};
+
+using vector4_array = vector4;
+using aligned_vector4 = vector4;
+using aligned_vector4_array = aligned_vector4;
+
 struct matrix {
-  vector rvec, uvec, fvec;
+constexpr static const size_t RIGHT_HAND = 0;
+constexpr static const size_t UP         = 1;
+constexpr static const size_t FORWARD    = 2;
+union {
+struct {
+  vector rvec;
+  vector uvec;
+  vector fvec;
+};
+  scalar a2d[3][3];
+  scalar a1d[9];
+};
+constexpr static inline const matrix id()
+{
+  return matrix{ vector::id(0), vector::id(1), vector::id(2) };
+}
+constexpr static inline const matrix ne()
+{
+  return matrix{ vector::ne(), vector::ne(), vector::ne() };
+}
 };
 
-struct matrix4 {
-  vector4 rvec, uvec, fvec;
+constexpr matrix IDENTITY_MATRIX = matrix::id();
+
+struct alignas(alignof(vector4) * 4) matrix4 {
+constexpr static const size_t RIGHT_HAND = 0;
+constexpr static const size_t UP         = 1;
+constexpr static const size_t FORWARD    = 2;
+constexpr static const size_t POSITION   = 3;
+union {
+struct {
+  vector4 rvec;
+  vector4 uvec;
+  vector4 fvec;
+  vector4 pos;
+};
+  scalar a2d[4][4];
+  scalar a1d[16];
+};
+constexpr static inline const matrix4 id()
+{
+  return matrix4{ vector4::id(0), vector4::id(1), vector4::id(2), vector4::id(3) };
+}
+constexpr static inline const matrix4 ne()
+{
+  return matrix4{ vector4::ne(), vector4::ne(), vector4::ne(), vector4::ne() };
+}
 };
 
 // Zero's out a vector
-static inline void vm_MakeZero(vector *v) { v->x = v->y = v->z = 0; }
+static inline void vm_MakeZero(vector *v) { *v = vector::ne(); }
 
 // Set an angvec to {0,0,0}
-static inline void vm_MakeZero(angvec *a) { a->p = a->h = a->b = 0; }
+static inline void vm_MakeZero(angvec *a) { *a = angvec::ne(); }
 
 // Checks for equality
 static inline bool operator==(vector a, vector b) {
@@ -271,4 +399,3 @@ static inline float vm_Dot3Vector(float x, float y, float z, vector *v) { return
 
 #define vm_GetSurfaceNormal vm_GetNormal
 
-#endif
diff --git a/vecmat/vector.cpp b/vecmat/vector.cpp
index 6ae13b24..88a54402 100644
--- a/vecmat/vector.cpp
+++ b/vecmat/vector.cpp
@@ -154,15 +154,15 @@
 #include "pserror.h"
 #include "psrand.h"
 
-const vector Zero_vector = {0.0, 0.0, 0.0};
+const vector Zero_vector = vector::ne();
 const matrix Identity_matrix = IDENTITY_MATRIX;
 
 void vm_AverageVector(vector *a, int num) {
   // Averages a vector.  ie divides each component of vector a by num
   ASSERT(num != 0);
-  a->x = a->x / (float)num;
-  a->y = a->y / (float)num;
-  a->z = a->z / (float)num;
+  a->x = a->x / (scalar)num;
+  a->y = a->y / (scalar)num;
+  a->z = a->z / (scalar)num;
 }
 
 void vm_AddVectors(vector *result, vector *a, vector *b) {
@@ -181,21 +181,21 @@ void vm_SubVectors(vector *result, const vector *a, const vector *b) {
   result->z = a->z - b->z;
 }
 
-float vm_VectorDistance(const vector *a, const vector *b) {
+scalar vm_VectorDistance(const vector *a, const vector *b) {
   // Given two vectors, returns the distance between them
 
   vector dest;
-  float dist;
+  scalar dist;
 
   vm_SubVectors(&dest, a, b);
   dist = vm_GetMagnitude(&dest);
   return dist;
 }
-float vm_VectorDistanceQuick(vector *a, vector *b) {
+scalar vm_VectorDistanceQuick(vector *a, vector *b) {
   // Given two vectors, returns the distance between them
 
   vector dest;
-  float dist;
+  scalar dist;
 
   vm_SubVectors(&dest, a, b);
   dist = vm_GetMagnitudeFast(&dest);
@@ -222,24 +222,24 @@ void vm_GetPerp(vector *n, vector *a, vector *b, vector *c) {
 //			v0,v1,v2 - three clockwise vertices
 // Returns the magnitude of the normal before it was normalized.
 // The bigger this value, the better the normal.
-float vm_GetNormal(vector *n, vector *v0, vector *v1, vector *v2) {
+scalar vm_GetNormal(vector *n, vector *v0, vector *v1, vector *v2) {
   vm_GetPerp(n, v0, v1, v2);
 
   return vm_NormalizeVector(n);
 }
 
 // Does a simple dot product calculation
-float vm_DotProduct(const vector *u, const vector *v) { return (u->x * v->x) + (u->y * v->y) + (u->z * v->z); }
+scalar vm_DotProduct(const vector *u, const vector *v) { return (u->x * v->x) + (u->y * v->y) + (u->z * v->z); }
 
 // Scales all components of vector v by value s and stores result in vector d
 // dest can equal source
-void vm_ScaleVector(vector *d, vector *v, float s) {
+void vm_ScaleVector(vector *d, vector *v, scalar s) {
   d->x = (v->x * s);
   d->y = (v->y * s);
   d->z = (v->z * s);
 }
 
-void vm_ScaleAddVector(vector *d, vector *p, vector *v, float s) {
+void vm_ScaleAddVector(vector *d, vector *p, vector *v, scalar s) {
   // Scales all components of vector v by value s
   // adds the result to p and stores result in vector d
   // dest can equal source
@@ -249,7 +249,7 @@ void vm_ScaleAddVector(vector *d, vector *p, vector *v, float s) {
   d->z = p->z + (v->z * s);
 }
 
-void vm_DivVector(vector *dest, vector *src, float n) {
+void vm_DivVector(vector *dest, vector *src, scalar n) {
   // Divides a vector into n portions
   // Dest can equal src
 
@@ -270,8 +270,8 @@ void vm_CrossProduct(vector *dest, vector *u, vector *v) {
 
 // Normalize a vector.
 // Returns:  the magnitude before normalization
-float vm_NormalizeVector(vector *a) {
-  float mag;
+scalar vm_NormalizeVector(vector *a) {
+  scalar mag;
 
   mag = vm_GetMagnitude(a);
 
@@ -286,8 +286,8 @@ float vm_NormalizeVector(vector *a) {
   return mag;
 }
 
-float vm_GetMagnitude(vector *a) {
-  float f;
+scalar vm_GetMagnitude(vector *a) {
+  scalar f;
 
   f = (a->x * a->x) + (a->y * a->y) + (a->z * a->z);
 
@@ -308,7 +308,7 @@ void vm_MakeInverseMatrix(matrix *dest) {
 void vm_TransposeMatrix(matrix *m) {
   // Transposes a matrix in place
 
-  float t;
+  scalar t;
 
   t = m->uvec.x;
   m->uvec.x = m->rvec.y;
@@ -402,7 +402,7 @@ matrix operator*=(matrix &src0, matrix src1) { return (src0 = src0 * src1); }
 // Parameters:	dest - filled in with the normalized direction vector
 //					start,end - the start and end points used to calculate the vector
 // Returns:		the distance between the two input points
-float vm_GetNormalizedDir(vector *dest, vector *end, vector *start) {
+scalar vm_GetNormalizedDir(vector *dest, vector *end, vector *start) {
   vm_SubVectors(dest, end, start);
   return vm_NormalizeVector(dest);
 }
@@ -412,31 +412,31 @@ float vm_GetNormalizedDir(vector *dest, vector *end, vector *start) {
 // Parameters:	dest - filled in with the normalized direction vector
 //					start,end - the start and end points used to calculate the vector
 // Returns:		the distance between the two input points
-float vm_GetNormalizedDirFast(vector *dest, vector *end, vector *start) {
+scalar vm_GetNormalizedDirFast(vector *dest, vector *end, vector *start) {
   vm_SubVectors(dest, end, start);
   return vm_NormalizeVectorFast(dest);
 }
 
-float vm_GetMagnitudeFast(vector *v) {
-  float a, b, c, bc;
+scalar vm_GetMagnitudeFast(vector *v) {
+  scalar a, b, c, bc;
 
   a = fabs(v->x);
   b = fabs(v->y);
   c = fabs(v->z);
 
   if (a < b) {
-    float t = a;
+    scalar t = a;
     a = b;
     b = t;
   }
 
   if (b < c) {
-    float t = b;
+    scalar t = b;
     b = c;
     c = t;
 
     if (a < b) {
-      float t = a;
+      scalar t = a;
       a = b;
       b = t;
     }
@@ -449,8 +449,8 @@ float vm_GetMagnitudeFast(vector *v) {
 
 // Normalize a vector using an approximation of the magnitude
 // Returns:  the magnitude before normalization
-float vm_NormalizeVectorFast(vector *a) {
-  float mag;
+scalar vm_NormalizeVectorFast(vector *a) {
+  scalar mag;
 
   mag = vm_GetMagnitudeFast(a);
 
@@ -471,7 +471,7 @@ float vm_NormalizeVectorFast(vector *a) {
 // Parms:	norm - the (normalized) surface normal of the plane
 //				planep - a point on the plane
 // Returns:	The signed distance from the plane; negative dist is on the back of the plane
-float vm_DistToPlane(vector *checkp, vector *norm, vector *planep) {
+scalar vm_DistToPlane(vector *checkp, vector *norm, vector *planep) {
   vector t;
 
   t = *checkp - *planep;
@@ -479,9 +479,9 @@ float vm_DistToPlane(vector *checkp, vector *norm, vector *planep) {
   return t * *norm;
 }
 
-float vm_GetSlope(float x1, float y1, float x2, float y2) {
+scalar vm_GetSlope(scalar x1, scalar y1, scalar x2, scalar y2) {
   // returns the slope of a line
-  float r;
+  scalar r;
 
   if (y2 - y1 == 0)
     return (0.0);
@@ -490,8 +490,8 @@ float vm_GetSlope(float x1, float y1, float x2, float y2) {
   return (r);
 }
 
-void vm_SinCosToMatrix(matrix *m, float sinp, float cosp, float sinb, float cosb, float sinh, float cosh) {
-  float sbsh, cbch, cbsh, sbch;
+void vm_SinCosToMatrix(matrix *m, scalar sinp, scalar cosp, scalar sinb, scalar cosb, scalar sinh, scalar cosh) {
+  scalar sbsh, cbch, cbsh, sbch;
 
   sbsh = (sinb * sinh);
   cbch = (cosb * cosh);
@@ -513,7 +513,7 @@ void vm_SinCosToMatrix(matrix *m, float sinp, float cosp, float sinb, float cosb
 }
 
 void vm_AnglesToMatrix(matrix *m, angle p, angle h, angle b) {
-  float sinp, cosp, sinb, cosb, sinh, cosh;
+  scalar sinp, cosp, sinb, cosb, sinh, cosh;
 
   sinp = FixSin(p);
   cosp = FixCos(p);
@@ -530,7 +530,7 @@ void vm_AnglesToMatrix(matrix *m, angle p, angle h, angle b) {
 //					v - the forward vector of the new matrix
 //					a - the angle of rotation around the forward vector
 void vm_VectorAngleToMatrix(matrix *m, vector *v, angle a) {
-  float sinb, cosb, sinp, cosp, sinh, cosh;
+  scalar sinb, cosb, sinp, cosp, sinh, cosh;
 
   sinb = FixSin(a);
   cosb = FixCos(a);
@@ -665,7 +665,7 @@ void vm_VectorToMatrix(matrix *m, vector *fvec, vector *uvec, vector *rvec) {
   }
 }
 
-void vm_SinCos(uint16_t a, float *s, float *c) {
+void vm_SinCos(uint16_t a, scalar *s, scalar *c) {
   if (s)
     *s = FixSin(a);
   if (c)
@@ -677,7 +677,7 @@ void vm_SinCos(uint16_t a, float *s, float *c) {
 
 // extract angles from a matrix
 angvec *vm_ExtractAnglesFromMatrix(angvec *a, matrix *m) {
-  float sinh, cosh, cosp, sinb, cosb;
+  scalar sinh, cosh, cosp, sinb, cosb;
 
   // Deal with straight up or straight down
   if (IS_ZERO(m->fvec.x) && IS_ZERO(m->fvec.z)) {
@@ -706,7 +706,7 @@ angvec *vm_ExtractAnglesFromMatrix(angvec *a, matrix *m) {
 }
 
 // returns the value of a determinant
-float calc_det_value(matrix *det) {
+scalar calc_det_value(matrix *det) {
   return det->rvec.x * det->uvec.y * det->fvec.z - det->rvec.x * det->uvec.z * det->fvec.y -
          det->rvec.y * det->uvec.x * det->fvec.z + det->rvec.y * det->uvec.z * det->fvec.x +
          det->rvec.z * det->uvec.x * det->fvec.y - det->rvec.z * det->uvec.y * det->fvec.x;
@@ -749,10 +749,10 @@ angle vm_DeltaAngVecNorm(vector *v0, vector *v1, vector *fvec) {
 
 // Gets the real center of a polygon
 // Returns the size of the passed in stuff
-float vm_GetCentroid(vector *centroid, vector *src, int nv) {
+scalar vm_GetCentroid(vector *centroid, vector *src, int nv) {
   ASSERT(nv > 2);
   vector normal;
-  float area, total_area;
+  scalar area, total_area;
   int i;
   vector tmp_center;
 
@@ -803,10 +803,10 @@ float vm_GetCentroid(vector *centroid, vector *src, int nv) {
 
 // Gets the real center of a polygon, but uses fast magnitude calculation
 // Returns the size of the passed in stuff
-float vm_GetCentroidFast(vector *centroid, vector *src, int nv) {
+scalar vm_GetCentroidFast(vector *centroid, vector *src, int nv) {
   ASSERT(nv > 2);
   vector normal;
-  float area, total_area;
+  scalar area, total_area;
   int i;
   vector tmp_center;
 
@@ -863,13 +863,13 @@ void vm_MakeRandomVector(vector *vec) {
 }
 
 // Given a set of points, computes the minimum bounding sphere of those points
-float vm_ComputeBoundingSphere(vector *center, vector *vecs, int num_verts) {
+scalar vm_ComputeBoundingSphere(vector *center, vector *vecs, int num_verts) {
   // This algorithm is from Graphics Gems I.  There's a better algorithm in Graphics Gems III that
   // we should probably implement sometime.
 
   vector *min_x, *max_x, *min_y, *max_y, *min_z, *max_z, *vp;
-  float dx, dy, dz;
-  float rad, rad2;
+  scalar dx, dy, dz;
+  scalar rad, rad2;
   int i;
 
   // Initialize min, max vars
@@ -923,14 +923,14 @@ float vm_ComputeBoundingSphere(vector *center, vector *vecs, int num_verts) {
   rad2 = rad * rad;
   for (i = 0, vp = vecs; i < num_verts; i++, vp++) {
     vector delta;
-    float t2;
+    scalar t2;
 
     delta = *vp - *center;
     t2 = delta.x * delta.x + delta.y * delta.y + delta.z * delta.z;
 
     // If point outside, make the sphere bigger
     if (t2 > rad2) {
-      float t;
+      scalar t;
 
       t = sqrt(t2);
       rad = (rad + t) / 2;