Blog

Do Your Units Match? An Embedded DSL for Physical Quantities

Aug 18, 2023 | 32 minutes read

Tags: pl, cpp

Computations of physical values in programming languages can quickly get out of hand due to the complex units that the programmer has to keep track of. Good variable naming and comments can help, but nothing is stopping double run_time_sec from storing a time measured in milliseconds. Furthermore, different functions might use the same value, but measured in different units, requiring manual conversions. Overall, it’s just an unpleasant burden placed on the programmer and a cause of many errors. For example, at an internship, due to data values being passed through multiple standard data structures, each requiring values in different units, a value originally measured in knots would have to be converted from knots (nautical miles per hour) to nautical miles per second, then to meters per second, and finally to feet per second. Suffice it to say, it was confusing keeping track of which structure required which units and making sure the correct conversions were occurring at the correct location.

One solution to this problem is Rust’s Newtype Pattern.

The gist of this idiom is creating structs to wrap (typically primitive) data values so that the programmer cannot accidentally use a value of one type where another is expected. In other words, it uses the type system to encode values with different semantics as different types. So instead of values measured in miles and meters both represented as f64, we can represent them with different structs that contain a f64.

For example:

struct Meters(f64);
struct Miles(f64);

fn foo(a: &Miles) {
    println!("{} miles", a.0);
}

fn main() {
    let miles = Miles(10);
    let meters = Meters(20);
    foo(miles);
    foo(meters); // Error
    foo(20.0); // Error
}

However, this pattern doesn’t enforce correct unit conversions:


fn speed(distance: Miles, time: Duration) -> MetersPerSecond {
    // Incorrect conversion from miles to meters !
    MetersPerSecond(distance.0 * 100 / time.as_secs_f64())
}

What if we could also encode conversions in the type system? So not only would the compiler be able to catch unit mismatch errors, but they would also catch conversion errors too? Going one step further, even values with the same unit might encode different information. For example, width, height, and depth can all be represented as meters, but we might like the compiler to prevent us from passing a width * width to a function expecting an area (width * height).

I explored all of these questions when developing a small embedded DSL (Domain Specific Language) for C++17 last year.

Let’s first define a few terms I’ll use:

Mismatch error - Error from using a value measured in one unit where a value measured in another was expected
Conversion error - Error from improperly converting units

Besides units, the system needed to handle different scalings (think kilometer, millimeter, etc. which are all just scalings of meters) and different semantic types within a given unit (such as the height and width example)

There were a few criteria I had when developing this:

I wanted the system to be automatic. A programmer would simply define a unit and the compiler would handle the rest.
It needed to catch both mismatch and conversion errors at compile time for units (meters vs. seconds), scalings (kilometers vs. centimeters), and semantic units (width vs. height).
Preferably, the system could handle scaling on its own. For example, passing a unit measured in kilometers to a function expecting a value in meters would work, and be correct with the compiler automatically scaling the value appropriately. This requirement comes from the way the time units in std::chrono work.
Be relatively easy on the user

The solution I came up with was to use extensive C++ template metaprogramming and encode the entire unit type system addon as template types.

The core idea is to encode a unit as a type list of powers of unit types. So something like $\frac{m^3}{s}$ could be encoded as

Pack<PowerType<Meters, 3, 1>, PowerType<Seconds, -1, 1>>

Where PowerType is defined as so:

/// A unit raised to a rational power
template<typename Unit, int32_t numerator, int32_t denominator>
struct PowerType;

Then, using template metaprogramming, we can implement transformations on these lists of power types to represent different computations such as multiplying, dividing, adding, subtracting, and raising values to a power.

By encoding units using template parameters of a template type, we get the compiler enforcement against mismatch and conversion errors “for free” because two different units would yield two different instantiations of a template which are considered different types by the type checker. That is, different units are represented as different (not implicitly convertible) types and incorrect conversions will not produce the same type.

Template Metaprogramming

This post will not walk through C++ template metaprogramming. But I will take the time to point out that this feature essentially relies on template specialization. The core idea is that we have a struct specialization as a base case and the struct definition as a main case. TMP in C++ does not have loops or variables, so our implementations will be purely functional.

Let’s look at an example of something we’ll need later: a simple conditional. Our inputs will be a constexpr boolean expression, a TrueType and a FalseType. If the condition is true, the TrueType will be “returned”, otherwise the FalseType will be “returned”. To “return” a type, we typically use a type alias named Type, which the user can use to “query” the result.

// False case
template<bool expr, typename TrueType, typename FalseType>
struct TypeConditional {
    // use of type aliases for "variables" and "returns"
    using Type = FalseType;
};

// Specialization for true case
template<typename TrueType, typename FalseType>
struct TypeConditional<true, TrueType, FalseType> {
    using Type = TrueType;
};

And then we can use this like so:

using i32 = typename TypeConditional<sizeof(int) == 4, int, long>::Type;

However, I typically like to create template type aliases to reduce all the boilerplate usages of Type and typename:

template<bool expr, typename TrueType, typename FalseType>
/// TrueType if expr is true, FalseType otherwise
using type_cond_t = typename TypeConditional<expr, TrueType, FalseType>::Type;

So our earlier example would become:

using i32 = type_cond_t<sizeof(int) == 4, int, long>;

But we can do a lot more than simple conditionals; we can implement arbitrary recursive functions! We do this by employing a pattern where we have a struct with a recursive definition and a template specialization to serve as a base case.

For example, we can compute the factorial like so:

