Class to store and manipulate Cartesian coordinates and vectors

Question

I have made a C++ class to manipulate 3D vectors in Cartesian coordinates.

However, the performance of my class is much slower (about 2.5x) than simply using something like double p [3] and then running for loops for additions, subtractions, etc.

My specific concerns are:

• Is there something fundamentally wrong that I am doing, that is causing this slowdown?
• Is there a way to achieve what I am trying, but without as much of a slowdown?
• Is this just a fundamental limitation of classes and operator overloading that I have to live with, in return for the convenience?

Here is the class:

#ifndef COORD_H
#define COORD_H
#include <iostream>
#include <fstream>
#include <iomanip>
#include <cstring>
#include <complex>
#include <vector>
#include <map>
#include <stdlib.h>
#include "math.h"
// Useful reference: http://courses.cms.caltech.edu/cs11/material/cpp/donnie/cpp-ops.html
/*! \brief Class to store and manipulate points or position vectors in 3D, in Cartesian coordinates.*/
class vect3D
{
public:
 double x = 0.0, y = 0.0, z = 0.0;
 double tol = 1.0e-15;
 // Constructors
 vect3D() {};
 vect3D(std::initializer_list<double> RHS) { x = *RHS.begin(); y = *(RHS.begin()+1); z = *(RHS.begin()+2); };
 vect3D(std::vector<double> &RHS) { x = RHS[0]; y = RHS[1]; z = RHS[2]; };
 vect3D(double *RHS) { x = RHS[0]; y = RHS[1]; z = RHS[2]; };
 // Assignment
 vect3D& operator=(const vect3D &RHS) { x = RHS.x; y = RHS.y; z = RHS.z; return *this; };
 // Addition and subtraction
 vect3D& operator+=(const vect3D &RHS) { x += RHS.x; y += RHS.y; z += RHS.z; return *this; };
 vect3D& operator-=(const vect3D &RHS) { x -= RHS.x; y -= RHS.y; z -= RHS.z; return *this; };
 vect3D& operator+=(const double &RHS) { x += RHS; y += RHS; z += RHS; return *this; };
 vect3D& operator-=(const double &RHS) { x -= RHS; y -= RHS; z -= RHS; return *this; };
 vect3D operator+(const vect3D &RHS) { return vect3D(*this) += RHS; };
 vect3D operator-(const vect3D &RHS) { return vect3D(*this) -= RHS; };
 vect3D operator+(const double &RHS) { return vect3D(*this) += RHS; };
 vect3D operator-(const double &RHS) { return vect3D(*this) -= RHS; };
 // Scalar product and division
 vect3D& operator*=(const double &RHS) { x *= RHS; y *= RHS; z *= RHS; return *this; };
 vect3D& operator/=(const double &RHS) { x /= RHS; y /= RHS; z /= RHS; return *this; };
 vect3D operator*(const double &RHS) { return vect3D(*this) *= RHS; };
 vect3D operator/(const double &RHS) { return vect3D(*this) /= RHS; };
 friend vect3D operator*(double c, vect3D &vec) { return vec*c; };
 friend vect3D operator/(double c, vect3D &vec) { return vec/c; };
 // Comparisons
 bool operator==(const vect3D &RHS) { return ((x - RHS.x < x*tol) && (y - RHS.y < y*tol) && (z - RHS.z < z*tol)); };
 bool operator!=(const vect3D &RHS) { return !(*this == RHS); };
 bool operator>=(const vect3D &RHS) { return ((x >= RHS.x) && (y >= RHS.y) && (z >= RHS.z)); };
 bool operator<=(const vect3D &RHS) { return ((x <= RHS.x) && (y <= RHS.y) && (z <= RHS.z)); };
 bool operator>(const vect3D &RHS) { return !(*this <= RHS); };
 bool operator<(const vect3D &RHS) { return !(*this >= RHS); };
 // Euclidean norm
 double norm2() { return std::sqrt(std::pow(x, 2) + std::pow(y, 2) + std::pow(z, 2)); };
 friend double norm2(vect3D const &a) { return std::sqrt(std::pow(a.x, 2) + std::pow(a.y, 2) + std::pow(a.z, 2)); };
 // Dot product
 friend double dot(vect3D const &a, vect3D const &b) { return a.x*b.x + a.y*b.y + a.z*b.z; };
 // Cross product
 friend vect3D cross(vect3D const &a, vect3D const &b) { return {(a.y*b.z - a.z*b.y), (a.z*b.x - a.x*b.z), (a.x*b.y - a.y*b.x)}; };
 // Print to stream
 friend std::ostream& operator<<(std::ostream &stream, vect3D const &p) { return stream << "(" << p.x << ", " << p.y << ", " << p.z << ")" << std::flush; };
 // Function to explicitly return coordinates as an array of doubles, if needed
 void DoubleArray(double *v) { v[0] = x; v[1] = y; v[2] = z; return; };
 // To access coordinates using square bracket notation, for convenience
 std::vector<double *> p = std::vector<double *> {&x, &y, &z};
 double operator [] (int ii) const { return *(p[ii]); };
};
#endif

