Compact C++ Maths APIs without copy-paste

Article

I've written dozens of maths APIs over the last 30 years and I've roundabout settled on designs that allow me to share C++ code with CUDA and HLSL. These are the properties that are important to me:

Straight structs for each type, no templates, as I want to forward declare them.
Vector types of varying precision.
All operations expressed as procedural, free functions.
Easy to maintain and extend, with as little duplication as possible.
Takes advantage of the target ABI to make the typical path as efficient as possible (e.g. __vectorcall).
Uses an entirely separate API for long streams of SIMD code.
Minimal reliance on compilers optimising code and very fast debug builds.
Difference between packed and unpacked types is explicit.
Compile times must be lightning fast, across all builds, and not drag in heavy STL dependencies.

Surprisingly (or not), C++ still makes this difficult. Were I designing a new language, these problems would be front and center, but why bother? Nobody would use it.

Anyway, the outcome is simple; code like this should work on C++, CUDA and HLSL:

float3 ProjectOntoAxis(float3 v, float3 axis)
{
    return axis * dot(v, axis);
}

void test()
{
    float3 a, b;

    float3 c = min(a, b);
    float  d = dot(a, b);
    float3 e = normalize(c);
    float3 f = a * 2.0f;
    float3 g = ProjectOntoAxis(f, e) * d;
}

CUDA has similarly-named vector types but you are required to implement the functions yourself, while HLSL has both the types and the functions.

This post will focus on everything above, excluding the SIMD vector streaming API, and attempt to remove all code duplication previous attempts have failed to do.

The final API

So we're in no doubt what this post achieves, here's a stripped down version of the API:

// Vector types are trivial, with no inheritance and allow for aggregate initialisation
struct float2  { float  x, y; };
struct float3  { float  x, y, z; };
struct float4  { float  x, y, z, w; };
struct double2 { double x, y; };
struct double3 { double x, y, z; };

// Public maths API with all operations explicitly declared
constexpr float2  __vectorcall operator*(float2,  float);
constexpr float3  __vectorcall operator*(float3,  float);
constexpr float4  __vectorcall operator*(float4,  float);
constexpr double2 __vectorcall operator*(double2, double);
constexpr double3 __vectorcall operator*(double3, double);
constexpr float  __vectorcall dot(float2,  float2);
constexpr float  __vectorcall dot(float3,  float3);
constexpr float  __vectorcall dot(float4,  float4);
constexpr double __vectorcall dot(double2, double2);
constexpr double __vectorcall dot(double3, double3);
constexpr float2  __vectorcall min(float2,  float2);
constexpr float3  __vectorcall min(float3,  float3);
constexpr float4  __vectorcall min(float4,  float4);
constexpr double2 __vectorcall min(double2, double2);
constexpr double3 __vectorcall min(double3, double3);

// Complete implementation of all functions above
template<typename V, typename S, auto... M>
struct VectorOps
{
    friend constexpr V __vectorcall operator*(V v, S s)
    {
        return {((v.*M) * s)...};
    }

    friend constexpr S __vectorcall dot(V a, V b)
    {
        return (((a.*M) * (b.*M)) + ...);
    }

    friend constexpr V __vectorcall min(V a, V b)
    {
        return {
            ((a.*M < b.*M) ? a.*M : b.*M)...
        };
    }
};

// Register all supported vector types
template<typename V>
struct VectorOpsMap;
template<> struct VectorOpsMap<float2> : VectorOps<float2, float, &float2::x, &float2::y> {};
template<> struct VectorOpsMap<float3> : VectorOps<float3, float, &float3::x, &float3::y, &float3::z> {};
template<> struct VectorOpsMap<float4> : VectorOps<float4, float, &float4::x, &float4::y, &float4::z, &float4::w> {};
template<> struct VectorOpsMap<double2> : VectorOps<double2, double, &double2::x, &double2::y> {};
template<> struct VectorOpsMap<double3> : VectorOps<double3, double, &double3::x, &double3::y, &double3::z> {};

We have five vector types, 15 functions and only three implementations. This is a spectacular win for an easy to maintain maths API:

If you want to add a new type, like double4 or int3, add your function declaration and specialise VectorOpsMap. VectorOps will automatically take care of the implementation for you.
If you want to add a new function, add it once and it will work for all supported vector types.
Pack expansion takes care of adapting to varying vector dimensions.

