C++ Vector3d class

Question

As part of a physics simulation written in C++ using SFML I needed a vector3d class.

Here is my implementation:

#include <cmath>
#include <ostream>

template <typename T>
class vector3d {

public:
  T x{};
  T y{};
  T z{};

  T norm() { return std::hypot(x, y, z); }

  vector3d<T> unit() { return vector3d<T>(*this) / norm(); }

  T dot(const vector3d<T>& RHS) { return x * RHS.x + y * RHS.y + z * RHS.z; }

  vector3d<T> cross(const vector3d<T>& RHS) {
    return vector3d<T>{y * RHS.z - z * RHS.y, z * RHS.x - x * RHS.z, x * RHS.y - y * RHS.x};
  }

  bool operator==(const vector3d<T>& rhs) {
    return std::abs(x - rhs.x) < std::numeric_limits<T>::epsilon() &&
           std::abs(y - rhs.y) < std::numeric_limits<T>::epsilon() &&
           std::abs(z - rhs.z) < std::numeric_limits<T>::epsilon();
  }

  bool operator!=(const vector3d<T>& rhs) { return !(*this == rhs); }

  // clang-format off
  vector3d<T>& operator+=(vector3d<T> RHS) { x += RHS.x; y += RHS.y; z += RHS.z; return *this; }
  vector3d<T>& operator-=(vector3d<T> RHS) { x -= RHS.x; y -= RHS.y; z -= RHS.z; return *this; }
  vector3d<T>& operator*=(vector3d<T> RHS) { x *= RHS.x; y *= RHS.y; z *= RHS.z; return *this; }
  vector3d<T>& operator/=(vector3d<T> RHS) { x /= RHS.x; y /= RHS.y; z /= RHS.z; return *this; }

  vector3d<T>& operator*=(T s) { x *= s; y *= s; z *= s; return *this; }
  vector3d<T>& operator/=(T s) { x /= s; y /= s; z /= s; return *this; }

  friend vector3d<T> operator+(const vector3d<T>& LHS, const vector3d<T>& RHS) { return vector3d<T>(LHS) += RHS; }
  friend vector3d<T> operator-(const vector3d<T>& LHS, const vector3d<T>& RHS) { return vector3d<T>(LHS) -= RHS; }
  friend vector3d<T> operator*(const vector3d<T>& LHS, const vector3d<T>& RHS) { return vector3d<T>(LHS) *= RHS; }
  friend vector3d<T> operator/(const vector3d<T>& LHS, const vector3d<T>& RHS) { return vector3d<T>(LHS) /= RHS; }

  friend vector3d<T> operator*(const vector3d<T>& v, T s) { return vector3d<T>(v) *= s; }
  friend vector3d<T> operator*(T s, const vector3d<T>& v) { return vector3d<T>(v) *= s; }
  friend vector3d<T> operator/(const vector3d<T>& v, T s) { return vector3d<T>(v) /= s; }
  // note that scalar / vector3d does not make any logical sense
  // clang-format on

  friend std::ostream& operator<<(std::ostream& os, vector3d<T> v) {
    return os << "[" << v.x << ", " << v.y << ", " << v.z << "]";
  }
};

using vec3  = vector3d<double>;
using vec3f = vector3d<float>;
using vec3l = vector3d<long double>;

And its potential usage:

#include "vector3d.h"
#include <iostream>

int main() {
  vec3 z;
  std::cout << "z = " << z << "\n";

  auto a = vec3{1, 2, 3};
  auto b = vec3{4, 5, 6};
  std::cout << "a = " << a << "\n";
  std::cout << "b = " << b << "\n";
  std::cout << "a + b = " << a + b << "\n";

  auto c = 2 * b;
  std::cout << "c = " << c << "\n";
  
  auto d = b * 2;
  std::cout << "d = " << d << "\n";
  
  auto e = b / 2;
  std::cout << "e = " << e << "\n";

  std::cout << "a.norm() = " << a.norm() << "\n";
  std::cout << "a.unit() = " << a.unit() << "\n";
  std::cout << "a.unit().norm() = " << a.unit().norm() << "\n";

  std::cout << "[4, 8, 10] . [9, 2, 7] = " << vec3{4, 8, 10}.dot(vec3{9, 2, 7}) << "\n";

  std::cout << "[1, 0, 0] x [0, 1, 0] = " << vec3{1, 0, 0}.cross(vec3{0, 1, 0}) << "\n";
  std::cout << "[2, 3, 4] x [5, 6, 7] = " << vec3{2, 3, 4}.cross(vec3{5, 6, 7}) << "\n";

}

