The Compiler as a Physics Engine

Using strong types for zero-overhead dimensional analysis at compile time

In our first two posts in this series, we discussed the need for strong types in C++ to improve type safety and better document APIs, all while adding zero runtime cost to our software beyond using simple primitives. We first defined these strong types, then added traits to describe the permitted operations on them using CRTP. From here though, we can hook further into C++’s compile time system to perform changes on the strong types as we compute, permitting us to write code that changes types, preventing incompatible type operations and generating new types on the fly.

Recall back to our first example, to compute the new position of some object after an elapsed period of time, using just our initial strong types

float determine_future_position(
  TimePoint t0, Position pos0, Velocity vel0, Acceleration acc0, TimePoint t1)
  {
    const auto time_delta = t1.value() - t0.value();
    return pos0.value() + vel0.value() * time_delta + 0.5f * acc0.value() * time_delta * time_delta;
  }

With our CRTP strong types, we were able to define Addable, Subtractable and Multiplyable traits on each of these types, but without further changes, they cannot be inter-operated. For example, a strong type Velocity cannot be multiplied by a strong type Time since they are different strong types. Clearly this is one of the basic goals of first creating strong types, making it so we can’t mix operations on different types. Additionally, what would the type of the output of (Velocity * Time) be? It shouldn’t be either velocity or time, it should be a new type of Distance as described by physics. Let’s now work through how we can implement the creation of this new type and enable our C++ compiler to understand physics computations.

Physical Laws to Templated Types

In order to represent our determine_future_position function with physically correct computations, we first need to take a step back and consider in the real world what each of these types actually represents. We’ll be using the SI units m to represent meters, s to represent seconds. As such we have

\(\text{distance} = m\)

\(\text{velocity} = m*s^{-1}\)

\(\text{acceleration} = m*s^{-2}\)

With this representation, we see our “types” for distance, velocity and acceleration actually have 2 dimensions that need to be represented, time and distance. Since these types are differentiated by the exponents of their units, conceptually we need to encode both the units and their exponents into our corresponding C++ types. Using the ability to template on a non-type parameter such as int, we can define our initial representation to be a templated type Dimension with template parameters to represent the exponent of both Meters and Seconds¹

template <int Meter, int Second>
struct Dimension
{
};

This allows us to create our above types as aliases

using DistanceDim = Dimension<1, 0>; // Meters^1, Seconds^0
using VelocityDim = Dimension<1, -1>; // Meters^1, Seconds^-1
using AccelerationDim = Dimension<1, -2>; // Meters^1, Seconds^-2
using TimeDim = Dimension<0, 1>; // Meters^0, Seconds^1
using ScalarDim = Dimension<0, 0>; // Dimensionless scalar quantity

Going back to our strong types, since those were templated on a type T for their tag, we can now template our strong types on a Dimension as their tag. Additionally, we will make our types Addable and Subtractable since we should support operations such as adding together two velocities, or subtracting two distances. With some additional type aliases, we now have the following

// Convenience alias for any dimensional type providing the
// underlying type (float), and the Addable/Subtractable traits
template <typename DimT>
using quantity = strong_type<float, DimT, Addable, Subtractable;

using Distance = quantity<DistanceDim>;
using Velocity = quantity<VelocityDim>;
using Acceleration = quantity<AccelerationDim>;
using Time = quantity<TimeDim>;
using Scalar = quantity<ScalarDim>;

We can now express our addition and subtraction on a single physics type cleanly without dropping down to the underlying type, but still don’t have the ability to do cross-type operations. Since we want our types to be able to be multiplied and divided, but not return the same type, we’re going to define our operator* and operator/ externally, simply performing multiplication and division on our underlying values.

template <typename DimLHS, typename DimRHS>
quantity<???> operator*(const quantity<DimLHS> lhs, const quantity<DimRHS> rhs)
{
  return quantity<???>(lhs.value() * rhs.value());
}

template <typename DimLHS, typename DimRHS>
quantity<???> operator/(const quantity<DimLHS> lhs, const quantity<DimRHS> rhs)
{
  return quantity<???>(lhs.value() / rhs.value());
}

These functions perform our multiply and divide operations on two arbitrary physics types, but I’ve explicitly left the resulting type as some unknown quantity type. To determine what that resulting type is going to be, we’re going to need to do some compile-time math.

Compile Time Math

We’re going to first extend our Dimension type we created above to allow easier access to the non-type template parameters in expressions, and then define a struct that exposes a type alias for when two dimensions are multiplied or divided

template <int Meter, int Second>
struct Dimension
{
  constexpr static auto meter_exponent = Meter;
  constexpr static auto second_exponent = Second;
};

// Helper structs to determine the type of multiply/divide of two different dimensional types
template <typename DimLHS, typename DimRHS>
struct DimMultiply
{
  // Add exponents together to represent multiplication
  // eg m^1 * m^1 = m^2 
  using type = Dimension<DimLHS::meter_exponent + DimRHS::meter_exponent,
                         DimLHS::second_exponent + DimRHS::second_exponent>;
};

