How to efficiently alias a float as both a named member and an element of an array?

Question

I am implementing a particle-based fluid simulation. To represent vectors such as velocity, acceleration etc I have defined a class that looks like this

class Vec3f {
public:
    float x, y, z;

    // ... bunch of constructors, operators and utility functions
}

I'm using the library nanoflann for kd-tree searches. To accommodate for arbitrary class designs, nanoflann requires a user-defined adaptor class that the kd-tree class then queries to get info about the particle dataset. One function that the adaptor has to offer, as described in the nanoflann documentation is the following.

// Must return the dim'th component of the idx'th point in the class:
inline T kdtree_get_pt(const size_t idx, int dim) const { ... }

The problem is that this interface does not work seamlessly with the x, y, z representation. Naively, it would need to do something like this

inline float kdtree_get_pt(const size_t idx, int dim) const {
    switch(dim) {
    case 0:
        return particlearray[idx].x;
    case 1:
        return particlearray[idx].y;
    case 2:
        return particlearray[idx].z;
    }
}

Building and querying the kd-tree consumes a significant portion of my app's runtime and kd_tree_get_pt gets queried multiple times in the process so I need it to be optimized. The following solution should be faster.

class Vec3f {
public:
    float values[3];
    // ...
}

// Then inside the adaptor class
inline float kdtree_get_pt(const size_t idx, int dim) const {
    return particlearrray[idx].values[dim];
}

However, I much prefer the x, y, z interface for my equations. Now that the problem is clear, my question is how can I keep the x, y, z notation for my equations without making kdtree_get_pt suboptimal.

Solutions I have considered:

• Vec3f has member float values[3] and getters in the form of float& x(). The function call should be optimized away completely so this almost works but I do not want to have to add the parentheses in my equations if I can avoid it. eg I want to be able to write vec1.x - vec2.x instead of vec1.x() - vec2.x(). As far as I know, C++ does not offer a way to "disguise" a function call as a member variable, excluding preprocessor macros which I do not consider a safe solution.

• Vec3f has members float values[3] and float& x, y, z where the latter are initialized to point to the corresponding floats in the array. I thought they would be optimized away as they are known at compile time and obviously cannot change value after initialization, but even with optimizations on, MSVC++ seems to actually store the float&s as can be seen by sizeof(Vec3f) growing by 12 bytes after their addition. This doubles the storage size of my dataset which raises a concern for cache misses when working with arbitrarily large datasets.

• kdtree_get_pt uses float& values[3] to point to x, y, z. This might eliminate the branching cost, but I don't believe the extra level of indirection, nor the need to initialize all 3 references can be optimized away so it is presumably slower than the return particlearrray[idx][dim]` version.

• kdtree_get_pt uses reinterpret_cast or pointer magic to directly point into Vec3f's members. Given a Vec3f object's address, I believe x, y, z are guaranteed to be stored in that order with the first one stored at the same address as the one given by the & operator on the Vec3f object, but even so I'm confused as to whether there exists a well-defined way to observe them.

Answer 1

From a software engineering standpoint, it's best to expose the data through accessor and modifier functions only.

I would suggest:

class Vec3f
{
   public:

      float& operator[](size_t index) { return values[index]; }
      float operator[](size_t index) const { return values[index]; }

      float& x() { return values[0]; }
      float x() const { return values[0]; }

      float& y() { return values[1]; }
      float y() const { return values[1]; }

      float& z() { return values[2]; }
      float z() const { return values[2]; }

   private:
      float values[3];
}

Re: kdtree_get_pt uses reinterpret_cast or pointer magic to directly point into Vec3f's members.

That's a bad idea in general. However, I don't see that being a problem with my suggestion.

Answer 2

You should always check if the switch statements will really introduce branchings in the final compiled output. A tool that might help you there is godbolt.

For both of those code snippets (random and cout have been added to prevent complete removeal of the code):

#include<cstddef>
#include<array>
#include<iostream>
#include <ctime>

class Vec3f {
public:
    float values[3];
};

struct Test {
    std::array<Vec3f,100> particlearray;

    float kdtree_get_pt(const size_t idx, int dim) const {
        return particlearray[idx].values[dim];
    }
};

int main() {
   Test t;

   std::srand(std::time(0));
   int random_variable = std::rand();

   std::cout << t.kdtree_get_pt(random_variable,0);
   std::cout << t.kdtree_get_pt(random_variable,1);
   std::cout << t.kdtree_get_pt(random_variable,2) << std::endl;

   return 0;
}

and

#include<iostream>
#include<array>
#include<ctime>
#include<cstdlib>
#include<cstddef>


class Vec3f {
public:
    float x, y, z;
};

struct Test {
    std::array<Vec3f,100> particlearray;

    float kdtree_get_pt(const size_t idx, int dim) const {
        switch(dim) {
        case 0:
            return particlearray[idx].x;
        case 1:
            return particlearray[idx].y;
        case 2:
            return particlearray[idx].z;
        }
    }
};