Why no inheritance is important

A typical way to achieve the above is to use inheritance on the vector types themselves:

struct VectorOps { };
struct float3 : VectorOps { float x, y, z; };

This is way more concise as you don't need to do any explicit registration, and it was the first approach I tried. However it has various trivial limitations, like messy aggregate initialisation or inability to extend a 3rd-party library. More importantly:

Inheritance breaks MSVC's ability to use __vectorcall register passing completely

__vectorcall is a Microsoft ABI calling convention that passes certain vector arguments to functions in SIMD registers, instead of passing them through memory. Compilers, especially in debug, can have a very hard time optimising around these boundaries and the memory traffic becomes surprisingly obese.

With a vector type inheriting from VectorOps, you lose it all, irrespective of how you annotate your functions. I did a whole bunch of tests to prove this, after first observing it in my own maths library:

https://godbolt.org/z/v9vzz6j79

Annoyingly, the ABI spec does not document this exception. And I could not tell you if this is a bug but, as you should already know: it really doesn't matter if it is a bug.

How these tricks combine

The code at first glance can look incredibly bizarre to somebody who's fed up with filtering through multiple C++ revisions to see if there's anything useful in there that won't inadvertently destroy your product¹. Specifically: how the hell can you call those class methods without a class instance, from the global namespace?

1. Friend functions that step out of the class

A feature of C++ I only became aware of recently is that friend functions defined in a class aren't actually methods; they're free functions that belong to the parent namespace:

struct MyStruct
{
    friend void Do()
    {
    }
};

// Make the function visible to ordinary name lookup
void Do();

void test()
{
    // This works!
    Do();
}

There is no class instance required here to call the function and the function has no implicit this pointer. What makes this particularly powerful is that the type can be templated and that it can make up for the lack of partial function specialisation support in C++:

enum Kind
{
    Value,
    Pointer,
};

template<typename T>
struct Data
{
    T value;

    friend constexpr Kind kind(Data)
    {
        return Value;
    }
};

// Specialise the implementation for all pointer types
template<typename T>
struct Data<T*>
{
    T* value;

    friend constexpr Kind kind(Data)
    {
        return Pointer;
    }
};

static_assert(kind(Data<int>{42}) == Value);
static_assert(kind(Data<float*>{nullptr}) == Pointer);
static_assert(kind(Data<double*>{nullptr}) == Pointer);

That's a feature for another day but it's nice to have in your box of tools.

2. Member pointers as variadic template arguments

Taking pointers to class members is an age-old technique that can be combined with the newer auto keyword to not care about the member type itself. While use of auto is generally discouraged where you want readers to better understand your code, use of auto in generic code actually simplifies it to a point of making it more understandable.

auto member = &float3::x;
float3 v{1, 2, 3};
float x_value = v.*member;

This becomes useful with variadic template arguments describing a pack:

template <typename V, typename S, auto... M>
struct VectorOps;

Registration supplies the addresses of the members of each type (e.g. &float3::x, &float3::y, &float3::z) and the following pack expansion:

return {((v.*M) * s)...};

expands into:

return {v.x * s, v.y * s, v.z * s};

It's this particular observation that allows you to avoid the use of std::tie or ordering the members and using std::index_sequence (or your own equivalent), or even providing an operator[] with a switch in it that does not optimise away in debug builds.

3. External type registration

As mentioned already, the typical combination of all these techniques ends with all the friend functions bundled into a class that each type inherits. We can get around this by moving the inheritance into a separate registration class:

template<typename V>
struct VectorOpsMap;

// Specialise for each type
template<> struct 
VectorOpsMap<float3> : VectorOps<float3, float, &float3::x, &float3::y, &float3::z>
{
};

The downside is that Argument Dependent Lookup will not find these functions because the vector types no longer inherit from VectorOps. But for me that's a price I'm willing to pay for __vectorcall, especially as all these core product types are in the global namespace anyway so that there is seamless interop between C++, CUDA and HLSL.

MSVC considerations

It's been this way for decades and no amount of pleading with people to not adopt the latest C++ features has ever discouraged the enthusiast: MSVC is always behind the curve with code generation quality and parser robustness.

