diff --git a/detail/coord.hpp b/detail/coord.hpp
index 302c9585..82c052ed 100644
--- a/detail/coord.hpp
+++ b/detail/coord.hpp
@@ -34,111 +34,111 @@ namespace util {
 #if defined(COMPILER_CLANG)
     #pragma GCC diagnostic ignored "-Wgnu"
 #endif
-        template <size_t S>
-        struct coord_data {
-            coord_data () { ; }
+        template <size_t S, typename T>
+        struct coord {
+            coord () { ; }
 
-            template <typename ...T>
-            coord_data (T ..._t): data{_t...}
+            template <typename ...U>
+            coord (U ..._u): data{_u...}
             { ; }
 
-            double data[S];
+            T data[S];
 
-            double& operator[] (size_t i)       { return data[i]; }
-            double  operator[] (size_t i) const { return data[i]; }
+            T& operator[] (size_t i)       { return data[i]; }
+            T  operator[] (size_t i) const { return data[i]; }
         };
 
-        template <>
-        struct coord_data<1> {
-            coord_data () { ; }
+        template <typename T>
+        struct coord<1,T> {
+            coord () { ; }
 
-            template <typename ...T>
-            coord_data (T ..._t): data{_t...}
+            template <typename ...U>
+            coord (U ..._u): data{_u...}
             { ; }
 
             union {
-                double data[1];
-                double x;
+                T data[1];
+                T x;
             };
 
-            double& operator[] (size_t i)       { return data[i]; }
-            double  operator[] (size_t i) const { return data[i]; }
+            T& operator[] (size_t i)       { return data[i]; }
+            T  operator[] (size_t i) const { return data[i]; }
         };
 
-        template <>
-        struct coord_data<2> {
-            coord_data () { ; }
+        template <typename T>
+        struct coord<2,T> {
+            coord () { ; }
 
-            template <typename ...T>
-            coord_data (T ..._t): data{_t...}
+            template <typename ...U>
+            coord (U ..._u): data{_u...}
             { ; }
 
             union {
-                double data[2];
+                T data[2];
                 struct {
-                    double x;
-                    double y;
+                    T x;
+                    T y;
                 };
                 struct {
-                    double r;
-                    double t;
+                    T r;
+                    T t;
                 };
             };
 
-            double& operator[] (size_t i)       { return data[i]; }
-            double  operator[] (size_t i) const { return data[i]; }
+            T& operator[] (size_t i)       { return data[i]; }
+            T  operator[] (size_t i) const { return data[i]; }
         };
 
-        template <>
-        struct coord_data<3> {
+        template <typename T>
+        struct coord<3,T> {
             union {
-                double data[3];
+                T data[3];
                 struct {
-                    double x;
-                    double y;
-                    double z;
+                    T x;
+                    T y;
+                    T z;
                 };
             };
 
-            coord_data () { ; }
+            coord () { ; }
 
-            template <typename... T>
-            coord_data (T... t): data{t...}
+            template <typename... U>
+            coord (U... u): data{u...}
             { ; }
 
-            double& operator[] (size_t i)       { return data[i]; }
-            double  operator[] (size_t i) const { return data[i]; }
+            T& operator[] (size_t i)       { return data[i]; }
+            T  operator[] (size_t i) const { return data[i]; }
         };
 
-        template <size_t S>
-        double dot (const coord_data<S> &a, const coord_data<S> &b)
+        template <size_t S, typename T>
+        T dot (const coord<S,T> &a, const coord<S,T> &b)
         {
-            double sum = 0;
+            T sum { 0 };
             for (size_t i = 0; i < S; ++i)
                 sum += a.data[i] * b.data[i];
             return sum;
         }
 
-        template <>
-        struct coord_data<4> {
+        template <typename T>
+        struct coord<4,T> {
             union {
-                double data[4];
+                T data[4];
                 struct {
-                    double x;
-                    double y;
-                    double z;
-                    double w;
+                    T x;
+                    T y;
+                    T z;
+                    T w;
                 };
             };
 
-            coord_data () { ; }
+            coord () { ; }
 
-            template <typename... T>
-            coord_data (T... t): data{t...}
+            template <typename... U>
+            coord (U... u): data{u...}
             { ; }
 
-            double& operator[] (size_t i)       { return data[i]; }
-            double  operator[] (size_t i) const { return data[i]; }
+            T& operator[] (size_t i)       { return data[i]; }
+            T  operator[] (size_t i) const { return data[i]; }
         };
 #pragma GCC diagnostic pop
     }
diff --git a/maths.hpp b/maths.hpp
index a3c8447d..294921f4 100644
--- a/maths.hpp
+++ b/maths.hpp
@@ -23,6 +23,12 @@
 #include <type_traits>
 #include <utility>
 
+template <typename T>
+T
+abs (T value)
+{ return value > 0 ? value : -value; }
+
+
 template <typename T>
 constexpr T
 pow2 (T value)
diff --git a/matrix.cpp b/matrix.cpp
index 63f05fa6..4f6403f9 100644
--- a/matrix.cpp
+++ b/matrix.cpp
@@ -326,9 +326,9 @@ matrix<T>::operator* (const matrix<T> &rhs) const {
 
 //-----------------------------------------------------------------------------
 template <typename T>
-vector<4>
-matrix<T>::operator* (const vector<4> &rhs) const {
-    return vector<4> {
+vector<4,T>
+matrix<T>::operator* (const vector<4,T> &rhs) const {
+    return vector<4,T> {
         values[0][0] * rhs.x + values[0][1] * rhs.y + values[0][2] * rhs.z + values[0][3] * rhs.w,
         values[1][0] * rhs.x + values[1][1] * rhs.y + values[1][2] * rhs.z + values[1][3] * rhs.w,
         values[2][0] * rhs.x + values[2][1] * rhs.y + values[2][2] * rhs.z + values[2][3] * rhs.w,
@@ -388,8 +388,8 @@ matrix<T>::operator== (const matrix<T> &rhs) const {
 
 //-----------------------------------------------------------------------------
 template <typename T>
-util::point<3>
-matrix<T>::to_local (const util::point<3> &p) const {
+util::point<3,T>
+matrix<T>::to_local (const util::point<3,T> &p) const {
     CHECK_SOFT (is_affine ());
 
     return { p.x * values[0][0] +
@@ -406,8 +406,8 @@ matrix<T>::to_local (const util::point<3> &p) const {
 
 //-----------------------------------------------------------------------------
 template <typename T>
-util::point<3>
-matrix<T>::to_global (const util::point<3> &p) const {
+util::point<3,T>
+matrix<T>::to_global (const util::point<3,T> &p) const {
     return inverse ().to_local (p);
 }
 
@@ -484,9 +484,9 @@ matrix<T>::perspective (T fov, T aspect, T near, T far)
 // Emulates gluLookAt
 template <typename T>
 matrix<T>
-matrix<T>::look_at (util::point<3> eye,
-                    util::point<3> centre,
-                    util::vector<3> up)
+matrix<T>::look_at (util::point<3,T> eye,
+                    util::point<3,T> centre,
+                    util::vector<3,T> up)
 {
     const auto f = eye.to (centre).normalise ();
     const auto s = cross (f, up).normalise ();
@@ -504,7 +504,7 @@ matrix<T>::look_at (util::point<3> eye,
 //-----------------------------------------------------------------------------
 template <typename T>
 matrix<T>
-matrix<T>::translate (util::vector<3> v)
+matrix<T>::translate (util::vector<3,T> v)
 {
     return { {
         { 1.f, 0.f, 0.f, v.x },
diff --git a/matrix.hpp b/matrix.hpp
index cd29ab68..7c5609fd 100644
--- a/matrix.hpp
+++ b/matrix.hpp
@@ -45,8 +45,8 @@ namespace util {
 
         T det (void) const;
 
-        matrix<T> operator* (const matrix<T>&) const;
-        vector<4> operator* (const vector<4>&) const;
+        matrix<T>   operator* (const matrix<T>&) const;
+        vector<4,T> operator* (const vector<4,T>&) const;
 
         matrix<T>& operator*= (T);
         matrix<T>  operator/  (T) const;
@@ -54,8 +54,8 @@ namespace util {
 
         bool operator== (const matrix<T>&) const;
 
-        point<3> to_local  (const point<3> &p) const;
-        point<3> to_global (const point<3> &p) const;
+        point<3,T> to_local  (const point<3,T> &p) const;
+        point<3,T> to_global (const point<3,T> &p) const;
 
         bool is_affine (void) const;
 
@@ -63,12 +63,12 @@ namespace util {
         static matrix<T> ortho   (T left, T right, T bottom, T top, T near, T far);
         static matrix<T> ortho2D (T left, T right, T bottom, T top);
         static matrix<T> perspective (T fov, T aspect, T near, T far);
-        static matrix<T> look_at (util::point<3> eye, util::point<3> centre, util::vector<3> up);
+        static matrix<T> look_at (point<3,T> eye, point<3,T> centre, vector<3,T> up);
 
         // Affine matrices
-        static matrix<T> translate (util::vector<3>);
-        static matrix<T> scale     (util::vector<3>);
-        static matrix<T> rotate    (util::vector<3> about, T angle);
+        static matrix<T> translate (util::vector<3,T>);
+        static matrix<T> scale     (util::vector<3,T>);
+        static matrix<T> rotate    (util::vector<3,T> about, T angle);
 
         // Constant matrices
         static const matrix<T> IDENTITY;
@@ -78,7 +78,7 @@ namespace util {
     typedef matrix<float> matrixf;
 
     template <typename T>
-    std::ostream& operator<< (std::ostream&, const util::matrix<T>&);
+    std::ostream& operator<< (std::ostream&, const matrix<T>&);
 }
 
 
diff --git a/noise/basis.cpp b/noise/basis.cpp
index fd3fcc04..05d5b007 100644
--- a/noise/basis.cpp
+++ b/noise/basis.cpp
@@ -46,12 +46,12 @@ generate (intmax_t x, intmax_t y, basis::seed_t seed) {
 
 
 template <>
-util::vector2
+util::vector2d
 generate (intmax_t x, intmax_t y, basis::seed_t seed) {
     auto u = permute (x, y, seed);
     auto v = permute (u ^ seed);
 
-    return util::vector2 (LUT[u], LUT[v]);
+    return util::vector2d (LUT[u], LUT[v]);
 }
 
 
@@ -165,10 +165,10 @@ gradient<L>::eval (double x, double y) const {
     // Generate the four corner values. It's not strictly necessary to
     // normalise the values, but we get a more consistent and visually
     // appealing range of outputs with normalised values.
-    vector2 p0 = generate<vector2> (x_int,     y_int,       this->seed).normalise ();
-    vector2 p1 = generate<vector2> (x_int + 1, y_int,       this->seed).normalise ();
-    vector2 p2 = generate<vector2> (x_int,     y_int + 1,   this->seed).normalise ();
-    vector2 p3 = generate<vector2> (x_int + 1, y_int + 1,   this->seed).normalise ();
+    vector2d p0 = generate<vector2d> (x_int,     y_int,       this->seed).normalise ();
+    vector2d p1 = generate<vector2d> (x_int + 1, y_int,       this->seed).normalise ();
+    vector2d p2 = generate<vector2d> (x_int,     y_int + 1,   this->seed).normalise ();
+    vector2d p3 = generate<vector2d> (x_int + 1, y_int + 1,   this->seed).normalise ();
 
     double  v0 = p0.x *        x_fac  + p0.y *  y_fac;
     double  v1 = p1.x * (x_fac - 1.0) + p1.y *  y_fac;
@@ -207,7 +207,7 @@ cellular::bounds (void) const
 
 double
 cellular::eval (double x, double y) const {
-    using util::point2;
+    using util::point2d;
 
     intmax_t x_int = static_cast<intmax_t> (x);
     intmax_t y_int = static_cast<intmax_t> (y);
@@ -228,14 +228,14 @@ cellular::eval (double x, double y) const {
     // | 6 | 7 | 8 |
     // +---+---+---+
 
-    point2 centre = { x_fac, y_fac };
+    point2d centre = { x_fac, y_fac };
     double distances[9] = { std::numeric_limits<double>::quiet_NaN () };
     double *cursor = distances;
 
     for (signed y_off = -1; y_off <= 1 ; ++y_off)
         for (signed x_off = -1; x_off <= 1; ++x_off) {
-            auto pos = point2 (double (x_off), double (y_off));
-            auto off = generate<vector2> (x_int + x_off, y_int + y_off, this->seed);
+            auto pos = point2d (double (x_off), double (y_off));
+            auto off = generate<vector2d> (x_int + x_off, y_int + y_off, this->seed);
             off += 1;
             off /= 2;
 
diff --git a/point.cpp b/point.cpp
index 873b731c..7cbc0a2e 100644
--- a/point.cpp
+++ b/point.cpp
@@ -20,8 +20,10 @@
 #include "point.hpp"
 
 #include "debug.hpp"
+#include "maths.hpp"
 
 #include <cmath>
+#include <cstdlib>
 
 using namespace std;
 
@@ -30,56 +32,61 @@ using namespace std;
 #endif
 
 //-----------------------------------------------------------------------------
-template <size_t S>
-util::point<S>::point ()
+template <size_t S, typename T>
+util::point<S,T>::point ()
 { ; }
 
 
 //-----------------------------------------------------------------------------
-template <size_t S>
-double
-util::point<S>::distance (const util::point<S> &other) const {
-    return sqrt (distance2 (other));
+template <size_t S, typename T>
+T
+util::point<S,T>::distance (const point<S,T> &other) const {
+    // TODO: this should not truncate on integral types
+    return static_cast<T> (
+        std::sqrt (distance2 (other))
+    );
 }
 
 
-template <size_t S>
-double
-util::point<S>::distance2 (const util::point<S> &other) const {
-    double total = 0.0;
+template <size_t S, typename T>
+T
+util::point<S,T>::distance2 (const point<S,T> &other) const {
+    T total { 0 };
 
     for (size_t i = 0; i < S; ++i)
         total += pow2 (this->data[i] - other.data[i]);
+
     return total;
 }
 
 
-template <size_t S>
-double
-util::point<S>::manhattan (const util::point<S> &other) const {
-    double total = 0.0;
+template <size_t S, typename T>
+T
+util::point<S,T>::manhattan (const point<S,T> &other) const {
+    T total { 0 };
 
     for (size_t i = 0; i < S; ++i)
-        total += fabs (this->data[i] - other.data[i]);
+        total += ::abs (this->data[i] - other.data[i]);
+
     return total;
 }
 
 
 //-----------------------------------------------------------------------------
-template <size_t S>
-util::point<S>&
-util::point<S>::operator*= (double f) {
-    for (double &i: this->data)
+template <size_t S, typename T>
+util::point<S,T>&
+util::point<S,T>::operator*= (T f) {
+    for (auto &i: this->data)
         i *= f;
 
     return *this;
 }
 
 
-template <size_t S>
-util::point<S>
-util::point<S>::operator* (double f) const {
-    util::point<S> out;
+template <size_t S, typename T>
+util::point<S,T>
+util::point<S,T>::operator* (T f) const {
+    util::point<S,T> out;
 
     for (size_t i = 0; i < S; ++i)
         out.data[i] = this->data[i] * f;
@@ -88,10 +95,10 @@ util::point<S>::operator* (double f) const {
 
 
 //-----------------------------------------------------------------------------
-template <size_t S>
-util::point<S>
-util::point<S>::operator- (const util::vector<S> &rhs) const {
-    util::point<S> out;
+template <size_t S, typename T>
+util::point<S,T>
+util::point<S,T>::operator- (const vector<S,T> &rhs) const {
+    util::point<S,T> out;
 
     for (size_t i = 0; i < S; ++i)
         out.data[i] = this->data[i] - rhs.data[i];
@@ -99,19 +106,19 @@ util::point<S>::operator- (const util::vector<S> &rhs) const {
 }
 
 
-template <size_t S>
-util::point<S>&
-util::point<S>::operator-= (const util::vector<S> &rhs) {
+template <size_t S, typename T>
+util::point<S,T>&
+util::point<S,T>::operator-= (const util::vector<S,T> &rhs) {
     for (size_t i = 0; i < S; ++i)
         this->data[i] -= rhs.data[i];
     return *this;
 }
 
 //-----------------------------------------------------------------------------
-template <size_t S>
-util::point<S>
-util::point<S>::operator+ (const util::vector<S> &rhs) const {
-    util::point<S> out;
+template <size_t S, typename T>
+util::point<S,T>
+util::point<S,T>::operator+ (const vector<S,T> &rhs) const {
+    util::point<S,T> out;
 
     for (size_t i = 0; i < S; ++i)
         out.data[i] = this->data[i] + rhs.data[i];
@@ -119,9 +126,9 @@ util::point<S>::operator+ (const util::vector<S> &rhs) const {
 }
 
 
-template <size_t S>
-util::point<S>&
-util::point<S>::operator+= (const util::vector<S> &rhs) {
+template <size_t S, typename T>
+util::point<S,T>&
+util::point<S,T>::operator+= (const util::vector<S,T> &rhs) {
     for (size_t i = 0; i < S; ++i)
         this->data[i] += rhs.data[i];
     return *this;
@@ -129,10 +136,10 @@ util::point<S>::operator+= (const util::vector<S> &rhs) {
 
 
 //-----------------------------------------------------------------------------
-template <size_t S>
-util::point<S>
-util::point<S>::operator- (const util::point<S> &rhs) const {
-    util::point<S> out;
+template <size_t S, typename T>
+util::point<S,T>
+util::point<S,T>::operator- (const point<S,T> &rhs) const {
+    util::point<S,T> out;
 
     for (size_t i = 0; i < S; ++i)
         out.data[i] = this->data[i] - rhs.data[i];
@@ -141,10 +148,10 @@ util::point<S>::operator- (const util::point<S> &rhs) const {
 
 
 //-----------------------------------------------------------------------------
-template <size_t S>
-util::vector<S>
-util::point<S>::to (const util::point<S> &rhs) const {
-    util::vector<S> out;
+template <size_t S, typename T>
+util::vector<S,T>
+util::point<S,T>::to (const point<S,T> &rhs) const {
+    util::vector<S,T> out;
 
     for (size_t i = 0; i < S; ++i)
         out.data[i] = rhs.data[i] - this->data[i];
@@ -153,9 +160,9 @@ util::point<S>::to (const util::point<S> &rhs) const {
 
 
 //-----------------------------------------------------------------------------
-template <size_t S>
+template <size_t S, typename T>
 void
-util::point<S>::sanity (void) const {
+util::point<S,T>::sanity (void) const {
     CHECK_SOFT (std::all_of (begin (this->data),
                              end   (this->data),
                              [] (double i) { return !std::isnan (i); }));
@@ -163,10 +170,10 @@ util::point<S>::sanity (void) const {
 
 
 //-----------------------------------------------------------------------------
-template <size_t S>
-util::point<S>
-util::operator* (const vector<S> &v, const point<S> &p) {
-    point<S> out;
+template <size_t S, typename T>
+util::point<S,T>
+util::operator* (const vector<S,T> &v, const point<S,T> &p) {
+    point<S,T> out;
     for (size_t i = 0; i < S; ++i)
         out.data[i] = p.data[i] + v.data[i];
 
@@ -174,26 +181,34 @@ util::operator* (const vector<S> &v, const point<S> &p) {
 }
 
 
-template util::point<1> util::operator* (const vector<1>&, const point<1>&);
-template util::point<2> util::operator* (const vector<2>&, const point<2>&);
-template util::point<3> util::operator* (const vector<3>&, const point<3>&);
+template util::point<1,float> util::operator* (const vector<1,float>&, const point<1,float>&);
+template util::point<2,float> util::operator* (const vector<2,float>&, const point<2,float>&);
+template util::point<3,float> util::operator* (const vector<3,float>&, const point<3,float>&);
+
+template util::point<1,double> util::operator* (const vector<1,double>&, const point<1,double>&);
+template util::point<2,double> util::operator* (const vector<2,double>&, const point<2,double>&);
+template util::point<3,double> util::operator* (const vector<3,double>&, const point<3,double>&);
 
 
-template <size_t S>
-util::point<S>
-util::operator* (const point<S> &p, const vector<S> &v)
+template <size_t S, typename T>
+util::point<S,T>
+util::operator* (const point<S,T> &p, const vector<S,T> &v)
     { return v * p; }
 
 
-template util::point<1> util::operator* (const point<1>&, const vector<1>&);
-template util::point<2> util::operator* (const point<2>&, const vector<2>&);
-template util::point<3> util::operator* (const point<3>&, const vector<3>&);
+template util::point<1,float> util::operator* (const point<1,float>&, const vector<1,float>&);
+template util::point<2,float> util::operator* (const point<2,float>&, const vector<2,float>&);
+template util::point<3,float> util::operator* (const point<3,float>&, const vector<3,float>&);
+
+template util::point<1,double> util::operator* (const point<1,double>&, const vector<1,double>&);
+template util::point<2,double> util::operator* (const point<2,double>&, const vector<2,double>&);
+template util::point<3,double> util::operator* (const point<3,double>&, const vector<3,double>&);
 
 
 //-----------------------------------------------------------------------------
-template <size_t S>
+template <size_t S, typename T>
 std::ostream&
-util::operator<< (std::ostream &os, const util::point<S> &p) {
+util::operator<< (std::ostream &os, const util::point<S,T> &p) {
     os << "point" << S << "(";
     os << p.data[0];
 
@@ -205,11 +220,23 @@ util::operator<< (std::ostream &os, const util::point<S> &p) {
 }
 
 
-template std::ostream& util::operator<< (std::ostream &os, const util::point<1>&);
-template std::ostream& util::operator<< (std::ostream &os, const util::point<2>&);
-template std::ostream& util::operator<< (std::ostream &os, const util::point<3>&);
+template std::ostream& util::operator<< (std::ostream &os, const util::point<1,float>&);
+template std::ostream& util::operator<< (std::ostream &os, const util::point<2,float>&);
+template std::ostream& util::operator<< (std::ostream &os, const util::point<3,float>&);
+
+template std::ostream& util::operator<< (std::ostream &os, const util::point<1,double>&);
+template std::ostream& util::operator<< (std::ostream &os, const util::point<2,double>&);
+template std::ostream& util::operator<< (std::ostream &os, const util::point<3,double>&);
 
 
-template struct util::point<1>;
-template struct util::point<2>;
-template struct util::point<3>;
+#define INSTANTIATE(T)              \
+template struct util::point<1,T>;   \
+template struct util::point<2,T>;   \
+template struct util::point<3,T>;
+
+INSTANTIATE(int32_t)
+INSTANTIATE(uint32_t)
+INSTANTIATE(int64_t)
+INSTANTIATE(uint64_t)
+INSTANTIATE(float)
+INSTANTIATE(double)
diff --git a/point.hpp b/point.hpp
index e3e4d5e0..a9ba0e8f 100644
--- a/point.hpp
+++ b/point.hpp
@@ -29,43 +29,51 @@
 
 namespace util {
     /// An n-dimensional position in space.
-    template <size_t S>
-    struct point : public detail::coord_data<S> {
+    template <size_t S, typename T>
+    struct point : public detail::coord<S,T> {
         static_assert (S > 0, "point dimensions must be strictly positive.");
 
         point ();
 
-        template <typename... T>
-        point (T ...t): detail::coord_data<S> {std::forward<T> (t)...} { ; }
+        template <typename... U>
+        point (U ...u): detail::coord<S,T> {std::forward<U> (u)...} { ; }
 
-        double distance  (const point &) const;
-        double distance2 (const point &) const;
-        double manhattan (const point &) const;
+        // point operators
+        T distance  (const point &) const;
+        T distance2 (const point &) const;
+        T manhattan (const point &) const;
 
-        point<S>& operator*= (double);
-        point<S>  operator*  (double) const;
-        point<S>  operator-  (const point<S>&) const;
+        vector<S,T> to (const point&) const;
 
-        point<S>  operator-  (const util::vector<S>&) const;
-        point<S>& operator-= (const util::vector<S>&);
-        point<S>  operator+  (const util::vector<S>&) const;
-        point<S>& operator+= (const util::vector<S>&);
+        // arithetic operators
+        point<S,T>& operator*= (T);
+        point<S,T>  operator*  (T) const;
+        point<S,T>  operator-  (const point<S,T>&) const;
 
-        util::vector<S> to (const point<S>&) const;
+        point<S,T>  operator-  (const vector<S,T>&) const;
+        point<S,T>& operator-= (const vector<S,T>&);
+        point<S,T>  operator+  (const vector<S,T>&) const;
+        point<S,T>& operator+= (const vector<S,T>&);
 
-        template <size_t D> point<D> redim (void) const;
+        template <size_t D> point<D,T> redim (void) const;
 
         void sanity (void) const;
     };
 
-    typedef point<2> point2;
-    typedef point<3> point3;
+    // free maths operators
+    template <size_t S, typename T> point<S,T> operator* (const vector<S,T>&, const point<S,T>&);
+    template <size_t S, typename T> point<S,T> operator* (const point<S,T>&, const vector<S,T>&);
 
-    template <size_t S> point<S> operator* (const vector<S>&, const point<S>&);
-    template <size_t S> point<S> operator* (const point<S>&, const vector<S>&);
+    // iostream operators
+    template <size_t S, typename T>
+    std::ostream& operator<< (std::ostream&, const point<S,T>&);
 
-    template <size_t S>
-    std::ostream& operator<< (std::ostream&, const util::point<S>&);
+    // Convenience typedefs
+    typedef point<2,float> point2f;
+    typedef point<3,float> point3f;
+
+    typedef point<2,double> point2d;
+    typedef point<3,double> point3d;
 }
 
 #include "point.ipp"
diff --git a/point.ipp b/point.ipp
index 9ff1a486..ed332335 100644
--- a/point.ipp
+++ b/point.ipp
@@ -20,10 +20,10 @@
 #include <algorithm>
 
 namespace util {
-    template<size_t S>
+    template<size_t S, typename T>
     template<size_t D>
-    point<D> point<S>::redim (void) const {
-        point<D> out;
+    point<D,T> point<S,T>::redim (void) const {
+        point<D,T> out;
         std::copy (std::begin (this->data), std::begin (this->data) + D, std::begin (out.data));
         return out;
     }
diff --git a/quaternion.cpp b/quaternion.cpp
index b675479a..cc9aa335 100644
--- a/quaternion.cpp
+++ b/quaternion.cpp
@@ -25,30 +25,34 @@ using util::quaternion;
 
 
 //-----------------------------------------------------------------------------
-const quaternion quaternion::IDENTITY = { 1.0, 0.0, 0.0, 0.0 };
+template<> const quaternion<float> quaternion<float>::IDENTITY = { 1, 0, 0, 0 };
+template<> const quaternion<double> quaternion<double>::IDENTITY = { 1, 0, 0, 0 };
 
 
 //-----------------------------------------------------------------------------
-quaternion
-quaternion::rotation (double radians, vector<3> axis) {
-    radians /= 2.0;
+template <typename T>
+quaternion<T>
+quaternion<T>::rotation (T radians, vector<3,T> axis) {
+    radians /= T{2};
     axis.normalise ();
 
     return {
-        cos (radians),
-        axis.x * sin (radians),
-        axis.y * sin (radians),
-        axis.z * sin (radians)
+        std::cos (radians),
+        axis.x * std::sin (radians),
+        axis.y * std::sin (radians),
+        axis.z * std::sin (radians)
     };
 }
 
 
-quaternion
-quaternion::rotation (vector<3> from, vector<3> to) {
+//-----------------------------------------------------------------------------
+template <typename T>
+quaternion<T>
+quaternion<T>::rotation (vector<3,T> from, vector<3,T> to) {
     auto v = util::cross (from, to);
 
     return {
-        acos (from.dot (to)),
+        std::acos (from.dot (to)),
         v.x,
         v.y,
         v.z
@@ -56,8 +60,10 @@ quaternion::rotation (vector<3> from, vector<3> to) {
 }
 
 
-quaternion
-quaternion::operator* (const quaternion &rhs) const {
+//-----------------------------------------------------------------------------
+template <typename T>
+quaternion<T>
+quaternion<T>::operator* (const quaternion<T> &rhs) const {
     return {
         w * rhs.w - (x * rhs.x + y * rhs.y + z * rhs.z),
         w * rhs.x + rhs.w * x + y * rhs.z - z * rhs.y,
@@ -66,3 +72,7 @@ quaternion::operator* (const quaternion &rhs) const {
     };
 }
 
+
+//-----------------------------------------------------------------------------
+template struct quaternion<float>;
+template struct quaternion<double>;
diff --git a/quaternion.hpp b/quaternion.hpp
index 66b99413..99b7eb36 100644
--- a/quaternion.hpp
+++ b/quaternion.hpp
@@ -23,13 +23,14 @@
 #include "vector.hpp"
 
 namespace util {
+    template <typename T>
     struct quaternion {
-        double w, x, y, z;
+        T w, x, y, z;
 
         static const quaternion IDENTITY;
 
-        static quaternion rotation (double radians, vector<3> axis);
-        static quaternion rotation (vector<3> from, vector<3> to);
+        static quaternion rotation (T radians, vector<3,T> axis);
+        static quaternion rotation (vector<3,T> from, vector<3,T> to);
 
         quaternion operator* (const quaternion&) const;
     };
diff --git a/region.cpp b/region.cpp
index f61212e0..1f545622 100644
--- a/region.cpp
+++ b/region.cpp
@@ -53,18 +53,18 @@ region<T>::area (void) const
 template <typename T>
 typename region<T>::size_type
 region<T>::diameter (void) const {
-    return static_cast<size_type> (sqrt (w * w + h * h));
+    return static_cast<size_type> (std::sqrt (w * w + h * h));
 }
 
 
 template <typename T>
 void
-region<T>::scale (double factor) {
-    x -= static_cast<T> ((w * factor - w) / 2.0);
-    y -= static_cast<T> ((h * factor - h) / 2.0);
+region<T>::scale (T factor) {
+    x -= (w * factor - w) / T{2};
+    y -= (h * factor - h) / T{2};
 
-    w = static_cast<T> (w * factor);
-    h = static_cast<T> (h * factor);
+    w = w * factor;
+    h = h * factor;
 }
 
 
@@ -76,25 +76,25 @@ region<T>::empty (void) const
 
 //-----------------------------------------------------------------------------
 template <typename T>
-point<2>
+point<2,T>
 region<T>::base (void) const {
-    return { static_cast<double> (x), static_cast<double> (y) };
+    return { x, y };
 }
 
 
 template <typename T>
-point<2>
+point<2,T>
 region<T>::centre (void) const {
-    double cx = x + static_cast<T>(w / 2.0),
-           cy = y + static_cast<T>(h / 2.0);
+    T cx = x + w / T{2},
+      cy = y + h / T{2};
 
-    return { cx, cy };
+    return point<2,T> { cx, cy };
 }
 
 
 template <typename T>
-point<2>
-region<T>::closest (point<2> p) const {
+point<2,T>
+region<T>::closest (point<2,T> p) const {
     return {
         p.x < x     ? x     :
         p.x > x + w ? x + w :
@@ -110,7 +110,7 @@ region<T>::closest (point<2> p) const {
 //-----------------------------------------------------------------------------
 template <typename T>
 bool
-region<T>::includes (const point<2> &p) const {
+region<T>::includes (const point<2,T> &p) const {
     return p.x >= x &&
            p.y >= y &&
            p.x - x <= w &&
@@ -120,7 +120,7 @@ region<T>::includes (const point<2> &p) const {
 
 template <typename T>
 bool
-region<T>::contains (const point<2> &p) const {
+region<T>::contains (const point<2,T> &p) const {
     return p.x > x &&
            p.y > y &&
            p.x - x < w &&
@@ -143,18 +143,19 @@ region<T>::overlaps (const region<T> &rhs) const {
 //-----------------------------------------------------------------------------
 template <typename T>
 void
-region<T>::constrain (point2 &p) const {
-    p.x = std::min (std::max (static_cast<T> (p.x), x), x + w);
-    p.y = std::min (std::max (static_cast<T> (p.y), y), y + h);
+region<T>::constrain (point<2,T> &p) const {
+    p.x = std::min (std::max (p.x, x), x + w);
+    p.y = std::min (std::max (p.y, y), y + h);
 }
 
 
 template <typename T>
-point2
-region<T>::constrained (const point2 &p) const {
-    point2 v;
-    v.x = std::min (std::max (static_cast<T> (p.x), x), x + w);
-    v.y = std::min (std::max (static_cast<T> (p.y), y), y + h);
+point<2,T>
+region<T>::constrained (const point<2,T> &p) const
+{
+    point<2,T> v;
+    v.x = std::min (std::max (p.x, x), x + w);
+    v.y = std::min (std::max (p.y, y), y + h);
 
     return v;
 }
@@ -305,5 +306,6 @@ namespace util {
     template std::ostream& operator<< (std::ostream&, const region< int64_t>&);
     template std::ostream& operator<< (std::ostream&, const region<uint32_t>&);
     template std::ostream& operator<< (std::ostream&, const region<uint64_t>&);
+    template std::ostream& operator<< (std::ostream&, const region<   float>&);
     template std::ostream& operator<< (std::ostream&, const region<  double>&);
 }
diff --git a/region.hpp b/region.hpp
index 4937b68e..947ae7b2 100644
--- a/region.hpp
+++ b/region.hpp
@@ -44,20 +44,20 @@ namespace util {
 
         size_type area     (void) const;
         size_type diameter (void) const;
-        void scale (double factor);
+        void scale (T factor);
 
         bool empty (void) const;
 
-        point2 base    (void) const;
-        point2 centre  (void) const;
-        point2 closest (point2) const;
+        point<2,T> base    (void) const;
+        point<2,T> centre  (void) const;
+        point<2,T> closest (point<2,T>) const;
 
-        bool includes (const point2&) const; // inclusive of borders
-        bool contains (const point2&) const; // exclusive of borders
+        bool includes (const point<2,T>&) const; // inclusive of borders
+        bool contains (const point<2,T>&) const; // exclusive of borders
         bool overlaps (const region<T>&) const; // exclusive of borders
 
-        void constrain (point2&) const;
-        point2 constrained (const point2&) const;
+        void constrain (point<2,T>&) const;
+        point<2,T> constrained (const point<2,T>&) const;
 
         region overlap (const region<T>&) const;
 
diff --git a/test/matrix.cpp b/test/matrix.cpp
index 052ef48e..c7abf3d6 100644
--- a/test/matrix.cpp
+++ b/test/matrix.cpp
@@ -8,7 +8,7 @@ int
 main (int, char **) {
     {
         // Identity matrix-vector multiplication
-        auto v = util::vector<4> { 1, 2, 3, 4 };
+        auto v = util::vector<4,float> { 1.f, 2.f, 3.f, 4.f };
         auto r = util::matrix<float>::IDENTITY * v;
         CHECK_EQ (r, v);
     }
@@ -22,7 +22,7 @@ main (int, char **) {
             { 13, 14, 15, 16 }
         } };
 
-        util::vector<4> v { 1, 2, 3, 4 };
+        util::vector<4,float> v { 1.f, 2.f, 3.f, 4.f };
 
         auto r = m * v;
 
diff --git a/test/point.cpp b/test/point.cpp
index e995d8cf..682b88ab 100644
--- a/test/point.cpp
+++ b/test/point.cpp
@@ -8,8 +8,8 @@ using namespace util;
 
 int
 main (int, char**) {
-    const point<3> p(0.0, 1.0, 2.0);
-    const point<2> q = p.redim<2> ();
+    const point3f p(0.f, 1.f, 2.f);
+    const point2f q = p.redim<2> ();
 
     CHECK_EQ (q.data[0], p.data[0]);
     CHECK_EQ (q.data[1], p.data[1]);
diff --git a/test/region.cpp b/test/region.cpp
index afd7427c..50e25a5f 100644
--- a/test/region.cpp
+++ b/test/region.cpp
@@ -3,7 +3,8 @@
 #include "../debug.hpp"
 
 using util::region;
-using util::point2;
+using util::point;
+using util::point2d;
 
 int
 main (int, char **) {
@@ -20,16 +21,16 @@ main (int, char **) {
     CHECK_EQ (region<double>::UNIT.area (), 1.0);
     CHECK_EQ (region< float>::UNIT.area (), 1.0f);
 
-    CHECK_HARD (region<int> (0, 0, 2, 2).includes (point2(1.0, 1.0)));
-    CHECK_HARD (region<int> (0, 0, 2, 2).includes (point2(0.0, 0.0)));
-    CHECK_HARD (region<int> (0, 0, 2, 2).includes (point2(2.0, 2.0)));
+    CHECK_HARD (region<int> (0, 0, 2, 2).includes (point<2,int>(1, 1)));
+    CHECK_HARD (region<int> (0, 0, 2, 2).includes (point<2,int>(0, 0)));
+    CHECK_HARD (region<int> (0, 0, 2, 2).includes (point<2,int>(2, 2)));
 
-    CHECK_HARD ( region<int> (0, 0, 2, 2).contains (point2(1.0, 1.0)));
-    CHECK_HARD (!region<int> (0, 0, 2, 2).contains (point2(0.0, 0.0)));
-    CHECK_HARD (!region<int> (0, 0, 2, 2).contains (point2(2.0, 2.0)));
+    CHECK_HARD ( region<int> (0, 0, 2, 2).contains (point<2,int>(1, 1)));
+    CHECK_HARD (!region<int> (0, 0, 2, 2).contains (point<2,int>(0, 0)));
+    CHECK_HARD (!region<int> (0, 0, 2, 2).contains (point<2,int>(2, 2)));
 
-    CHECK_HARD (region<intmax_t> (0, 0, 10, 10).includes (point2 (0.4, 0.01)));
-    CHECK_HARD (region<intmax_t> (0, 0, 10, 10).contains (point2 (0.4, 0.01)));
+    //CHECK_HARD (region<intmax_t> (0, 0, 10, 10).includes (point2d (0.4, 0.01)));
+    //CHECK_HARD (region<intmax_t> (0, 0, 10, 10).contains (point2d (0.4, 0.01)));
 
     return 0;
 }
diff --git a/vector.cpp b/vector.cpp
index 8dc9428d..63201679 100644
--- a/vector.cpp
+++ b/vector.cpp
@@ -39,46 +39,49 @@ using std::end;
 
 
 //-----------------------------------------------------------------------------
-template <size_t S>
-util::vector<S>::vector ()
+template <size_t S, typename T>
+util::vector<S,T>::vector ()
 { ; }
 
 
 //-----------------------------------------------------------------------------
-template <size_t S>
-util::vector<S>
-util::vector<S>::operator* (double rhs) const {
-    util::vector<S> out;
+template <size_t S, typename T>
+util::vector<S,T>
+util::vector<S,T>::operator* (T rhs) const {
+    util::vector<S,T> out;
+
     for (size_t i = 0; i < S; ++i)
         out.data[i] = this->data[i] * rhs;
+
     return out;
 }
 
 
-template <size_t S>
-util::vector<S>&
-util::vector<S>::operator*= (double rhs) {
-    for (double &i: this->data)
+template <size_t S, typename T>
+util::vector<S,T>&
+util::vector<S,T>::operator*= (T rhs) {
+    for (auto &i: this->data)
         i *= rhs;
 
     return *this;
 }
 
 
-template <size_t S>
-util::vector<S>
-util::vector<S>::operator* (const util::vector<S> &rhs) const {
-    util::vector<S> out;
+template <size_t S, typename T>
+util::vector<S,T>
+util::vector<S,T>::operator* (const vector<S,T> &rhs) const {
+    util::vector<S,T> out;
 
     for (size_t i = 0; i < S; ++i)
         out.data[i] = this->data[i] * rhs.data[i];
+
     return out;
 }
 
 
-template <size_t S>
-util::vector<S>&
-util::vector<S>::operator*= (const util::vector<S> &rhs) {
+template <size_t S, typename T>
+util::vector<S,T>&
+util::vector<S,T>::operator*= (const vector<S,T> &rhs) {
     for (size_t i = 0; i < S; ++i)
         this->data[i] *= rhs.data[i];
 
@@ -87,10 +90,10 @@ util::vector<S>::operator*= (const util::vector<S> &rhs) {
 
 
 //-----------------------------------------------------------------------------
-template <size_t S>
-util::vector<S>
-util::vector<S>::operator/ (double rhs) const {
-    util::vector<S> out;
+template <size_t S, typename T>
+util::vector<S,T>
+util::vector<S,T>::operator/ (T rhs) const {
+    util::vector<S,T> out;
 
     for (size_t i = 0; i < S; ++i)
         out.data[i] = this->data[i] / rhs;
@@ -98,9 +101,9 @@ util::vector<S>::operator/ (double rhs) const {
 }
 
 
-template <size_t S>
-util::vector<S>&
-util::vector<S>::operator/= (double rhs) {
+template <size_t S, typename T>
+util::vector<S,T>&
+util::vector<S,T>::operator/= (T rhs) {
     for (size_t i = 0; i < S; ++i)
         this->data[i] /= rhs;
     return *this;
@@ -108,10 +111,10 @@ util::vector<S>::operator/= (double rhs) {
 
 
 //-----------------------------------------------------------------------------
-template <size_t S>
-util::vector<S>
-util::vector<S>::operator+ (const util::vector<S> &rhs) const {
-    util::vector<S> out;
+template <size_t S, typename T>
+util::vector<S,T>
+util::vector<S,T>::operator+ (const util::vector<S,T> &rhs) const {
+    util::vector<S,T> out;
 
     for (size_t i = 0; i < S; ++i)
         out.data[i] = this->data[i] + rhs.data[i];
@@ -119,10 +122,10 @@ util::vector<S>::operator+ (const util::vector<S> &rhs) const {
 }
 
 
-template <size_t S>
-util::vector<S>
-util::vector<S>::operator+ (double rhs) const {
-    util::vector<S> out;
+template <size_t S, typename T>
+util::vector<S,T>
+util::vector<S,T>::operator+ (T rhs) const {
+    util::vector<S,T> out;
 
     for (size_t i = 0; i < S; ++i)
         out.data[i] = this->data[i] + rhs;
@@ -130,9 +133,9 @@ util::vector<S>::operator+ (double rhs) const {
 }
 
 
-template <size_t S>
-util::vector<S>&
-util::vector<S>::operator+= (const util::vector<S> &rhs) {
+template <size_t S, typename T>
+util::vector<S,T>&
+util::vector<S,T>::operator+= (const util::vector<S,T> &rhs) {
     for (size_t i = 0; i < S; ++i)
         this->data[i] += rhs.data[i];
 
@@ -140,9 +143,9 @@ util::vector<S>::operator+= (const util::vector<S> &rhs) {
 }
 
 
-template <size_t S>
-util::vector<S>&
-util::vector<S>::operator+= (double rhs) {
+template <size_t S, typename T>
+util::vector<S,T>&
+util::vector<S,T>::operator+= (T rhs) {
     for (size_t i = 0; i < S; ++i)
         this->data[i] += rhs;
     return *this;
@@ -150,10 +153,10 @@ util::vector<S>::operator+= (double rhs) {
 
 
 //-----------------------------------------------------------------------------
-template <size_t S>
-util::vector<S>
-util::vector<S>::operator- (void) const {
-    util::vector<S> out;
+template <size_t S, typename T>
+util::vector<S,T>
+util::vector<S,T>::operator- (void) const {
+    util::vector<S,T> out;
 
     for (size_t i = 0; i < S; ++i)
         out.data[i] = -this->data[i];
@@ -161,10 +164,10 @@ util::vector<S>::operator- (void) const {
 }
 
 
-template <size_t S>
-util::vector<S>
-util::vector<S>::operator- (const util::vector<S> &rhs) const {
-    util::vector<S> out;
+template <size_t S, typename T>
+util::vector<S,T>
+util::vector<S,T>::operator- (const util::vector<S,T> &rhs) const {
+    util::vector<S,T> out;
 
     for (size_t i = 0; i < S; ++i)
         out.data[i] = this->data[i] - rhs.data[i];
@@ -172,9 +175,9 @@ util::vector<S>::operator- (const util::vector<S> &rhs) const {
 }
 
 
-template <size_t S>
-util::vector<S>&
-util::vector<S>::operator-= (const util::vector<S> &rhs) {
+template <size_t S, typename T>
+util::vector<S,T>&
+util::vector<S,T>::operator-= (const util::vector<S,T> &rhs) {
     for (size_t i = 0; i < S; ++i)
         this->data[i] -= rhs.data[i];
 
@@ -182,10 +185,10 @@ util::vector<S>::operator-= (const util::vector<S> &rhs) {
 }
 
 
-template <size_t S>
-util::vector<S>
-util::vector<S>::operator- (double rhs) const {
-    util::vector<S> out;
+template <size_t S, typename T>
+util::vector<S,T>
+util::vector<S,T>::operator- (T rhs) const {
+    util::vector<S,T> out;
 
     for (size_t i = 0; i < S; ++i)
         out.data[i] = this->data[i] - rhs;
@@ -193,9 +196,9 @@ util::vector<S>::operator- (double rhs) const {
 }
 
 
-template <size_t S>
-util::vector<S>&
-util::vector<S>::operator-= (double rhs) {
+template <size_t S, typename T>
+util::vector<S,T>&
+util::vector<S,T>::operator-= (T rhs) {
     for (size_t i = 0; i < S; ++i)
         this->data[i] -= rhs;
 
@@ -204,18 +207,18 @@ util::vector<S>::operator-= (double rhs) {
 
 
 //-----------------------------------------------------------------------------
-template <size_t S>
-util::vector<S>&
-util::vector<S>::operator= (const util::vector<S> &rhs) {
+template <size_t S, typename T>
+util::vector<S,T>&
+util::vector<S,T>::operator= (const util::vector<S,T> &rhs) {
     std::copy (begin (rhs.data), end (rhs.data), begin (this->data));
 
     return *this;
 }
 
 
-template <size_t S>
+template <size_t S, typename T>
 bool
-util::vector<S>::operator== (const util::vector<S> &rhs) const {
+util::vector<S,T>::operator== (const util::vector<S,T> &rhs) const {
     for (size_t i = 0; i < S; ++i)
         if (!almost_equal (this->data[i], rhs.data[i]))
             return false;
@@ -225,27 +228,28 @@ util::vector<S>::operator== (const util::vector<S> &rhs) const {
 
 
 //-----------------------------------------------------------------------------
-template <size_t S>
-double
-util::vector<S>::magnitude (void) const {
-    return sqrt (magnitude2 ());
+template <size_t S, typename T>
+T
+util::vector<S,T>::magnitude (void) const {
+    // TODO: this should not truncate for integral types
+    return static_cast<T> (std::sqrt (magnitude2 ()));
 }
 
 
-template <size_t S>
-double
-util::vector<S>::magnitude2 (void) const {
-    double total = 0.0;
+template <size_t S, typename T>
+T
+util::vector<S,T>::magnitude2 (void) const {
+    T total { 0 };
     for (size_t i = 0; i < S; ++i)
         total += pow2 (this->data[i]);
     return total;
 }
 
 
-template <size_t S>
-util::vector<S>&
-util::vector<S>::normalise (void) {
-    double mag = magnitude ();
+template <size_t S, typename T>
+util::vector<S,T>&
+util::vector<S,T>::normalise (void) {
+    T mag = magnitude ();
 
     for (size_t i = 0; i < S; ++i)
         this->data[i] /= mag;
@@ -254,11 +258,11 @@ util::vector<S>::normalise (void) {
 }
 
 
-template <size_t S>
-util::vector<S>
-util::vector<S>::normalised (void) const {
-    double mag = magnitude ();
-    util::vector<S> out;
+template <size_t S, typename T>
+util::vector<S,T>
+util::vector<S,T>::normalised (void) const {
+    T mag = magnitude ();
+    util::vector<S,T> out;
 
     for (size_t i = 0; i < S; ++i)
         out.data[i] = this->data[i] / mag;
@@ -267,39 +271,47 @@ util::vector<S>::normalised (void) const {
 }
 
 //-----------------------------------------------------------------------------
-util::vector<2>
-util::polar_to_cartesian (const util::vector<2> &v) {
-    return util::vector<2> {
+template <typename T>
+util::vector<2,T>
+util::polar_to_cartesian (const util::vector<2,T> &v) {
+    return util::vector<2,T> {
         v.r * std::cos (v.t),
         v.r * std::sin (v.t)
     };
 }
 
 //-----------------------------------------------------------------------------
-template <size_t S>
-double
-util::vector<S>::dot (const util::vector<S> &rhs) const {
-    double total = 0.0;
+template <size_t S, typename T>
+T
+util::vector<S,T>::dot (const util::vector<S,T> &rhs) const {
+    T total { 0 };
     for (size_t i = 0; i < S; ++i)
         total += this->data[i] * rhs.data[i];
     return total;
 }
 
 
-util::vector<3>
-util::cross (const util::vector<3> &a, const util::vector<3> &b) {
-    return util::vector<3> {
+template <typename T>
+util::vector<3,T>
+util::cross (const util::vector<3,T> &a,
+             const util::vector<3,T> &b)
+{
+    return util::vector<3,T> {
         a.y * b.z - a.z * b.y,
         a.z * b.x - a.x * b.z,
         a.x * b.y - a.y * b.x
     };
 }
 
+template util::vector3f util::cross(const util::vector3f&, const util::vector3f&);
+template util::vector3d util::cross(const util::vector3d&, const util::vector3d&);
+
 
 //-----------------------------------------------------------------------------
-util::vector<3>
-util::spherical_to_cartesian (const util::vector<3> &s) {
-    return util::vector<3> {
+template <typename T>
+util::vector<3,T>
+util::spherical_to_cartesian (const util::vector<3,T> &s) {
+    return util::vector<3,T> {
         s.x * sin (s.y) * cos (s.z),
         s.x * sin (s.y) * sin (s.z),
         s.x * cos (s.y),
@@ -307,11 +319,12 @@ util::spherical_to_cartesian (const util::vector<3> &s) {
 }
 
 
-util::vector<3>
-util::cartesian_to_spherical (const util::vector<3> &c) {
-    double mag = c.magnitude ();
+template <typename T>
+util::vector<3,T>
+util::cartesian_to_spherical (const util::vector<3,T> &c) {
+    T mag = c.magnitude ();
 
-    return util::vector<3> {
+    return util::vector<3,T> {
         mag,
         acos  (c.z / mag),
         atan2 (c.y, c.x)
@@ -320,63 +333,77 @@ util::cartesian_to_spherical (const util::vector<3> &c) {
 
 
 //-----------------------------------------------------------------------------
-template <size_t S>
+template <size_t S, typename T>
 bool
-util::vector<S>::is_zero (void) const {
+util::vector<S,T>::is_zero (void) const {
     return std::all_of (begin (this->data),
                         end   (this->data),
-                        [] (double i) { return almost_zero (i); });
+                        [] (T i) { return almost_zero (i); });
 }
 
 
 //-----------------------------------------------------------------------------
-template <size_t S>
+template <size_t S, typename T>
 void
-util::vector<S>::sanity (void) const {
+util::vector<S,T>::sanity (void) const {
     CHECK_SOFT (std::all_of (begin (this->data),
                              end   (this->data),
-                             [] (double i) { return !std::isnan (i); }));
+                             [] (T i) { return !std::isnan (i); }));
 }
 
 
 //-----------------------------------------------------------------------------
-template <size_t S>
-util::vector<S>
-util::operator* (double a, const util::vector<S> &b)
+template <size_t S, typename T>
+util::vector<S,T>
+util::operator* (T a, const util::vector<S,T> &b)
     { return b * a; }
 
 
-template util::vector<1> util::operator* (double, const util::vector<1>&);
-template util::vector<2> util::operator* (double, const util::vector<2>&);
-template util::vector<3> util::operator* (double, const util::vector<3>&);
+template util::vector<1,float> util::operator* (float, const util::vector<1,float>&);
+template util::vector<2,float> util::operator* (float, const util::vector<2,float>&);
+template util::vector<3,float> util::operator* (float, const util::vector<3,float>&);
 
-
-template <size_t S>
-util::vector<S>
-util::operator+ (double a, const util::vector<S> &b)
-    { return b + a; }
-
-
-template util::vector<1> util::operator+ (double, const util::vector<1>&);
-template util::vector<2> util::operator+ (double, const util::vector<2>&);
-template util::vector<3> util::operator+ (double, const util::vector<3>&);
-
-
-template <size_t S>
-util::vector<S>
-util::operator- (double a, const util::vector<S> &b)
-    { return a + (-b); }
-
-
-template util::vector<1> util::operator- (double, const util::vector<1>&);
-template util::vector<2> util::operator- (double, const util::vector<2>&);
-template util::vector<3> util::operator- (double, const util::vector<3>&);
+template util::vector<1,double> util::operator* (double, const util::vector<1,double>&);
+template util::vector<2,double> util::operator* (double, const util::vector<2,double>&);
+template util::vector<3,double> util::operator* (double, const util::vector<3,double>&);
 
 
 //-----------------------------------------------------------------------------
-template <size_t S>
+template <size_t S, typename T>
+util::vector<S,T>
+util::operator+ (T a, const util::vector<S,T> &b)
+    { return b + a; }
+
+
+template util::vector<1,float> util::operator+ (float, const util::vector<1,float>&);
+template util::vector<2,float> util::operator+ (float, const util::vector<2,float>&);
+template util::vector<3,float> util::operator+ (float, const util::vector<3,float>&);
+
+template util::vector<1,double> util::operator+ (double, const util::vector<1,double>&);
+template util::vector<2,double> util::operator+ (double, const util::vector<2,double>&);
+template util::vector<3,double> util::operator+ (double, const util::vector<3,double>&);
+
+
+//-----------------------------------------------------------------------------
+template <size_t S, typename T>
+util::vector<S,T>
+util::operator- (T a, const util::vector<S,T> &b)
+    { return a + (-b); }
+
+
+template util::vector<1,float> util::operator- (float, const util::vector<1,float>&);
+template util::vector<2,float> util::operator- (float, const util::vector<2,float>&);
+template util::vector<3,float> util::operator- (float, const util::vector<3,float>&);
+
+template util::vector<1,double> util::operator- (double, const util::vector<1,double>&);
+template util::vector<2,double> util::operator- (double, const util::vector<2,double>&);
+template util::vector<3,double> util::operator- (double, const util::vector<3,double>&);
+
+
+//-----------------------------------------------------------------------------
+template <size_t S, typename T>
 std::ostream&
-util::operator<< (std::ostream &os, const util::vector<S> &v) {
+util::operator<< (std::ostream &os, const util::vector<S,T> &v) {
     os << "vec" << S << "(" << v.data[0];
     for (size_t i = 1; i < S; ++i)
         os << ", " << v.data[i];
@@ -385,16 +412,21 @@ util::operator<< (std::ostream &os, const util::vector<S> &v) {
 }
 
 
-template std::ostream& util::operator<< (std::ostream&, const util::vector<1> &v);
-template std::ostream& util::operator<< (std::ostream&, const util::vector<2> &v);
-template std::ostream& util::operator<< (std::ostream&, const util::vector<3> &v);
-template std::ostream& util::operator<< (std::ostream&, const util::vector<4> &v);
+template std::ostream& util::operator<< (std::ostream&, const util::vector<1,float> &v);
+template std::ostream& util::operator<< (std::ostream&, const util::vector<2,float> &v);
+template std::ostream& util::operator<< (std::ostream&, const util::vector<3,float> &v);
+template std::ostream& util::operator<< (std::ostream&, const util::vector<4,float> &v);
+
+template std::ostream& util::operator<< (std::ostream&, const util::vector<1,double> &v);
+template std::ostream& util::operator<< (std::ostream&, const util::vector<2,double> &v);
+template std::ostream& util::operator<< (std::ostream&, const util::vector<3,double> &v);
+template std::ostream& util::operator<< (std::ostream&, const util::vector<4,double> &v);
 
 
 //-----------------------------------------------------------------------------
-template <size_t S>
+template <size_t S, typename T>
 const json::node&
-util::operator>> (const json::node &node, util::vector<S> &v) {
+util::operator>> (const json::node &node, util::vector<S,T> &v) {
     const json::array &array = node.as_array ();
     if (array.size () != S)
         throw std::runtime_error ("Invalid dimensionality for json-to-vector");
@@ -403,30 +435,47 @@ util::operator>> (const json::node &node, util::vector<S> &v) {
     // compiler error at this point in release mode, so we dumb it down a
     // little.
     for (size_t i = 0; i < array.size (); ++i)
-        v.data[i] = array[i].as_number ().native ();
+        v.data[i] = static_cast<T> (array[i].as_number ().native ());
 
     return node;
 }
 
 
-template const json::node& util::operator>> (const json::node&, util::vector<1>&);
-template const json::node& util::operator>> (const json::node&, util::vector<2>&);
-template const json::node& util::operator>> (const json::node&, util::vector<3>&);
-template const json::node& util::operator>> (const json::node&, util::vector<4>&);
+template const json::node& util::operator>> (const json::node&, util::vector<1,float>&);
+template const json::node& util::operator>> (const json::node&, util::vector<2,float>&);
+template const json::node& util::operator>> (const json::node&, util::vector<3,float>&);
+template const json::node& util::operator>> (const json::node&, util::vector<4,float>&);
+
+template const json::node& util::operator>> (const json::node&, util::vector<1,double>&);
+template const json::node& util::operator>> (const json::node&, util::vector<2,double>&);
+template const json::node& util::operator>> (const json::node&, util::vector<3,double>&);
+template const json::node& util::operator>> (const json::node&, util::vector<4,double>&);
 
 
 //-----------------------------------------------------------------------------
-template struct util::vector<1>;
-template struct util::vector<2>;
-template struct util::vector<3>;
-template struct util::vector<4>;
+#define INSTANTIATE(T)              \
+template struct util::vector<1,T>;  \
+template struct util::vector<2,T>;  \
+template struct util::vector<3,T>;  \
+template struct util::vector<4,T>;
 
+INSTANTIATE(uint32_t)
+INSTANTIATE(int32_t)
+INSTANTIATE(uint64_t)
+INSTANTIATE(int64_t)
+INSTANTIATE(float)
+INSTANTIATE(double)
 
 //-----------------------------------------------------------------------------
 namespace util {
-    template <> vector<1> random (void) { util::vector<1> out; randomise (out.data); return out; }
-    template <> vector<2> random (void) { util::vector<2> out; randomise (out.data); return out; }
-    template <> vector<3> random (void) { util::vector<3> out; randomise (out.data); return out; }
-    template <> vector<4> random (void) { util::vector<4> out; randomise (out.data); return out; }
+    template <> vector<1,float> random (void) { util::vector<1,float> out; randomise (out.data); return out; }
+    template <> vector<2,float> random (void) { util::vector<2,float> out; randomise (out.data); return out; }
+    template <> vector<3,float> random (void) { util::vector<3,float> out; randomise (out.data); return out; }
+    template <> vector<4,float> random (void) { util::vector<4,float> out; randomise (out.data); return out; }
+
+    template <> vector<1,double> random (void) { util::vector<1,double> out; randomise (out.data); return out; }
+    template <> vector<2,double> random (void) { util::vector<2,double> out; randomise (out.data); return out; }
+    template <> vector<3,double> random (void) { util::vector<3,double> out; randomise (out.data); return out; }
+    template <> vector<4,double> random (void) { util::vector<4,double> out; randomise (out.data); return out; }
 }
 
diff --git a/vector.hpp b/vector.hpp
index 76d372a7..5e5d49b0 100644
--- a/vector.hpp
+++ b/vector.hpp
@@ -28,74 +28,83 @@
 #include <initializer_list>
 
 namespace util {
-    template <size_t S>
-    struct vector : public detail::coord_data<S> {
+    template <size_t S, typename T>
+    struct vector : public detail::coord<S, T> {
         static_assert (S > 0, "vector dimensions must be strictly positive");
 
         vector ();
 
-        template <typename... T>
-        explicit vector (T ...t): detail::coord_data<S> {std::forward<T> (t)...} { ; }
+        template <typename... U>
+        explicit vector (U ...u): detail::coord<S, T> {std::forward<U> (u)...} { ; }
 
-        util::vector<S>  operator* (double) const;
-        util::vector<S>& operator*=(double);
+        // arithmetic operators
+        vector<S,T>  operator* (T) const;
+        vector<S,T>& operator*=(T);
 
-        util::vector<S>  operator/ (double) const;
-        util::vector<S>& operator/=(double);
+        vector<S,T>  operator/ (T) const;
+        vector<S,T>& operator/=(T);
 
-        util::vector<S>  operator+ (double) const;
-        util::vector<S>& operator+=(double);
+        vector<S,T>  operator+ (T) const;
+        vector<S,T>& operator+=(T);
 
-        util::vector<S>  operator- (double) const;
-        util::vector<S>& operator-=(double);
+        vector<S,T>  operator- (T) const;
+        vector<S,T>& operator-=(T);
 
+        // element operators
+        vector<S,T>  operator* (const vector<S,T>&) const;
+        vector<S,T>& operator*=(const vector<S,T>&);
 
-        util::vector<S>  operator* (const util::vector<S>&) const;
-        util::vector<S>& operator*=(const util::vector<S>&);
+        vector<S,T>  operator+ (const vector<S,T>&) const;
+        vector<S,T>& operator+=(const vector<S,T>&);
 
-        util::vector<S>  operator+ (const util::vector<S>&) const;
-        util::vector<S>& operator+=(const util::vector<S>&);
+        vector<S,T>  operator- (void) const;
+        vector<S,T>  operator- (const vector<S,T>&) const;
+        vector<S,T>& operator-=(const vector<S,T>&);
 
-        util::vector<S>  operator- (void) const;
-        util::vector<S>  operator- (const util::vector<S>&) const;
-        util::vector<S>& operator-=(const util::vector<S>&);
+        vector<S, T>& operator =(const vector<S,T>&);
 
-        util::vector<S>& operator =(const util::vector <S>&);
-
-        bool operator== (const util::vector <S>&) const;
-
-        double magnitude  (void) const;
-        double magnitude2 (void) const;
-
-        double dot (const util::vector<S>&) const;
-
-        util::vector<S>& normalise  (void);
-        util::vector<S>  normalised [[gnu::warn_unused_result]] (void) const;
+        // logical operators
+        bool operator== (const vector<S,T>&) const;
 
         bool is_zero (void) const;
 
-        template <size_t D> vector<D> redim (void) const;
+        // vector operators
+        T magnitude  (void) const;
+        T magnitude2 (void) const;
+
+        T dot (const vector<S,T>&) const;
+
+        vector<S,T>& normalise  (void);
+        vector<S,T>  normalised [[gnu::warn_unused_result]] (void) const;
+
+        template <size_t D> vector<D,T> redim (void) const;
 
         void sanity (void) const;
     };
 
-    vector<2> polar_to_cartesian (const vector<2>&);
+    // free vector operators
+    template <typename T> vector<2,T> polar_to_cartesian (const vector<2,T>&);
 
-    vector<3> cross (const vector<3>&, const vector<3>&);
-    vector<3> spherical_to_cartesian (const util::vector <3>&);
-    vector<3> cartesian_to_spherical (const util::vector <3>&);
+    template <typename T> vector<3,T> cross (const vector<3,T>&, const vector<3,T>&);
+    template <typename T> vector<3,T> spherical_to_cartesian (const vector<3,T>&);
+    template <typename T> vector<3,T> cartesian_to_spherical (const vector<3,T>&);
 
-    typedef vector<2> vector2;
-    typedef vector<3> vector3;
+    template <size_t S, typename T> vector<S,T> operator* (T, const vector<S,T>&);
+    template <size_t S, typename T> vector<S,T> operator+ (T, const vector<S,T>&);
+    template <size_t S, typename T> vector<S,T> operator- (T, const vector<S,T>&);
 
-    template <size_t S> util::vector<S> operator* (double, const util::vector<S>&);
-    template <size_t S> util::vector<S> operator+ (double, const util::vector<S>&);
-    template <size_t S> util::vector<S> operator- (double, const util::vector<S>&);
+    // output and serialisation operators
+    template <size_t S, typename T> std::ostream& operator<< (std::ostream&, const vector<S,T>&);
 
-    template <size_t S> std::ostream& operator<< (std::ostream&, const util::vector<S>&);
+    template <size_t S, typename T>
+    const json::node& operator>> (const json::node&, vector<S,T>&);
 
-    template <size_t S>
-    const json::node& operator>> (const json::node&, util::vector<S>&);
+    // convenience typedefs
+    typedef vector<2,float> vector2f;
+    typedef vector<3,float> vector3f;
+
+    typedef vector<2,double> vector2d;
+    typedef vector<3,double> vector3d;
 }