// recursive case
template <unsigned v>
struct Factorial {
    constexpr static value = v * Factorial<v - 1>::value;
};

// base case
template <>
struct Factorial<0> {
    constexpr static value = 1;
}

static_assert(Factorial<10>::value == 3628800);

In C++17, constexpr functions can replace template metaprogramming that computes values such as this Factorial example, however, it cannot replace algorithms that “compute” a type, which is what we’ll focus on for the rest of this post.

Another core structure we’ll use is the type list. Traditionally, this can be represented as a linked list of template types. For example:

template <typename Head, typename Tail>
struct List;

List<int, List<bool, List<char, EndType>>> f;

In my implementation, I chose to use parameter packs for my type lists.

/// A type list
/// This type list is represented just as a parameter pack
/// The type list does not store a tail node, but operations on the type list
/// will return `Pack<EmptyPack>` when there are no more elements
template<typename ... Ts>
struct Pack {};

// usage: Pack<int, short, long, std::string>

In the following discussion, I will use “type list” and “pack” fairly interchangeably. I use the term “pack” since it is not a traditional linked list and because it relies on C++ parameter packs.

To represent an empty type list, I used a Pack consisting only of a special EmptyPack type. While this enabled a less bulky notation for defining and working with type lists, it required implementing special cases to handle empty lists.

The Unit

We define our unit type as so:

using scale_t = Rational;


/// @brief A semantically typed unit class
/// @tparam T The underlying type of the unit
/// @tparam scale the scale of the unit
/// @tparam Semantic the Pack of subcategory PowerTypes of the unit to 
///     differentiate different values of the same unit. 
///     Ex: `Pack<PowerType<Length, 2, 1>>` which is Length^2
///     Requires every type in the Pack are subtypes of UnitBase
/// @tparam UnitPowerPack the pack of PowerTypes of Units that the value has
///     Ex. `Pack<PowerType<Meters, 2, 1>>`
///     Requires every type in the Pack are subtypes of UnitBase
template<typename T, scale_t scale, typename SemanticPowerPack, typename UnitPowerPack>
struct Unit;

Each unit will have a base numeric representation such as double or int. Then each unit will have a scale, which is how we can differentiate something like kilometers from meters. The scale will be a rational number. The Rational class has all constexpr members so it can be used in a constexpr context, and its constructor will ensure that all rational numbers are represented in simplest form with a negative numerator for negative rationals. Units will have attached semantic types which provide further semantic information such as width and height. This is encoded in SemanticPowerPack, which is a Pack<> of PowerTypes of semantic unit types (like width and height). For example, we could have: Pack<PowerType<Width, 1, 1>, PowerType<Length, 1, 1>>. Finally, each unit will have a UnitPowerPack which is the Pack of PowerTypes of regular unit types (like meters and seconds). For example, we could have: Pack<PowerType<Seconds, 1, 1>, PowerType<Meters, 2, 1>>.

You may be wondering what is the difference between unit types (the UnitPowerPack) and SemanticUnitTypes (the SemanticPowerPack). From an implementation standpoint, they will behave the same way. However, the semantic types provide a subcategory within the unit. The intended usage is for a unit to match up with physical units in the real world such as meters or miles. However, even among values measured in the same unit, they could represent different things and not be interchangeable. Thus, the semantic unit type provides more fine-grained type checking to prevent mismatch errors between values with the same unit that measures different quantities.

In this discussion, I will refer to types that are part of the SemanticPowerPack as “semantic unit types”. These include things like Length and Width from my previous example. I will refer to types that are part of the UnitPowerPack as “unit types”. These include Seconds and Meters. I will refer to the Unit class as Unit and refer to a unit in the general sense as “unit”. Finally, I will refer to a unit or semantic unit type, combined with an exponent in a PowerType, as a “power type”.

As an example, an area measured in meters squared might be encoded in this system as

using area_t = Unit<double, Rational{1}, 
    // Semantic Power Pack
    sort_unit_pack_t<Pack<PowerType<Length, 1, 1>, PowerType<Width, 1, 1>>>,
    // Unit Power Pack
    Pack<PowerType<Meters, 2, 1>>>;

This is a tad unwieldy of notation, so let’s break it down. The base numeric type is double, and its scale factor is 1. Then, semantically, an area is composed of a length value times a width value. This is denoted by Pack<PowerType<Length, 1, 1>, PowerType<Width, 1, 1>>, however, we sort this list to ensure that each Unit has a unique representation. In terms of units, an area is measured in meters squared, denoted by Pack<PowerType<Meters, 2, 1>>.

Some Units do not have a semantic unit time. A user might only have one notion of time, for example. Therefore, SemanticPowerPack can be empty (Pack<EmptyPack>).

We’d like to define operations on Unit such as +, -, /, *, and pow that “do the right thing”. For example, multiplication a * b should return a value of type Unit with a unit power pack that is the concatenation of the unit types of a and b. Likewise, the addition should only work on Units of the same type.

We can start by defining addition and subtraction fairly easily. The operands must have the same unit and semantic types (or one of the semantic types must be empty). Then, we simply add the values taking the scale factor into account:

/// True if `Self` and `Other` are convertible semantic power packs 
template<typename Self, typename Other>
constexpr bool is_semantic_convertable_v = std::is_same_v<Self, Other> || 
                                           std::is_same_v<Self, Pack<EmptyPack>> ||
                                           std::is_same_v<Other, Pack<EmptyPack>>;


template<typename T, scale_t scale, typename Semantic, typename Units,
    scale_t otherScale, typename OtherSemantic>