• I tried to "do as the ints do" and make all basic operators available.
• I realise that +-*/ operators could be free functions but I preferred to write them as friends to avoid having to template every one of them. A valid style choice?
• I realise that the dot() and cross() could be free functions, but preferred the a.cross(b) usage syntax.
• Minor point: I chose to format the code to wider than normal as I believe it makes it more readable / makes the patterns more obvious.
• I added some type aliases for convenience.
• Is relying on aggregate initialisation (without any constructor) and the default copy constructor a good choice?

Any comments on the code and use of member / friend / free functions please?

Answer 1

This looks like a pretty good start. It compiles cleanly with a pedantic compiler, and the test program seems clean under Valgrind.

The first few members should not be modifying the object, so declare them const:

T norm() const { return std::hypot(x, y, z); }

vector3d unit() const { return vector3d{*this} / norm(); }

T dot(const vector3d& rhs) const { return x * rhs.x + y * rhs.y + z * rhs.z; }

vector3d cross(const vector3d& rhs) const {
    return {y * rhs.z - z * rhs.y, z * rhs.x - x * rhs.z, x * rhs.y - y * rhs.x};
}

Note also that unqualified vector3d within the scope of vector3d<T> means vector3d<T>, so we get to remove some clutter here.

Do the * and / operations of two vectors have any useful mathematical meaning? I'd just omit these:

  vector3d& operator*=(vector3d RHS) { x *= RHS.x; y *= RHS.y; z *= RHS.z; return *this; }
  vector3d& operator/=(vector3d RHS) { x /= RHS.x; y /= RHS.y; z /= RHS.z; return *this; }

(Later) Having been informed that these are the rarely-needed Hadamard product and its inverse, it may be appropriate to retain the functionality, but as an ordinary (named) member function rather than as an operator. This will reduce the chance of accidental misuse.

Also, don't use all-caps names for locals or arguments - it's idiomatic to reserve those names for preprocessor macros.

The friend functions are probably okay for now, but lean towards non-friend non-member functions, as this won't require any change when you add implicit conversions (e.g. from vector3d<float> to vector3d<double>) in future.

We can avoid the explicit copy if we take one argument by value, which will cause an implicit copy:

template<typename T>
auto operator+(vector3d<T> lhs, const vector3d<T>& rhs) { return lhs += rhs; }
template<typename T>
auto operator-(vector3d<T> lhs, const vector3d<T>& rhs) { return lhs -= rhs; }

template<typename T>
auto operator*(vector3d<T> v, T s) { return v *= s; }
template<typename T>
auto operator*(T s, vector3d<T> v) { return v *= s; }
template<typename T>
auto operator/(vector3d<T> v, T s) { return v /= s; }

The streaming operator << should not be modifying its vector argument, so pass that as a const-ref. We can micro-optimise by passing single characters rather than one-character literal strings.

Consider adding == and !=. Since C++20, we can kill both birds with a single stone:

bool operator==(const vector3d<T>&) const = default;

Consider constraining the template so that only arithmetic types can be used as the type T. I don't think a vector of std::complex makes much sense, for example.

---

Modified code

#include <cmath>
#include <concepts>
#include <ostream>

template<typename T>
    requires std::floating_point<T> or std::integral<T>
class vector3d
{
public:
    T x{};
    T y{};
    T z{};

    T norm() const { return std::hypot(x, y, z); }

    vector3d unit() const { return vector3d{*this} / norm(); }

    T dot(const vector3d& rhs) const { return x * rhs.x + y * rhs.y + z * rhs.z; }

    vector3d cross(const vector3d& rhs) const {
        return {y * rhs.z - z * rhs.y, z * rhs.x - x * rhs.z, x * rhs.y - y * rhs.x};
    }

    bool operator==(const vector3d&) const = default;

    vector3d& operator+=(const vector3d& rhs) { x += rhs.x; y += rhs.y; z += rhs.z; return *this; }
    vector3d& operator-=(const vector3d& rhs) { x -= rhs.x; y -= rhs.y; z -= rhs.z; return *this; }

    vector3d& operator*=(T s) { x *= s; y *= s; z *= s; return *this; }
    vector3d& operator/=(T s) { x /= s; y /= s; z /= s; return *this; }
};

template<typename T> auto operator+(const vector3d<T>& v) { return v; }
template<typename T> auto operator-(const vector3d<T>& v) { return vector3d<T>{} - v; }

template<typename T> auto operator+(vector3d<T> lhs, const vector3d<T>& rhs) { return lhs += rhs; }
template<typename T> auto operator-(vector3d<T> lhs, const vector3d<T>& rhs) { return lhs -= rhs; }

template<typename T, typename U> auto operator*(vector3d<T> v, U s) { return v *= s; }
template<typename T, typename U> auto operator*(U s, vector3d<T> v) { return v *= s; }
template<typename T, typename U> auto operator/(vector3d<T> v, U s) { return v /= s; }