int main() {
   Test t;

   std::srand(std::time(0));
   int random_variable = std::rand();

   std::cout << t.kdtree_get_pt(random_variable,0);
   std::cout << t.kdtree_get_pt(random_variable,1);
   std::cout << t.kdtree_get_pt(random_variable,2) << std::endl;


   return 0;
}

The access to x, y and z or values[dim] will be compiled (by gcc 7) to:

cvtss2sd xmm0, DWORD PTR [rsp+rbx]
cvtss2sd xmm0, DWORD PTR [rsp+4+rbx]
cvtss2sd xmm0, DWORD PTR [rsp+8+rbx]

Without any branching.

Answer 3

However, I much prefer the x, y, z interface for my equations. Now that the problem is clear, my question is how can I keep the x, y, z

Declare x,y,z as local references before calculation:

auto& [x1, y1, z1] = v1.values;
auto& [x2, y2, z2] = v2.values;
return x1*x2 + y1*y2 + z1*z2;

For pre-C++17, you need more verbose:

auto& x = values[0];
auto& y = values[1];
auto& z = values[2];

The compiler will not need to use any storage for these references.

This of course introduces some repetition; One line (in C++17) per vector per function.

Extra parentheses introduced by your first suggestion is another good way to go. Whether the introduction of parentheses is better or worse than local reference declaration boiler plate depends on the use case and personal preference.

---

Edit: Another alternative: Define operator[] and use named constants for indices.

namespace axes {
    enum axes {
        x, y, z
    };
}

struct Vec3f {
    float values[3];
    float& operator[](size_t index)       { return values[index]; }
    float  operator[](size_t index) const { return values[index]; }
};

// usage
using namespace axes;
return v1[x]*v2[x] + v1[y]*v2[y] + v1[z]*v2[z];

Answer 4

There is known technique for mixing up access via x, y, z and array indices using union of identical data types. Resolves problem with UB, sizeof() is 12 bytes, access time is as fast as it can be, one could use SIMD vector in very similar fashion. Code below tested with VS2017

#include <iostream>
#include <type_traits>

template <int Size> struct VectorBase {
    float _data[Size];

    float operator[](int Index) {
        return _data[Index];
    }
};

template <typename VectorType, int Index> struct ScalarAccessor {
    VectorType _v;

    operator float() const {
        return _v._data[Index];
    }

    float operator = (float x) {
        _v._data[Index] = x;
        return *this;
    }
};

union uuu {
        VectorBase<3>                    xyz;
        ScalarAccessor<VectorBase<3>, 0> x;
        ScalarAccessor<VectorBase<3>, 1> y;
        ScalarAccessor<VectorBase<3>, 2> z;
};

template <int Size> struct Vector {
    union
    {
        VectorBase<3>                    xyz;
        ScalarAccessor<VectorBase<3>, 0> x;
        ScalarAccessor<VectorBase<3>, 1> y;
        ScalarAccessor<VectorBase<3>, 2> z;
    };

    float operator[](int Index) {
        return xyz[Index];
    }
};

using Vec3f = Vector<3>;

int main() {
    Vec3f a;

    a.x = 1.0f;
    a.y = a.x + 3.0f;
    a.z = a.x * 3.0f;

    std::cout << sizeof(a) << "\n";
    std::cout << a.x << " " << a.y << " " << a.z << "\n";
    std::cout << a[0] << " " << a[1] << " " << a[2] << "\n";

    std::cout << std::is_standard_layout<VectorBase<3>>::value << "\n";
    std::cout << std::is_standard_layout<ScalarAccessor<VectorBase<3>, 0>>::value << "\n";
    std::cout << std::is_standard_layout<ScalarAccessor<VectorBase<3>, 1>>::value << "\n";
    std::cout << std::is_standard_layout<ScalarAccessor<VectorBase<3>, 2>>::value << "\n";
    std::cout << std::is_standard_layout<Vec3f>::value << "\n";
    std::cout << std::is_standard_layout<uuu>::value << "\n";

    return 0;
}

UPDATE

Here some C++ standard reading

I'm relying on the definition of standard-layout type 12.7 Classes

A class S is a standard-layout class if it: (7.1) — has no non-static data members of type non-standard-layout class (or array of such types) or reference, (7.2) — has no virtual functions (13.3) and no virtual base classes (13.1), (7.3) — has the same access control (Clause 14) for all non-static data members, (7.4) — has no non-standard-layout base classes, (7.5) — has at most one base class subobject of any given type ...

It is easy to check if all proposed classes are standard-layout - I've changed the code to check for that.

They are all layout-compatible, I believe

Also If a standard-layout class object has any non-static data members, its address is the same as the address of its first non-static data member.

Union is a standard-layout class as well, so we have classes aligned in union with only data member being array of the same type and size, and looks like standard requires it to be byte-by-byte compatible