Thanks!

• Graipher: Hey, welcome to Code Review! If your code works but you want it to be better (in terms of speed or maintainability), you have come to the right place. Stack Overflow is a better place if you are looking for why something does not work or how to do something in the first place.
• Incomputable: Have you thought about inverting the access? E.g. store std::array<double, 3> and provide x(), y(), z() functions returning a reference to their corresponding elements.
• josh_eime: @Incomputable Now that you mention it, I refer to the x, y, z members almost as frequently through [ ] as through .x, .y etc., so I'll try out what you suggest and assess the impact on performance.

Answer 1

This is a general review, not a specific performance review.

---

Many of the operators ought to be const-qualified. Relational operators are of questionable value, given there isn't a natural ordering for coordinate vectors.

---

The constructors are fragile because they don't enforce the requirement that the caller provide at least three values. The initializer_list constructor is pointless really - just provide an ordinary 3-argument constructor that will be properly checked at the call site, and prefer initialization to assignment (this last may be a performance issue):

vect3D(double x, double y, double z)
 : x{x}, y{y}, z{z}
{}

---

I don't think there's any need for the assignment operator, unless you really don't want to copy tol as the compiler-generated one does - in which case, I'd recommend an explanatory comment, because it just looks like you forgot that member.

---

You seem to missing an include of <cmath> for std::pow() and std::sqrt(). However, you don't want to be using those functions, given that <cmath> provides a (much better behaved and possibly more efficient) std::hypot() function:

// Euclidean norm
double norm2() const { return std::hypot(x, y, z); };
friend double norm2(vect3D const &a) { return a.norm2(); };

(Note: in C++14 and earlier, std::hypot() only took two arguments, so you had to write std::hypot(std::hypot(x, y), z) or similar.)

---

The extra std::vector of pointers that each vect3D carries around with it is likely to be a performance drag. If you really want that help for indexing, consider a static constexpr array of member-pointers instead (which won't take space in the object, nor need to allocate when it's copied):

 // To access coordinates using square bracket notation, for convenience
 double operator[](int i) const { return this->*p[i]; };
 double& operator[](int i) { return this->*p[i]; };
private:
 static constexpr double vect3D::*p[] = { &vect3D::x, &vect3D::y, &vect3D::z };

---

Don't stream std::flush when writing output. If callers want output to be flushed, they will do this themselves; if not, they certainly don't want it to be imposed upon them.

---

Header rationalisation - in conjunction with the other improvements I recommend, we only need

#include <cmath>
#include <ostream>

• josh_eime: These are extremely helpful; thank you. In case it is helpful to others, I want to note that implementing your replacement of std::vector with the static constexpr improved performance by ~2x, so now the class is almost on par with using raw pointers.
• Toby Speight: I'm not sure whether you ever give different coordinates different tol values; if not, then making that member a static const would reduce the object size by 25%, probably with a corresponding speed boost.
• Incomputable: What do you think about this: a free standing template function component, which has component’s name as first argument and vector as second, input to function argument? That would allow crazy extensibility. Here is an example component<vec3d::x>(line3d). One could also write delegate like x_component(line3d) to make it less verbose.

Answer 2

Is this just a fundamental limitation of classes and operator overloading that I have to live with, in return for the convenience?

Not at all! C++ is all about abstractions with zero overhead. The design behind C++ is that it shouldn't leave any place for another language closer to the machine, while providing powerful abstractions to leverage native speed.

Further down a semi-complete review of your code (@TobySpeight already covered much ground).

#ifndef COORD_H
#define COORD_H
#include <iostream>
#include <fstream>
#include <iomanip>

That seems a bit heavy on I/O headers

#include <cstring>
#include <complex>
#include <vector>
#include <map>
#include <stdlib.h>

There's a bit of clean up to do here

//#include "math.h"
// Useful reference: http://courses.cms.caltech.edu/cs11/material/cpp/donnie/cpp-ops.html
/*! \brief Class to store and manipulate points or position vectors in 3D, in Cartesian coordinates.*/
class vect3D

Why not a struct? Since your data is public anyway