Look at this lovely bit of code:

friend constexpr S __vectorcall dot(V a, V b)
{
    return (((a.*M) * (b.*M)) + ...);
}

Fully generic to any vector type or rank. But it doesn't work in MSVC which complains about the + fold expression, with error C3520. You instead have to do this:

friend constexpr S __vectorcall dot(V a, V b)
{
    S result = a.x * b.x + a.y * b.y;
    if constexpr (sizeof...(M) >= 3)
        result += a.z * b.z;
    if constexpr (sizeof...(M) >= 4)
        result += a.w * b.w;
    return result;
}

I mean it's great that we can use if constexpr as a fallback here, but we've gone from an implementation supporting any dimension to explicitly handling 2, 3 or 4. To be clear: there are no runtime branches here, this is all compile time.

You will have the same problem with operators like += that can be solved in a similar way.

More complex implementations

So far the API implementations have been pretty trivial so I'll show some more complicated examples.

As with the simpler template techniques you may have been using already, you can generate a single function for all types that reuse existing functions that use this new technique:

friend V __vectorcall slerp(V a, V b, S t)
{
    // Find the component of b perpendicular to a
    S cosine = min(max(dot(a, b), S(-1)), S(1));
    V relative = b - a * cosine;
    S sine = length(relative);

    // Preserve the endpoints
    if (t == S(0)) return a;
    if (t == S(1)) return b;
    if (sine < S(1e-5)) return a;

    // Rotate a towards b along their common great circle
    S theta = atan2(sine, cosine) * t;
    return a * cos(theta) + relative * (sin(theta) / sine);
}

You can also embed the fold expressions to reduce function calls in debug, as with this safer normalise function which uses a Scaled Sum of Squares approach to avoid overflow and gradual underflow:

friend V __vectorcall normalize(V v)
{
    // Find the largest absolute component
    S scale = 0;
    ((scale = max(scale, abs(v.*M))), ...);
    if (scale == 0) return {};

    // Scale first to prevent overflow or underflow when calculating the squared length
    V u{((v.*M) / scale)...};

    // Calculate the length without calling dot() in debug builds
    S sum = 0;
    ((sum += (u.*M) * (u.*M)), ...);

    // Normalise the scaled vector
    S inv_len = S(1) / sqrt(sum);
    return {((u.*M) * inv_len)...};
}

The shared implementation isn't limited to trivial component-wise operations.

Prior Art

Before trying to solve this, I asked all the LLMs to fix the problem for me. They couldn't. They enjoyed throwing around words like std::tie and std::index_sequence which I learned long ago to associate with people who don't like their debug builds (or even, shock horror, debugging itself). This unfortunately implies we are in the minority, as the training data must be overwhelmingly stacked in the non-debuggerers' favour.

So I asked the LLMs to dig up prior art instead. Highlighting progress in the community is important, none of us stand alone. And we found some that were pretty close!

Arthur O'Dwyer's An example of the Barton–Nackman trick is a useful explanation of how class templates can generate concrete, non-template friend functions.
Valentin Galea's Rapid Prototyping of Graphics Shaders in Modern C++ uses hidden friend functions, generated operations and fold expressions while paying particular attention to debug performance.
Järvi, Willcock and Lumsdaine's Concept-Controlled Polymorphism is a delightful early exploration of external type registration over inheritance.
The “private member robber” idiom, from Johannes Schaub’s Access to private members, uses the same member-pointer and friend injection mechanics to break into a class from outside.
CUDA's helper_math.h is a great example of a very clean looking API that contains the very pains we're trying to solve.
Then of course there is pycgen that I wrote years ago to automate the generation of repetitive code blocks like this. It was way more powerful than this, though, able to write arbitrary Python code to generate, for example, Spherical Harmonic expressions directly rather than relying on magic constants.
But it all started in 2006 for me with an old research paper titled Generative Programming and Matrix Multiplication, exploring similar problems in C++ that had to use terribly complicated template constructions.
I also created a proof-of-concept procedural math API back then called pM.

Conclusion

Relief is a word that best describes my discovery of this technique.

I spent 10+ years of my 25 year C++ dependence as a self-professed "language lawyer". It was a waste of my life and it's a waste of yours. Now get off my lawn. ↩︎