constexpr auto operator+(Unit<T, scale, Semantic, Units> a, 
                         const Unit<T, otherScale, OtherSemantic, Units>& b)
    -> std::enable_if_t<is_semantic_convertable_v<Semantic, OtherSemantic>, decltype(a)>
{
    a.val += static_cast<T>(b.val * static_cast<double>(otherScale / scale));
    return a;
}

template<typename T, scale_t scale, typename Semantic, typename Units,
    scale_t otherScale, typename OtherSemantic>
constexpr auto operator-(Unit<T, scale, Semantic, Units> a, 
                         const Unit<T, otherScale, OtherSemantic, Units>& b)
    -> std::enable_if_t<is_semantic_convertable_v<Semantic, OtherSemantic>, decltype(a)>
{
    a.val -= static_cast<T>(b.val * static_cast<double>(otherScale / scale));
    return a;
}

We use a trailing return type that will cause a compiler error if is_semantic_convertable_v<Semantic, OtherSemantic> evaluates to false, otherwise, it will return a Unit with the same type as the left operand.

Notice, we pass the first operand by value since we want to make a copy of it. However, the second operand we only read from so we can pass by const reference.

To provide implicit conversion between Units of different scales, we also overload operator=:

template<typename OtherSemantic, scale_t otherScale>
constexpr auto operator=(const Unit<T, otherScale, OtherSemantic, UnitPowerPack>& other)
    -> std::enable_if_t<is_semantic_convertable_v<SemanticPowerPack, OtherSemantic>, Unit&> 
{
    val = static_cast<T>(other.val * static_cast<double>(otherScale / scale));
    return *this;
}

To handle the other operations, we’ll need to develop a few utilities for manipulating type lists.

The problems we need to solve can be summed up as follows:

Type list concatenation: a: Meters * b: Seconds should result in c: Pack<Meters, Seconds>
Unique ordering: we want Meters * Seconds to correspond to the same type as Seconds * Meters
Type list removal: a: Pack<Meters, Seconds> / b: Pack<Seconds> should result in c: Pack<Meters>
Raising a power pack to a rational power and the uniqueness involved with this: we want PowerType<Meters, -1, 2> to be the same as PowerType<Meters, 1, -2> and we want to be able to recognize that
```
Pack<PowerType<Meters, 1, 1>, PowerType<Seconds, 1, 1>> * Pack<PowerType<Meters, 2, 1>>
```
should result in a value of the type Pack<PowerType<Meters, 3, 1>, PowerType<Seconds, 1, 1>>
Zero power units: $m^3s^0N^0$ should be the same as $s^0m^3$ and $m^3$

Let’s go through the solutions I developed for each one:

Type List Manipulations

Let’s start with the first problem: type list concatenation. For inspiration, let’s write a simple list concatenation function in OCaml:

(** [concat a b] is the list [a] followed by the list [b] *)
let rec concat a b =
    match a with
    | hd :: tl -> hd :: (concat tl b)
    | [] -> b;

We’re going to translate this same thing into C++ TMP. First, we need to destruct a type list into its head and tail:

/// Gets the head of a type list. Returns `Pack<EmptyPack>` if the list is empty
template<typename Head, typename ... Tail>
constexpr auto pack_head(Pack<Head, Tail...> p) {
    return Head{};
}

constexpr auto pack_head(Pack<EmptyPack> p) {
    return EmptyPack{};
}

/// Gets the tail of the type list. Returns `Pack<EmptyPack>` if there is no tail
template<typename Head, typename ... Tail>
constexpr auto pack_tail(Pack<Head, Tail...> p, float)
    -> std::enable_if_t<sizeof...(Tail) != 0, Pack<Tail...>> 
{
    return Pack<Tail...>{};
}

template<typename T>
constexpr auto pack_tail(Pack<T> p, int) {
    return Pack<EmptyPack>{};
}

// Usage type aliases:

template<typename Pack>
using pack_head_t = decltype(pack_head(std::declval<Pack>()));

template<typename Pack>
using pack_tail_t = decltype(pack_tail(std::declval<Pack>(), 0));

We use constexpr functions to destruct Pack and get out the parameter pack that it encapsulates. We use std::declval to create a value of type Pack, and decltype to get the return type of the constexpr function call.

The unused float and int parameters on pack_tail are for SFINAE purposes. Essentially, it will first try and instantiate pack_tail(Pack<Head, Tail...>, float) but if that fails, it will revert to pack_tail<Pack<T>, int).

Next, we need to write a cons helper template to cons a type onto a type list. This one is more straightforward:

/// The return type is a Pack of T consed onto Ts
/// Consing an element of type T onto a Pack<EmptyPack> will result in Pack<T>
/// @{
template<typename T, typename ... Ts>
constexpr auto pack_cons(T a, Pack<Ts...> p) {
    return Pack<T, Ts...>{};
}

template<typename T>
constexpr auto pack_cons(T a, Pack<EmptyPack> p) {
    return Pack<T>{};
}
/// @}

template<typename T, typename Pack>
/// Gets the Pack of T consed onto Pack
/// Consing an element of type T onto a Pack<EmptyPack> will result in Pack<T>
using pack_cons_t = decltype(pack_cons(std::declval<T>(), std::declval<Pack>()));

With these out of the way, we can finally implement our recursive concatenation function:

/// Returns a pack of types of the types in PackA followed by the types in PackB
/// @{
template<typename PackA, typename PackB>
struct PackConcat;

