2D/3D Vector Classes

Most of my projects end up using vector mathematics, so I've written a couple of classes for handling the 2D and 3D cases. They have evolved a bit over time but have now reached a point where they do most everything I need.

The classes are templatized on the component type, so that I can support vectors of float and vectors of double with the same piece of code. I even use vectors of int sometimes.

There are two different notations for referring to the components of two- and three-dimensional vectors. Sometimes people use the letters x, y, and z since the coordinate system axes are traditionally named the same way. Sometimes people use numeric subscripts instead.

Numeric subscripts are superior to letter subscripts in at least one situation: when you need to choose components programmatically. For instance, let's say you've got a plane in 3D space, and you want to project that plane to 2D (to triangulate a polygon in the plane, say). An easy way to do this is to determine which of the three components of the plane's normal has the largest absolute value, and then use the axes corresponding to the other two components as the 2D projection plane. This type of projection is cheap since it just involves selecting two of the three components of a 3D point; a more general projection would involve dot products against basis vectors, which is computationally more expensive.

Because of the scenario above, I've gone with numeric subscripts in my vector classes. Some vector classes use a union of an array and members named x, y, and z to support both notations; I didn't bother.

Functions that operate on vectors can be standalone, or defined as methods of the vector class. There isn't too much difference either way. I implemented the dot, cross, len, etc. functions as methods of the class because the a.dot(b) notation puts the dot between the operands, like you would do in written math. It also keeps those names out of the global namespace; they are inside the vector class.

The one situation where you can't use methods is when you have an operator that doesn't take a vector as its first argument, for instance if you want to be able to multiply a scalar times a vector with the scalar on the left. This type of function has to be standalone.

The < operator gives vectors an ordering, which allows them to be used as keys in a lookup data structure such as std::map.

The empty constructor allows an uninitialized vector to be constructed, which is nice when you have an if-statement that initializes a vector in one of two different ways.

The +=, *= and similar operators which mutate a vector all return a reference to the vector itself, mimicking the behavior of these operators for standard types such as int. It's not typical or recommended, but you can write a = b += 5;. The vector class supports this behavior in order to be as much like a standard numeric type as possible.

This code is completely standalone, with one exception: square root finding. I use the standard C library functions sqrt and sqrtf for finding square roots of 64- and 32-bit floating point values, respectively. In order to choose the right square-root function for the scalar type at hand, I use template specialization.

These vector classes don't do anything to make use of the SIMD instruction sets that most processors have these days. Using special machine instructions from within C++ is platform-specific, and I have wanted to keep my code simple. Using SIMD badly can result in code that is slower than plain floating-point, too.

Here's the code:

#ifndef VEC_H
#define VEC_H

#include <cmath>

// Square root template, specialized for float or double.

template<typename T> T sqrt_template(T);
template<> inline float sqrt_template<float>(float value)
{
    return std::sqrtf(value);
}
template<> inline double sqrt_template<double>(double value)
{
    return std::sqrt(value);
}

// A basic 2D vector.
// vec2 = float version
// dvec2 = double version

template<typename T> class vec2_template
{
public:
    vec2_template() {}
    vec2_template(T x, T y)
    {
        m_v[0] = x;
        m_v[1] = y;
    }

    bool operator < (const vec2_template & rhs) const
    {
        if (m_v[0] < rhs.m_v[0]) return true;
        if (rhs.m_v[0] < m_v[0]) return false;
        return m_v[1] < rhs.m_v[1];
    }

    bool operator == (const vec2_template & rhs) const
    {
        return m_v[0] == rhs.m_v[0] && m_v[1] == rhs.m_v[1];
    }

    bool operator != (const vec2_template & rhs) const
    {
        return m_v[0] != rhs.m_v[0] || m_v[1] != rhs.m_v[1];
    }

    T & operator [] (int dim)
    {
        return m_v[dim];
    }

    const T & operator [] (int dim) const
    {
        return m_v[dim];
    }

    vec2_template operator + (const vec2_template & rhs) const
    {
        return vec2_template(
            m_v[0] + rhs.m_v[0],
            m_v[1] + rhs.m_v[1]);
    }

    vec2_template operator - (const vec2_template & rhs) const
    {
        return vec2_template(
            m_v[0] - rhs.m_v[0],
            m_v[1] - rhs.m_v[1]);
    }

    vec2_template operator / (T f) const
    {
        return vec2_template(m_v[0] / f, m_v[1] / f);
    }

    vec2_template & operator += (const vec2_template & rhs)
    {
        m_v[0] += rhs.m_v[0];
        m_v[1] += rhs.m_v[1];
        return *this;
    }

    vec2_template & operator -= (const vec2_template & rhs)
    {
        m_v[0] -= rhs.m_v[0];
        m_v[1] -= rhs.m_v[1];
        return *this;
    }

    vec2_template & operator *= (T f)
    {
        m_v[0] *= f;
        m_v[1] *= f;
        return *this;
    }

    vec2_template & operator /= (T f)
    {
        m_v[0] /= f;
        m_v[1] /= f;
        return *this;
    }

    vec2_template operator - () const
    {
        return vec2_template(-m_v[0], -m_v[1]);
    }

    vec2_template operator * (T f) const
    {
        return vec2_template(m_v[0]*f, m_v[1]*f);
    }

    T dot(const vec2_template & rhs) const
    {
        return m_v[0]*rhs.m_v[0] + m_v[1]*rhs.m_v[1];
    }

    // "Perp dot" product
    T perpdot(const vec2_template & rhs) const
    {
        return m_v[0]*rhs.m_v[1] - m_v[1]*rhs.m_v[0];
    }

    T sqlen() const
    {
        return dot(*this);
    }

    T len() const
    {
        return sqrt_template<T>(sqlen());
    }

    vec2_template normalized() const
    {
        return *this / len();
    }

    void normalize()
    {
        *this /= len();
    }

private:
    T m_v[2];
};