{
public:
 double x = 0.0, y = 0.0, z = 0.0;
 double tol = 1.0e-15;

What you're looking for is standardized here https://en.cppreference.com/w/cpp/types/numeric_limits/epsilon

 // Constructors
 vect3D() {};

Why would you define such a useless constructor? Just default it: vect3D() = default

By the way, there's a useless superfluous semi-colon at the end of all your function definitions.

 vect3D(std::initializer_list<double> RHS) { x = *RHS.begin(); y = *(RHS.begin()+1); z = *(RHS.begin()+2); };

If you conceive your vect3D as an aggregate (no user provided constructor, no virtual functions, and some other restrictions) you can use aggregate initialization directly (vect3D v{1., 2., 3.};) without having to define a brace-list constructor

 vect3D(std::vector<double> &RHS) { x = RHS[0]; y = RHS[1]; z = RHS[2]; };

So what if I provide a vector with one or two elements? Or even with six elements: that's clearly a mistake that won't be called out by the compiler.

 vect3D(double *RHS) { x = RHS[0]; y = RHS[1]; z = RHS[2]; };

Same here. The advice is to provide an interface that is hard to misuse

By the way, did you wonder why you got a lot of warnings about a deprecated implicit copy constructor? do you see why? Think of what will happen when the default copy constructor will copy your vector member by member, your vector of pointers included

 // Assignment
 vect3D& operator=(const vect3D &RHS) { x = RHS.x; y = RHS.y; z = RHS.z; return *this; };
 // Addition and subtraction
 vect3D& operator+=(const vect3D &RHS) { x += RHS.x; y += RHS.y; z += RHS.z; return *this; };
 vect3D& operator-=(const vect3D &RHS) { x -= RHS.x; y -= RHS.y; z -= RHS.z; return *this; };
 vect3D& operator+=(const double &RHS) { x += RHS; y += RHS; z += RHS; return *this; };
 vect3D& operator-=(const double &RHS) { x -= RHS; y -= RHS; z -= RHS; return *this; };
 vect3D operator+(const vect3D &RHS) { return vect3D(*this) += RHS; };
 vect3D operator-(const vect3D &RHS) { return vect3D(*this) -= RHS; };
 vect3D operator+(const double &RHS) { return vect3D(*this) += RHS; };
 vect3D operator-(const double &RHS) { return vect3D(*this) -= RHS; };
 // Scalar product and division
 vect3D& operator*=(const double &RHS) { x *= RHS; y *= RHS; z *= RHS; return *this; };
 vect3D& operator/=(const double &RHS) { x /= RHS; y /= RHS; z /= RHS; return *this; };
 vect3D operator*(const double &RHS) { return vect3D(*this) *= RHS; };
 vect3D operator/(const double &RHS) { return vect3D(*this) /= RHS; };
 friend vect3D operator*(double c, vect3D &vec) { return vec*c; };
 friend vect3D operator/(double c, vect3D &vec) { return vec/c; };
 // Comparisons
 bool operator==(const vect3D &RHS) { return ((x - RHS.x < x*tol) && (y - RHS.y < y*tol) && (z - RHS.z < z*tol)); };
 bool operator!=(const vect3D &RHS) { return !(*this == RHS); };
 bool operator>=(const vect3D &RHS) { return ((x >= RHS.x) && (y >= RHS.y) && (z >= RHS.z)); };
 bool operator<=(const vect3D &RHS) { return ((x <= RHS.x) && (y <= RHS.y) && (z <= RHS.z)); };
 bool operator>(const vect3D &RHS) { return !(*this <= RHS); };
 bool operator<(const vect3D &RHS) { return !(*this >= RHS); };
 // Euclidean norm
 double norm2() { return std::sqrt(std::pow(x, 2) + std::pow(y, 2) + std::pow(z, 2)); };
 friend double norm2(vect3D const &a) { return std::sqrt(std::pow(a.x, 2) + std::pow(a.y, 2) + std::pow(a.z, 2)); };

You don't need all those friends since x, y and z are public. That only clutters your class interface.

 // Dot product
 friend double dot(vect3D const &a, vect3D const &b) { return a.x*b.x + a.y*b.y + a.z*b.z; };
 // Cross product
 friend vect3D cross(vect3D const &a, vect3D const &b) { return {(a.y*b.z - a.z*b.y), (a.z*b.x - a.x*b.z), (a.x*b.y - a.y*b.x)}; };
 // Print to stream
 friend std::ostream& operator<<(std::ostream &stream, vect3D const &p) { return stream << "(" << p.x << ", " << p.y << ", " << p.z << ")" << std::flush; };
 // Function to explicitly return coordinates as an array of doubles, if needed
 void DoubleArray(double *v) { v[0] = x; v[1] = y; v[2] = z; return; };
 // To access coordinates using square bracket notation, for convenience
 std::vector<double *> p = std::vector<double *> {&x, &y, &z};

I'm amazed that someone with a sane mind (at least I hope so) would have such an extravagant idea. I'm not sure that square bracket notation is a convenience here -it might as well persuade the user that vect3D behaves like an array or a pointer- but even if you want to provide it, creating a vector of pointers is crazy, without even mentioning what will happen when a user naively calls for the third dimension with vec[3] (and they will). If you really want the bracket operator, implement it with a switch: case 0: return x; ...; default: throw std::out_of_bonds();

 double operator [] (int ii) const { return *(p[ii]); };
};
#endif