// base case for when we have `[] @ pack`
template<typename PackB>
struct PackConcat<Pack<EmptyPack>, PackB> {
    using Type = PackB;
};

// base case for `pack @ []`
template<typename PackA>
struct PackConcat<PackA, Pack<EmptyPack>> {
    using Type = PackA;
};

// Base case for `[] @ []`. We need a special case for this to
// resolve ambiguity between the two other base cases
template<>
struct PackConcat<Pack<EmptyPack>, Pack<EmptyPack>> {
    using Type = Pack<EmptyPack>;
};

// recursive case: `hd :: (concat tl b)`
template<typename PackA, typename PackB>
struct PackConcat {
    using Type = pack_cons_t<pack_head_t<PackA>, typename PackConcat<pack_tail_t<PackA>, PackB>::Type>;
};


// user interface
template<typename PackA, typename PackB>
using pack_concat_t = typename PackConcat<PackA, PackB>::Type;
/// @}

Observe the redundant base cases. We need this to prevent ambiguity if we ever come across pack_concat_t<Pack<EmptyPack>, Pack<EmptyPack>>. In this case, the compiler wouldn’t know whether to choose the struct PackConcat<PackA, Pack<EmptyPack>> or the struct PackConcat<Pack<EmptyPack>, PackB> specialization since both are equally good matches. To avoid this, we provide a specific struct PackConcat<Pack<EmptyPack>, Pack<EmptyPack>> specialization. I will note that we won’t use pack_concat_t, but we’ll build upon the idea since combining exponents takes a little more logic than just concatenating power packs. This is nowhere near as pretty as our OCaml reference implementation but with this, we can do the following:

static_assert(std::is_same_v<pack_cons_t<Meters, Pack<Meters, Meters, Seconds, int>>, 
        Pack<Meters, Meters, Meters, Seconds, int>>);
static_assert(std::is_same_v<pack_concat_t<Pack<Meters, int>, Pack<Meters, Seconds>>,
         Pack<Meters, int, Meters, Seconds>>);

The next thing we want to do is to be able to determine when two unit-type lists represent the same thing. We can do this by comparing two sorted type lists. But we need some way to uniquely identify each type at compile time. To do this, every unit and semantic unit type will be a subtype of UnitBase. This class will provide a public static id member that we can use to order the types.

/// The base class for all unit types
/// Gives the inheriter a unique constexpr id
template<uint32_t major, uint32_t minor>
struct UnitBase {
    static constexpr uint64_t id = uint64_t(major) << 32 | minor;
};

We can then use a macro to give each unit type a unique ID based on the file and line number at which it is defined.

#define SEMANTIC_UNIT_TYPE UnitBase<file_hash(__FILE__), __LINE__>

struct Grams : SEMANTIC_UNIT_TYPE {};

In my demo implementation, I simply made file_hash a function that returns the length of the filename, but an actual implementation would want to use a good compile-time hash function.

With this, we can sort a type list based on the id of each type.

Again, let’s start with a reference implementation in OCaml. I chose to use Selection Sort:

(** [min x] is the minimum element of the list [x] *)
let rec min x = 
  match x with
  | hd :: tl -> 
      let tl_min = min tl in
      if hd < tl_min then hd else tl_min
  | [] -> max_int;

(** [remove lst x] is the list [lst] without the first occurence of the 
    element [x] *)
let rec remove lst x =
  match lst with
  | hd :: tl -> 
      if hd = x then tl else hd :: (remove tl x)
  | [] -> [];

(** [sort ls] is the list [ls] in sorted, least to greatest order *)
let rec sort ls = 
  match ls with 
  | [] -> []
  | lst ->
      let min = min lst in
      min :: sort (remove lst min);

We can test this out in UTOP to convince ourselves we haven’t made a mistake.

Let’s start translating:

We begin with the min helper function. I chose to abstract the comparator away by using a template-template parameter. This is a template parameter this is itself a template.

/// Comparator to order unit types
template<typename A, typename B>
struct UnitIdComparator {
    static constexpr bool value = A::id < B::id;
};

/// Gets the minimum type in the parameter pack as denoted by the comparator
/// @tparam Comparator a type with a `value` member that is `true` if its first type
///     parameter is less than its second type, `false` otherwise.
//      Comparator should establish a partial ordering
/// @tparam Head head of the parameter pack
/// @tparam Tail tail of the parameter pack
template<template <typename, typename> class Comparator, typename Head, typename ... Tail>
struct MinType {
private:
    // recursive call
    using TailMinType = typename MinType<Comparator, Tail...>::Type;
public:
    // if Head < TailMinType then Head else TailMinType
    using Type = type_cond_t<Comparator<Head, TailMinType>::value, Head, TailMinType>;
};

// base case for when the pack has one element
template<template <typename, typename> class Comparator, typename Last>
struct MinType<Comparator, Last> {
    using Type = Last;
};

template<typename ... Ts>
/// The minimum of all the unit types in the template list
/// Requires all types to have an id member, which is used to order them
using min_unit_type_t = typename MinType<UnitIdComparator, Ts...>::Type;

Note that MinType is operating on raw C++ parameter packs and not Pack<...>. We use constexpr functions to convert the two representations. Next up is the remove helper function:

/// Removes the first occurrence of T from the list
/// Returns a Pack of types without T
/// @{
template<typename T, typename Head, typename ... Tail>
struct RemoveType {
    // if hd = x then
    //  tl
    // else
    //  hd :: remove tl x
    using Type = type_cond_t<std::is_same_v<T, Head>, Pack<Tail...>, 
        pack_cons_t<Head, typename RemoveType<T, Tail...>::Type>>;
};