typedef vec2_template<float> vec2;
typedef vec2_template<double> dvec2;

// A basic 3D vector.
// vec3 = float version
// dvec3 = double version

template<typename T> class vec3_template
{
public:
    vec3_template() {}
    vec3_template(T x, T y, T z)
    {
        m_v[0] = x;
        m_v[1] = y;
        m_v[2] = z;
    }

    bool operator < (const vec3_template & rhs) const
    {
        if (m_v[0] < rhs.m_v[0]) return true;
        if (rhs.m_v[0] < m_v[0]) return false;
        if (m_v[1] < rhs.m_v[1]) return true;
        if (rhs.m_v[1] < m_v[1]) return false;
        return m_v[2] < rhs.m_v[2];
    }

    bool operator == (const vec3_template & rhs) const
    {
        return
            m_v[0] == rhs.m_v[0] &&
            m_v[1] == rhs.m_v[1] &&
            m_v[2] == rhs.m_v[2];
    }

    bool operator != (const vec3_template & rhs) const
    {
        return
            m_v[0] != rhs.m_v[0] ||
            m_v[1] != rhs.m_v[1] ||
            m_v[2] != rhs.m_v[2];
    }

    T & operator [] (int dim)
    {
        return m_v[dim];
    }

    const T & operator [] (int dim) const
    {
        return m_v[dim];
    }

    vec3_template operator + (const vec3_template & rhs) const
    {
        return vec3_template(
            m_v[0] + rhs.m_v[0],
            m_v[1] + rhs.m_v[1],
            m_v[2] + rhs.m_v[2]);
    }

    vec3_template operator - (const vec3_template & rhs) const
    {
        return vec3_template(
            m_v[0] - rhs.m_v[0],
            m_v[1] - rhs.m_v[1],
            m_v[2] - rhs.m_v[2]);
    }

    vec3_template operator / (T f) const
    {
        return vec3_template(m_v[0] / f, m_v[1] / f, m_v[2] / f);
    }

    vec3_template & operator += (const vec3_template & rhs)
    {
        m_v[0] += rhs.m_v[0];
        m_v[1] += rhs.m_v[1];
        m_v[2] += rhs.m_v[2];
        return *this;
    }

    vec3_template & operator -= (const vec3_template & rhs)
    {
        m_v[0] -= rhs.m_v[0];
        m_v[1] -= rhs.m_v[1];
        m_v[2] -= rhs.m_v[2];
        return *this;
    }

    vec3_template & operator *= (T f)
    {
        m_v[0] *= f;
        m_v[1] *= f;
        m_v[2] *= f;
        return *this;
    }

    vec3_template & operator /= (T f)
    {
        m_v[0] /= f;
        m_v[1] /= f;
        m_v[2] /= f;
        return *this;
    }

    vec3_template operator - () const
    {
        return vec3_template(-m_v[0], -m_v[1], -m_v[2]);
    }

    vec3_template operator * (T f) const
    {
        return vec3_template(m_v[0]*f, m_v[1]*f, m_v[2]*f);
    }

    T dot(const vec3_template & rhs) const
    {
        return
            m_v[0]*rhs.m_v[0] +
            m_v[1]*rhs.m_v[1] +
            m_v[2]*rhs.m_v[2];
    }

    vec3_template cross(const vec3_template & rhs) const
    {
        return vec3_template(
            m_v[1]*rhs.m_v[2] - m_v[2]*rhs.m_v[1],
            m_v[2]*rhs.m_v[0] - m_v[0]*rhs.m_v[2],
            m_v[0]*rhs.m_v[1] - m_v[1]*rhs.m_v[0]);
    }

    T sqlen() const
    {
        return dot(*this);
    }

    T len() const
    {
        return sqrt_template<T>(sqlen());
    }

    vec3_template normalized() const
    {
        return *this / len();
    }

    void normalize()
    {
        *this /= len();
    }

private:
    T m_v[3];
};

typedef vec3_template<float> vec3;
typedef vec3_template<double> dvec3;

#endif

At the moment I am taking another crack at learning how to implement rigid body dynamics. This is a fairly well-trod path. I'm working from the highly-influential tutorials by Andrew Witkin and David Baraff. In my case, I'm developing it in 2D, not necessarily because it's any simpler, but because that's what I need it for. I hope to begin writing about it next week.