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!

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>

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