// base case
template<typename T, typename Head>
struct RemoveType<T, Head> {
    using Type = type_cond_t<std::is_same_v<T, Head>, Pack<EmptyPack>, Pack<Head>>;
    // note that pack_cons_t appropriately handles T :: Pack<EmptyPack> by
    // resulting in Pack<T> and NOT Pack<T, EmptyPack>
};
/// @}

/// Removes the first occurrence of T from the the Pack
/// @{
template<typename T, typename ... Ts>
constexpr auto remove_from_pack(T a, Pack<Ts...> p) {
    // constexpr function to extract Ts from Pack<Ts...>
    return typename RemoveType<T, Ts...>::Type{};
}

template<typename T, typename Pack>
using remove_pack_t = decltype(remove_from_pack(std::declval<T>(), std::declval<Pack>()));

And finally, the top-level sort “function”:

/// Returns a sorted pack of types, ordered by the comparator via Selection Sort
/// @tparam Comparator - a type with a `value` member that is `true` if its first type
///     parameter should be ordered before its second
/// @tparam ... Ts - parameter pack of types to sort
/// @returns a PACK of the ordered types. So this type constructs a sorted type list
///     from a parameter pack
/// @{
template<template<typename, typename> class Comparator, typename ... Ts>
struct SortTypes;

template<template<typename, typename> class Comparator, typename ... Ts>
constexpr auto sort_pack(Pack<Ts...> p) {
    return typename SortTypes<Comparator, Ts...>::Type{};
}

// destruct Pack<...>
template<template<typename, typename> class Comparator, typename Pack>
using sort_pack_t = decltype(sort_pack<Comparator>(std::declval<Pack>()));

// user "interface" for sorting a pack of unit or semantic unit types
template<typename Pack>
using sort_unit_pack_t = sort_pack_t<UnitIdComparator, Pack>;

// recursive case
template<template<typename, typename> class Comparator, typename ... List>
struct SortTypes {
private:
    // let min_elem = min lst
    using Min = typename MinType<Comparator, List...>::Type;

    // let tail_pack = remove lst min_elem
    using TailPack = typename RemoveType<Min, List...>::Type;
public:
    // min_elem :: tail_pack
    using Type = pack_cons_t<Min, sort_pack_t<Comparator, TailPack>>;
};

// base case
template<template<typename, typename> class Comparator, typename Last>
struct SortTypes<Comparator, Last> {
    using Type = Pack<Last>;
};

/// @}

Observe that SortTypes is effectively a higher-order function on types. It takes a list of types and a comparator that orders two types, and then it “returns” the sorted list of types.

Power List Manipulation

With these basic tools out of the way, we can build up to manipulating lists of PowerTypes, which is a part of the representation of Unit. Multiplication of power types is not as simple as concatenating lists, we want to concatenate lists and add the powers of repeated units. So we want $m^2 \cdot m^3$ to become $m^5$. If we exponentiate a power type by another power type, we want to do the same thing except multiply powers. The library currently doesn’t support this, but I abstracted away how two PowerTypes with the same base unit type are combined with a template-template parameter to make this feature easier to add later. This allows us to use the same procedure for adding exponents and multiplying exponents.

We define a new power type cons helper. Calling this cons is a poor naming choice as what it is doing is searching through the list and appending the type if it does not already exist. Otherwise, it replaces the existing type with a new version according to PowerCombiner

/// Adds the PowerType to the end of the PowerPack unless the PowerPack contains 
/// a PowerType with the same base unit type
/// if so, then combines the exponents of both power types via the PowerCombiner
/// @tparam PowerType power type to add to the pack
/// @tparam PowerPack existing pack of power types
/// @tparam PowerCombiner operation on exponents to combine powers of the same base unit
template<typename PowerType, typename PowerPack, template <typename, typename> class PowerCombiner>
struct ConsPowerType {
private:
    using Head = pack_head_t<PowerPack>;
    using Tail = pack_tail_t<PowerPack>;
public:
    /*
        if head = x then
            (combine head x) :: tail
        else
            head :: (append_power_type tail x)
    */
    using Type = type_cond_t<same_unit_v<PowerType, Head>,
        pack_cons_t<PowerCombiner<PowerType, Head>, Tail>,
        pack_cons_t<Head, typename ConsPowerType<PowerType, Tail, PowerCombiner>::Type>>;
};

// base case for `append_power_type [] x`
template<typename PowerType, template <typename, typename> class PowerCombiner>
struct ConsPowerType<PowerType, Pack<EmptyPack>, PowerCombiner> {
    using Type = Pack<PowerType>;
};

/// ---------------
/// SameBaseUnit 

template<typename A, typename B>
struct same_power_unit : std::false_type {};

// specialization for units of the same base type
template<typename Unit, int32_t numA, int32_t numB, int32_t denA, int32_t denB>
struct same_power_unit<PowerType<Unit, numA, denA>, PowerType<Unit, numB, denB>> 
    : std::true_type {};

/// True if two power types have the same base unit. Their powers may not match
template<typename UnitA, typename UnitB>
constexpr bool same_unit_v = same_power_unit<UnitA, UnitB>::value;

I choose the convention where all rational exponents are in simplest form and if the rational is negative, the numerator is represented as a negative number and not the denominator. Enforcing this allows unique representations for a single power, so we don’t need to worry about recognizing things like $\frac{1}{1} = \frac{2}{2} = \frac{-3}{-3}$.