template <typename DimLHS, typename DimRHS>
struct DimDivide
{
  // Subtract exponents to represent division
  using type = Dimension<DimLHS::meter_exponent - DimRHS::meter_exponent,
                         DimLHS::second_exponent - DimRHS::second_exponent>;
};

With these new helper structs, we can validate at compile time that our dimensional operations work as expected

using maybe_velocity = DimDivide<DistanceDim, TimeDim>::type;
static_assert(std::is_same_v<maybe_velocity, VelocityDim>);
using maybe_also_velocity = DimMultiply<AccelerationDim, TimeDim>::type;
static_assert(std::is_same_v<maybe_also_velocity, VelocityDim>);
// Multiply by a scalar leaves the type unchanged
using maybe_position = DimMultiply<DistanceDim, ScalarDim>::type;
static_assert(std::is_same_v<maybe_position, DistanceDim>);

The static asserts all pass indicating that our dimensional multiplication and division is indeed creating the correct types.

Putting it all together

Now all we need to do is to plug that into our multiplication and division functions for the tag on the strong type and we should be good to go.

template <typename DimLHS, typename DimRHS>
auto operator*(quantity<DimLHS> lhs, const quantity<DimRHS> rhs)
{
  return quantity<typename DimMultiply<DimLHS, DimRHS>::type>(lhs.value() * rhs.value());
}

template <typename DimLHS, typename DimRHS>
auto operator/(quantity<DimLHS> lhs, const quantity<DimRHS> rhs)
{
  return quantity<typename DimDivide<DimLHS, DimRHS>::type>(lhs.value() / rhs.value());
}

Since we’ve defined our dimensional multiplication and division, we are now able to fully write out our determine future position function without ever using the underlying types

Distance determine_future_position(
  Time t0, Distance pos0, Velocity vel0, Acceleration acc0, Time t1)
  {
    const auto time_delta = t1 - t0;
    return pos0 + vel0 * time_delta + Scalar{0.5f} * acc0 * time_delta * time_delta;
  }

By writing our multiplication and division operators to return new types, we’re able to fully express our mathematical equation without resorting to using primitive operations. Of course, we still have the exact same assembly output for this function as when we used primitive types, since all the code to generate the new types is executed at compile time.

.LCPI0_0:
        .long   0x3f000000
determine_future_position(strong_type<float, Dimension<0, 1>, Addable, Subtractable>, strong_type<float, Dimension<1, 0>, Addable, Subtractable>, strong_type<float, Dimension<1, -1>, Addable, Subtractable>, strong_type<float, Dimension<1, -2>, Addable, Subtractable>, strong_type<float, Dimension<0, 1>, Addable, Subtractable>):
        subss   xmm4, xmm0
        mulss   xmm2, xmm4
        addss   xmm2, xmm1
        mulss   xmm3, dword ptr [rip + .LCPI0_0]
        mulss   xmm3, xmm4
        mulss   xmm3, xmm4
        addss   xmm3, xmm2
        movaps  xmm0, xmm3
        ret

Furthermore, the compiler now checks us to ensure we’re not performing invalid operations. For example, if we made a mistake in our function, such as the following

Distance determine_future_position(
  Time t0, Distance pos0, Velocity vel0, Acceleration acc0, Time t1)
  {
    const auto time_delta = t1 - t0;
    return pos0 + vel0 /*forgot to multiply velocity by time delta*/ + 
                  Scalar{0.5f} * acc0 * time_delta * time_delta;
  }

We get a fantastic error message telling us explicitly that we cannot add ‘Distance’ and ‘Velocity’.

Beyond Written Types

One of the other powerful features of this technique is that we can create new types for the compiler without ever having to write them out in code. We never defined ‘Jerk’ (the first derivative of acceleration) as a concept, but if we run the following code, we can have the compiler generate jerk as a type.

auto calculate_jerk(Acceleration a0, Time t0, Acceleration a1, Time t1)
{
  return (a1 - a0) / (t1 - t0);
}

Acceleration a0{0.1}, a1{0.2};
Time t0{0}, t1{1};
auto jerk = calculate_jerk(a0, t0, a1, t1);
// jerk is of type quantity<Dimension<1, -3>>

With our dimensional types, we’ve showed how we can encode physical concepts into our C++ compiler, allowing us to further utilize strong types and generate safer code with zero overhead. This can also be extended into many other domains, such as encoding options greek exposures with concepts like delta and theta being encoded as dimensional quantities with respect to the underlying or time respectively.

A huge thank you to anyone who read through this first series of articles detailing strong types and their extensions. Future posts will likely cover more performance related concepts as well as software engineering principles. If you’ve made it this far, please consider subscribing to stay up to date on all new posts. Thank you for reading, until next time!

1. Note we could use the full set of SI units and include mass, current, luminosity, temperature and substance, but we’re keeping it simple for this article. It would be straightforward to add these additional units following the same framework laid out here. Please also check out Nic Holthaus’ fantastic units library which utilizes and expands far beyond what this post describes