template<typename T> auto& operator<<(std::ostream& os, const vector3d<T>& v) {
    return os << '[' << v.x << ", " << v.y << ", " << v.z << ']';
}

Answer 2

If you have an up-to-date compiler, you don't need to write operator!=.

std::abs(x - rhs.x) < std::numeric_limits<T>::epsilon()
The type T is a template argument. Does it have an epsilon? You should use if constexpr (std::numeric_limits<T>::is_exact) to conditionally use this logic.
But, does what you wrote even work? In fact, there is an example on the cppreference page that explains:

the machine epsilon has to be scaled to the magnitude of the values used
and multiplied by the desired precision in ULPs (units in the last place)

and you are missing that logic.

---

friend vector3d<T> operator+(const vector3d<T>& LHS, const vector3d<T>& RHS) { return vector3d<T>(LHS) += RHS; }
You are passing LHS by reference and then forcing a copy to be made with more verbage. It would be simpler to just pass the left argument by value:
friend vector3d<T> operator+(vector3d<T> LHS, const vector3d<T>& RHS) { return LHS += RHS; }

But, you can further optimize it for the case where the left argument is a temporary rvalue — something that is common when you write expressions. Consider z=a+b+c+d;. If the result of a+b can be modified in-place and not re-copied, you save a lot. This happens again with (temp)+d.

The easiest way to do that is to add a second form that takes an rvalue reference: friend vector3d<T> operator+(vector3d<T>&& LHS, const vector3d<T>& RHS) { return LHS += RHS; }
Notice that the body is the same. But I think it's more optimal to write {LHS+=RHS; return LHS;} The return is implicitly taken as return std::move{LHS};} because it is a simple name of a local variable or parameter.

---

friend std::ostream& operator<<(std::ostream& os, vector3d<T> v) {
What's the need to pass v by value in this function? It seems to be an unnecessary copying.

---

std::cout << "a.norm() = " << a.norm() << "\n";
Writing '\n' instead is more efficient. A string of one character is a lot more work to process.

---

There are a lot of places where const is missing.

---

I realise that +-*/ operators could be free functions but I preferred to write them as friends to avoid having to template every one of them. A valid style choice?

You did better than you realize! This is called the "hidden friend" idiom, and it means that these operators are only found during ADL. This reduces the size of the overload set for the first pass of every + (etc.) in the program, and reduces conflicts introduced by any implicit conversions.

Is relying on aggregate initialisation (without any constructor) and the default copy constructor a good choice?

Yes!

Answer 3

Floating point equality

  bool operator==(const vector3d<T>& rhs) {
    return std::abs(x - rhs.x) < std::numeric_limits<T>::epsilon() &&
           std::abs(y - rhs.y) < std::numeric_limits<T>::epsilon() &&
           std::abs(z - rhs.z) < std::numeric_limits<T>::epsilon();
  }

We've probably all heard about the dangers of floating point equality, and seen the suggestion to compare them something like that. But this is probably not really what you want, for some different reasons.

• Whether two vectors are "close", is a different question than whether each of the coordinates are "close". For example { 1e99, 1, 0 } and { 1e99, 0, 0 } would be considered different this way. That may be appropriate, but it may also not be. Those vectors have approximately the same length and have an angle of approximately zero degrees between them. On the other hand they are also legitimately different vectors by another measure: the difference between them is a nice vector with length 1.
• This measure of closeness considers 1.0 and nextafter(1.0, 10) to be different. The distance between them is (by definition) epsilon, and epsilon is not less than epsilon, so the test fails. Should that be how it works? Maybe! That depends on what you want to happen. It's a suspicious behaviour though: it makes this test a test for exact equality for numbers that are 1.0 or higher, which looks unintentional.
• This measure of closeness considered any pair of sufficiently small numbers to be equal, even they are relatively far apart. Again, maybe that's what you want, but it's suspicious.
• The expected amount of numerical error for a given value depends not only on that value, but also on how it was computed.

Since you are using these vectors in a specific context, there may be an appropriate implementation for operator== (in which case, I don't expect that the current one does what you need it to do), but there may also not be. It's entirely possible that some of your vectors need to be compared for closeness one way, and some other vectors in another way. Just to make up some random examples: 1 meter above the ground is still 1 meter above the ground when the horizontal distance from the origin is large, that doesn't suddenly become "pretty much on the ground" just because we're located in a different corner of the map. On the other hand, two directions with an angle of almost 0 between them are close, regardless of any differences between the coordinates.

It seems reasonable to me to not implement operator== at all, and always explicitly state how vectors are to be compared (length of their difference, angle, or something else).