With this, we can define an implementation of a type that satisfies the concept of a PowerCombiner, for adding exponents.

To do this, we use constexpr functions and type deduction by determining what the resulting power should be and returning a default-constructed value with that type. In C++, if an exception is thrown during compile-time evaluation of a constexpr function, that exception becomes a compiler error. Typically, it will say <function-name> failed to produce a constant value.

In C++, the most specific overload is prioritized, so the implementation of power_add that throws will only be chosen if a user tries to call power_add on two power types with different base unit types.

template<typename Unit, int32_t numA, int32_t denA, int32_t numB, int32_t denB>
constexpr auto power_add(PowerType<Unit, numA, denA>, PowerType<Unit, numB, denB>) {
    constexpr auto r = Rational(numA, denA) + Rational(numB, denB);
    return PowerType<Unit, r.num, r.den>{};
}

// This overload is needed for type checking, but should never get called as two 
// types of different base unit types should not be added together
template<typename A, typename B>
constexpr A power_add(A a, B) {
    throw std::invalid_argument("Cannot add two units of different base units");
    return a;
};

/// Combines the powers of two power types by adding them
/// Requires the power types share the same base unit
template<typename UnitA, typename UnitB>
using add_unit_powers_t = decltype(
    power_add(std::declval<UnitA>(), std::declval<UnitB>()));

/// Adds a power type to a power pack, combining it by adding exponents
/// if the unit type is already present in the pack
/// If the unit type is not present, adds the new power type to the END of the pack
template<typename PowerType, typename PowerPack>
using cons_power_add_t = typename ConsPowerType<PowerType, PowerPack, add_unit_powers_t>::Type;

Now, we can define a way for us to concatenate power packs while also merging any repeated units:

/// Returns a pack of power types of the types in PackA followed by the types in PackB
/// Combines the powers of the unit types that are the same in PackA and PackB
/// @tparam PackA first pack of power types
/// @tparam PackB second pack of power types
/// @tparam PowerCombiner operation on exponents to cons the powers of the same base unit
/// @{
template<typename PackA, typename PackB, template<typename, typename> class PowerCombiner>
struct PowerPackCombine;

// Base Case: [] @ _
template<typename PackB, template<typename, typename> class PowerCombiner>
struct PowerPackCombine<Pack<EmptyPack>, PackB, PowerCombiner> {
    using Type = PackB;
};

// Base Case: _ @ []
template<typename PackA, template<typename, typename> class PowerCombiner>
struct PowerPackCombine<PackA, Pack<EmptyPack>, PowerCombiner> {
    using Type = PackA;
};

// Base Case: [] @ []
// again, this is technically redundant, but we need it to avoid ambiguity
template<template<typename, typename> class PowerCombiner>
struct PowerPackCombine<Pack<EmptyPack>, Pack<EmptyPack>, PowerCombiner> {
    using Type = Pack<EmptyPack>;
};

/// Merges the power packs, combining exponents if the unit types are the same
/// @tparam PackA first pack of power types
/// @tparam PackB second pack of power types
/// @tparam PowerCombiner operation to cons a power type onto a power pack
template<typename PackA, typename PackB, template<typename, typename> class PowerCombiner>
struct PowerPackCombine {
    using Type = PowerCombiner<pack_head_t<PackA>, 
        typename PowerPackCombine<pack_tail_t<PackA>, PackB, PowerCombiner>::Type>;
};

/// Combines the power packs, merging any power types that have the same
/// base unit type in both packs by adding exponents
template<typename PackA, typename PackB>
using power_pack_add_t = 
    typename PowerPackCombine<PackA, PackB, cons_power_add_t>::Type;

There’s one hiccup. Suppose we add the exponents of PowerType<Meters, 2, 1> with PowerType<Meters, -2, 1>. We’d end up with PowerType<Meters, 0, 1>. To deal with zero-powered power types, we always clean the type list after combining. This clean operation removes all zero-power types. If we didn’t do this, then a power list like Pack<PowerType<Meters, 0, 1>, PowerType<Seconds, 1, 1>> would be a different type from Pack<PowerTyper<Seconds, 1, 1>> when they are the same.

template<typename PowerType>
struct IsZeroPowerType : std::false_type {};

template<typename Unit, int32_t den>
struct IsZeroPowerType<PowerType<Unit, 0, den>> : std::true_type {};

/// Determines if a power type is raided to the zeroth power zero
template<typename PowerType>
constexpr bool is_zero_power_v = IsZeroPowerType<PowerType>::value;

/// Removes any power types in the pack with a zero exponent
/// @{
template<typename PowerPack>
struct CleanPowerPack {
private:
    using Head = pack_head_t<PowerPack>;
    using Tail = pack_tail_t<PowerPack>;
public:
    /*
        if is_zero_power hd then
            clean_power_pack tl
        else
            hd :: (clean_power_pack tl)
    */
    using Type = type_cond_t<is_zero_power_v<Head>,
        typename CleanPowerPack<Tail>::Type,
        pack_cons_t<Head, typename CleanPowerPack<Tail>::Type>>;
};

template<>
struct CleanPowerPack<Pack<EmptyPack>> {
    using Type = Pack<EmptyPack>;
};
/// @}

/// Removes all zero-powered power types from the pack
template<typename PowerPack>
using clean_power_pack_t = typename CleanPowerPack<PowerPack>::Type;

Instead of creating a whole new operation for division, we can represent division as multiplication by the negated powers. So PowerType<Meters, 1, 1> / PowerType<Seconds, 1, 1> = PowerType<Meters, 1, 1> * PowerType<Seconds, -1, 1>. Then we just need to define a type to negate all of the powers in the type list:

template<typename PowerType>
struct NegatePowerType {};

// negates a single rational number
template<typename Unit, int32_t num, int32_t den>
struct NegatePowerType<PowerType<Unit, num, den>> {
    using Type = PowerType<Unit, -num, den>;
};

/// Negates the power that a unit type is raised to
template<typename PowerType>
using negate_power_t = typename NegatePowerType<PowerType>::Type;

// negation recursive case
// lst matches with hd :: tl
template<typename PowerPack>
struct NegatePowerPack {
    // -hd :: (negate_power_list tl)
    using Type = pack_cons_t<negate_power_t<pack_head_t<PowerPack>>, 
        typename NegatePowerPack<pack_tail_t<PowerPack>>::Type>;
};

// negation base case
// lst matches with []
template<>
struct NegatePowerPack<Pack<EmptyPack>> {
    using Type = Pack<EmptyPack>;
};

/// Negates all the powers of the power types in the power pack
template<typename PowerPack>
using negate_power_pack_t = typename NegatePowerPack<PowerPack>::Type;

And the final piece of the puzzle is being able to exponentiate unit types by rational powers. We won’t worry about being able to exponentiate a unit by another unit.

/// Exponentiates a power type by `powerNum / powerDen`
/// @{
template<typename PowerUnit, int32_t powerNum, int32_t powerDen>
struct RaiseUnitPower {};


template<typename Unit, int32_t unitNum, int32_t unitDen, 
    int32_t powerNum, int32_t powerDen>
struct RaiseUnitPower<PowerType<Unit, unitNum, unitDen>, powerNum, powerDen> {
private:
    // multiplies the exponents
    // taking care of putting the rational into its unique form
    constexpr static auto r = 
        Rational(unitNum, unitDen) * Rational(powerNum, powerDen);
public:
    using Type = PowerType<Unit, r.num, r.den>;
};
/// @}

// recursive case
template<typename PowerPack, int32_t powerNum, int32_t powerDen>
struct RaisePowerPack {
private:
    using Head = pack_head_t<PowerPack>;
    using Tail = pack_tail_t<PowerPack>;
public:
    /*
        (pow hd power_num power_den) :: raise_power tl power_num power_den
    */
    using Type = pack_cons_t<typename RaiseUnitPower<Head, powerNum, powerDen>::Type,
        typename RaisePowerPack<Tail, powerNum, powerDen>::Type>;

};

// base case
template<int32_t powerNum, int32_t powerDen>
struct RaisePowerPack<Pack<EmptyPack>, powerNum, powerDen> {
    using Type = Pack<EmptyPack>;
};


/// Exponentiates all power types in a power pack by `powerNum / powerDen`
/// @tparam PowerPack The power pack to raise
/// @tparam powerNum The numerator of the power to raise to
/// @tparam powerDen The denominator of the power to raise to
template<typename PowerPack, int32_t powerNum, int32_t powerDen>
using raise_power_pack_t = typename RaisePowerPack<PowerPack, powerNum, powerDen>::Type;

The complexity and verbosity come from needing to encode all of these things in C++ templates. I find it can be easier to understand what is going on by trying to write an implementation of the algorithm in your favorite functional language and matching up the functional and metaprogramming implementations.

I should further mention that by writing this post retrospectively, I am skipping a very crucial step: testing. In reality, I implemented this with a form of TDD. Every piece was immediately tested and I moved on to the next thing only after I trusted my implementation was mostly correct.

We are finally able to move on to the last step.

Unit Operations

We already discussed addition and subtraction. So let’s now look at multiplication and division:


template<typename T, scale_t scaleA, typename SemanticA, typename UnitPowerPackA,
scale_t scaleB, typename SemanticB, typename UnitPowerPackB>
constexpr auto operator*(const Unit<T, scaleA, SemanticA, UnitPowerPackA>& a,
    const Unit<T, scaleB, SemanticB, UnitPowerPackB>& b)
{
    // only requirement is that the units' base numeric type is the same


    using NewPowerPack = clean_power_pack_t<
        sort_unit_pack_t<power_pack_add_t<UnitPowerPackA, UnitPowerPackB>>>;
    // concatenate power lists, merging duplicates
    // we sort and clean the lists so that each power list is unique
    using NewSemantic = clean_power_pack_t<
        sort_unit_pack_t<power_pack_add_t<SemanticA, SemanticB>>>;
    if constexpr (std::is_same_v<NewPowerPack, Pack<EmptyPack>>) {
        // if all powers cancel out, return a scaler

        // `constexpr if` is determined at compile time
        return static_cast<T>(a.val * static_cast<double>(scaleA) 
            * b.val * static_cast<double>(scaleB));
    } else {
        return Unit<T, scaleA * scaleB, NewSemantic, NewPowerPack>
            (a.val * b.val);
    }
}

template<typename T, scale_t scale, typename Semantic, typename UnitPowerPack>
constexpr auto operator/(T b, const Unit<T, scale, Semantic, UnitPowerPack>& a)
{
    // inverse() member function of `Rational` inverts the number
    return Unit<T, scale.inverse(), negate_power_pack_t<Semantic>, 
        negate_power_pack_t<UnitPowerPack>>(b / a.val);
}

template<typename T, scale_t scaleA, typename Semantic, typename UnitPowerPackA,
    scale_t scaleB, typename SemanticB, typename UnitPowerPackB>
constexpr auto operator/(const Unit<T, scaleA, Semantic, UnitPowerPackA>& a,
    const Unit<T, scaleB, SemanticB, UnitPowerPackB>& b)
{
    return a * (T{1} / b);
}

As discussed earlier, the division is implemented in terms of multiplication. Now let’s look at exponentiation:

template<int32_t powerNum, int32_t powerDen, typename T, scale_t scale, typename Semantic, typename Units>
auto pow(const Unit<T, scale, Semantic, Units>& a)
{
    return Unit<T, Rational(1), raise_power_pack_t<Semantic, powerNum, powerDen>, 
        raise_power_pack_t<Units, powerNum, powerDen>>
            (std::pow(a.val * static_cast<double>(scale), 
                static_cast<double>(powerNum) / powerDen));
}

Our unit pow function is slightly more restrictive than a general pow function. In particular, it requires that the power is known at compile time. For most equations where this library could be applied, such exponents would be known ahead of time and the restriction is reasonable.

Note that because this relies on std::pow, pow is the only unit operation that cannot be performed at compile time.

Usage

Let’s recap what we did:

Created a Unit type with a rational scale factor, a list of rational powers of semantic types, and a list of rational powers of unit types
Ensured Units are represented uniquely
Provided arithmetic operations on Units
Provided implicit conversion between Units with the same base type but different scales

To see how this works, suppose we define volume and area like so:

// volume has units m^3 and is a length * width * height
using volume_t = Unit<double, Rational(1), 
    sort_unit_pack_t<Pack<PowerType<Length, 1, 1>, PowerType<Width, 1, 1>, PowerType<Depth, 1, 1>>>,
    Pack<PowerType<Meters, 3, 1>>>;


// area has units m^2 and is a length * width
using area_t = Unit<double, Rational(1), 
    sort_unit_pack_t<Pack<PowerType<Length, 1, 1>, PowerType<Width, 1, 1>>>,
    Pack<PowerType<Meters, 2, 1>>>;

// a distance is measured in m and has some semantic type
template<typename Semantic>
using dist_t = Unit<double, Rational(1), Pack<PowerType<Semantic, 1, 1>>, Pack<PowerType<Meters, 1, 1>>>;

Then we can create a function that takes a volume, divides it by a depth, and returns an area:

constexpr area_t foo(volume_t vol) {
    constexpr auto val = dist_t<Depth>(7);
    return vol / val + vol / val;
}

Note that, if we divided by a length, width, or any other unit, we’d get a compiler error!

constexpr area_t foo(volume_t vol) {
    constexpr auto val = dist_t<Width>(7);
    return vol / val + vol / val; // error: no viable conversion
}

Now we can call the function foo, but only by passing a volume:

constexpr auto x = dist_t<Length>(20);
constexpr auto y = dist_t<Width>(10);
constexpr auto z = dist_t<Depth>(5);

constexpr auto w = foo(x * y * z); // works!

constexpr auto w = foo(z * x * y); // also works!

If we make a typo:

constexpr auto w = foo(x * x * z); // error!: no viable conversion

We get a compiler error! This is because a Pack<PowerType<Length, 2, 1>, PowerType<Depth, 1, 1>> is not a volume. This not only works for semantic unit subcategories but also the regular unit types.

Now let’s define some more units like so:

using sec_t = Unit<double, Rational(1), Pack<EmptyPack>, Pack<PowerType<Seconds, 1, 1>>>;

using hour_t = Unit<double, Rational(3600), Pack<EmptyPack>, Pack<PowerType<Seconds, 1, 1>>>;

using km_t = Unit<double, Rational(1000), Pack<EmptyPack>, Pack<PowerType<Meters, 1, 1>>>;

using meter_t = Unit<double, Rational(1), Pack<EmptyPack>, Pack<PowerType<Meters, 1, 1>>>;

using meter_per_sec = decltype(std::declval<meter_t>() / std::declval<sec_t>());

constexpr auto get_speed(meter_per_sec speed) {
    return speed;
}

One thing to note is that for each dimension (distance, velocity, etc.), we require a single unit that all other units of that dimension are based on. This isn’t too unreasonable since the “cheeseburger per bald eagle air-speed” system (or as a Chemistry professor I had liked to call it: “stupid American units”) is defined in terms of the metric system.

Then we have:

constexpr auto f = meter_t{10} / sec_t{2};
constexpr auto t = km_t{10} / hour_t{2};

static_assert(get_speed(f).val == 5);

// km / hour implicitly converted to meters / sec
static_assert(get_speed(t).val == 5 / 3.6);

We see that Units that differ only in scale can be used interchangeably and the conversion functions handle everything for us! Therefore, this system solves the two main problems we aimed to tackle: mismatch errors and conversion errors. Mismatches in scale are handled automatically, and mismatches in unit or semantic type cause compiler errors. This ensures that conversions are correct (assuming the type definitions are correct) since an incorrect conversion will cause a mismatch error at compile time or be automatically handled if the mismatch is just a difference of scale.

There is a downside to this, however. Extensive usage of many different units can cause code bloat since each version of the Unit class must be instantiated. Furthermore, unless explicit template instantiation is used, this can increase compile time.

Despite this, I think this was an interesting experiment using the type system to eliminate a class of programmer errors. This also demonstrates the feasibility of achieving this in a common general-purpose programming language with nothing more than a header-only library.

Source