Attempt at templates by creating a class for N-dimensional mathematical vectors

Question

Originally I had tried to implement this topic after having learned about Inheritence and posted it on here, but that forced complexity and the main suggestion was that this was ideal for templates. Now that I've learnt about templates I was hoping someone could review my code for a basic implementation of mathematical vectors.

One thing I will note: There probably is a way to generalise the cross-product to N-dimensions but I wasn't too interested in finding out how since that's a more mathematical query, I just wanted to code so I only have cross product for three dimensions.

My main concerns are the uses of headers and source files, is it okay to have all of the implementation in one header file and the few necessary non-template functions defined in source files (like cross_product)?

And how should I treat friendship for cross_product? It only needs to be friends with Vector<3>, but I couldn't see a way of making exclusive Vector<3> friendship without using template specialisation, which seemed a bit too much. So it is friends with all Vector<N> but only accesses members of Vector<3>s.

Inside Vector.h:

#ifndef MATHSVectors
#define MATHSVectors

#include <string>
#include <iostream>
#include <initializer_list>
#include <stdexcept> //for the exceptions
#include <cmath> //sqrt, sin, cos, abs

//main template forward declaration
template <unsigned N> class Vector;
template <unsigned N> std::ostream& operator<< (std::ostream&, const Vector<N>&);
template <unsigned N> std::istream& operator>> (std::istream&, Vector<N>&);
template <unsigned N> bool operator==(const Vector<N>&, const Vector<N>&);
template <unsigned N> bool operator!=(const Vector<N>&, const Vector<N>&);
template <unsigned N> Vector<N> operator+(const Vector<N>&, const Vector<N>&);
template <unsigned N> Vector<N> operator-(const Vector<N>&, const Vector<N>&);
template <unsigned N> Vector<N> operator*(const Vector<N>&, double);
template <unsigned N> Vector<N> operator*(double, const Vector<N>&);
template <unsigned N> Vector<N> operator/(const Vector<N>&, double);
template <unsigned N> double dot_product(const Vector<N>&, const Vector<N>&);


//in Vector.cpp
bool double_equals(double, double);
double approximate(double, double);
Vector<3> cross_product(const Vector<3>&, const Vector<3>&);


template <unsigned N> class Vector{
    friend std::ostream& operator<< <N>(std::ostream&, const Vector&);
    friend std::istream& operator>> <N>(std::istream&, Vector&);
    friend bool operator== <N>(const Vector&, const Vector&);
    friend Vector operator+ <N>(const Vector&, const Vector&);
    friend Vector operator- <N>(const Vector&, const Vector&);
    friend Vector operator* <N>(const Vector&, double);
    friend Vector operator/ <N>(const Vector&, double);
    friend double dot_product<N>(const Vector&, const Vector&);
    friend Vector<3> cross_product(const Vector<3>&, const Vector<3>&);
public:
    Vector() = default;
    Vector(std::initializer_list<double>); //implicit conversion means we can assign from an initializer_list<double> 

    explicit operator bool() const;

    double& operator[](size_t p);
    const double& operator[](size_t p) const;

    Vector& operator+=(const Vector&);
    Vector& operator-=(const Vector&);
    Vector& operator*=(double);
    Vector& operator/=(double);

    double length() const;
    Vector& normalise();
    Vector& rotateCoordinates(size_t, size_t, double); //radians
private:
    double x[N] = {};
};





template <unsigned N> std::ostream& operator<<(std::ostream& os, const Vector<N>& rhs){
    os << "[";

    for(unsigned it = 0; it != N; ++it){
        os << approximate(rhs.x[it], 0);
        if(it != N-1) os << ", ";
    }

    os << "]";
    return os;
}

template <unsigned N> std::istream& operator>>(std::istream& is, Vector<N>& rhs){
    Vector<N> errorRet = rhs;

    for(unsigned it = 0; it != N; ++it)
        is >> rhs.x[it];

    if(!is)
        rhs = errorRet;

    return is;
}

template <unsigned N> bool operator==(const Vector<N>& lhs, const Vector<N>& rhs){
    for(unsigned it = 0; it != N; ++it)
        if(!double_equals(lhs.x[it], rhs.x[it])) return false;
    return true;
}

template <unsigned N> bool operator!=(const Vector<N>& lhs, const Vector<N>& rhs){
    return !(lhs == rhs);
}

template <unsigned N> Vector<N> operator+(const Vector<N>& lhs, const Vector<N>& rhs){
    Vector<N> sum = lhs;
    sum += rhs;
    return sum;
}

template <unsigned N> Vector<N> operator-(const Vector<N>& lhs, const Vector<N>& rhs){
    Vector<N> sum = lhs;
    sum -= rhs;
    return sum;
}

template <unsigned N> Vector<N> operator*(const Vector<N>& rhs, double d){
    Vector<N> product = rhs;
    product *= d;
    return product;
}

template <unsigned N> Vector<N> operator*(double d, const Vector<N>& rhs){
    return rhs*d;
}

template <unsigned N> Vector<N> operator/(const Vector<N>& rhs, double d){
    Vector<N> remain = rhs;
    remain /= d;
    return remain;
}

template <unsigned N> double dot_product(const Vector<N>& lhs, const Vector<N>& rhs) {
    double sum = 0;
    for(unsigned it = 0; it != N; ++it)
        sum += lhs.x[it]*rhs.x[it];
    return sum;
}

template <unsigned N> Vector<N>::Vector(std::initializer_list<double> li){
    if(N != li.size()) throw std::length_error("Attempt to initialise Vector with an initializer_list of different size.");

    for(unsigned it = 0; it != li.size(); ++it)
        x[it] = *(li.begin()+it);
}

template <unsigned N> Vector<N>::operator bool() const {
    return !(*this == Vector<N>());
}

template <unsigned N> double& Vector<N>::operator[](size_t p){
    if(p >= N) throw std::out_of_range(std::string("Invalid coordinate specified for function ") + __func__);
    return x[p];
}

template <unsigned N> const double& Vector<N>::operator[](size_t p) const {
    if(p >= N) throw std::out_of_range(std::string("Invalid coordinate specified for function ") + __func__);
    return x[p];
}

template <unsigned N> Vector<N>& Vector<N>::operator+=(const Vector<N>& rhs){
    for(unsigned it = 0; it != N; ++it)
        x[it] += rhs.x[it];
    return *this;
}

template <unsigned N> Vector<N>& Vector<N>::operator-=(const Vector<N>& rhs){
    for(unsigned it = 0; it != N; ++it)
        x[it] -= rhs.x[it];
    return *this;
}

template <unsigned N> Vector<N>& Vector<N>::operator*=(double d){
    for(unsigned it = 0; it != N; ++it)
        x[it] *= d;
    return *this;
}

template <unsigned N> Vector<N>& Vector<N>::operator/=(double d){
    if(d == 0) throw std::domain_error(std::string("Division by zero in function ") + __func__);
    for(unsigned it = 0; it != N; ++it)
        x[it] /= d;
    return *this;
}

template <unsigned N> double Vector<N>::length() const {
    double sum = 0;
    for(unsigned it = 0; it != N; ++it)
        sum += x[it]*x[it];
    return sqrt(sum);
}


template <unsigned N> Vector<N>& Vector<N>::normalise() {
    if(!(*this)) return (*this); //null vector
    return (*this)/=length();
}

template <unsigned N> Vector<N>& Vector<N>::rotateCoordinates(size_t i, size_t j, double angle){
    if(i >= N || j >= N) throw std::out_of_range(std::string("Invalid coordinate specified for function ") + __func__);

    double newI = x[i]*cos(angle) - x[j]*sin(angle);
    double newJ = x[j]*cos(angle) + x[i]*sin(angle);
    x[i] = newI;
    x[j] = newJ;
    return *this;
}


#endif

Inside Vector.cpp:

/* A few notes about dealing with doubles

    1 i. The double compare function needs to be changed. It is not a transitive equality operation. A method to
         fix this would be to snap the doubles on to a grid and return true if two doubles snap on to the same section
         of the grid. I don't know how to implement this just yet - wait until I've read more about float comparisons.
    2. Doubles can get stored as negative zero, so adding +.0 when outputting the vector prevents displaying "-0".

*/

#include "Vector.h"

using namespace std;



const double epsilon = 1e-6; //double tolerance

bool double_equals(double a, double b){
    return abs(a-b) < epsilon;
}

double approximate(double a, double b){
    return double_equals(a,b) ? b : a;
}

Vector<3> cross_product(const Vector<3>& lhs, const Vector<3>& rhs){
    double newX = (lhs.x[1]*rhs.x[2]) - (lhs.x[2]*rhs.x[1]);
    double newY = (lhs.x[2]*rhs.x[0]) - (lhs.x[0]*rhs.x[2]);
    double newZ = (lhs.x[0]*rhs.x[1]) - (lhs.x[1]*rhs.x[0]);

    return Vector<3>({newX, newY, newZ});

}

Answer 1

This looks pretty good. I don't have any deal-breaker comments, just a bunch of minor ones.

Prefer std::array

Just straight up std::array<double, N> x instead of double x[N]. Raw arrays are broken. There's not really any real difference for the purposes of your use-case, but it's just a more pleasant type to deal with in general.

Throwing and non-throwing

For functions that throw, there's an extra mechanism that needs to exist to support that use-case. Plus even an always-predicted branch is going to be more code that no branch. To that end the standard containers offer throwing and non-throwing functions. If you rewrite your two operator[]s to be non-throwing:

template <unsigned N>
double& Vector<N>::operator[](size_t p) {
    return x[p];
}

And introduce a throwing alternative:

template <unsigned N>
double& Vector<N>::at(size_t p) {
    if(p >= N) throw std::out_of_range(std::string("Invalid coordinate specified for function ") + __func__); 
    return (*this)[p];
}

That will give you two benefits. First, users of your class would expect operator[] to not throw since that's pretty typical. But second, you no longer need to friend any of the free functions. cross_product only needed to be a friend because you wanted to avoid the throwing, now that's just a non-issue.

Use some standard algorithms

Your equality operator can be reimplemented with std::equal:

template <unsigned N> 
bool operator==(const Vector<N>& lhs, const Vector<N>& rhs) {
    return std::equal(lhs.x.begin(), lhs.x.end(), rhs.x.begin(), double_equals);
}

Similarly, operator bool doesn't need to construct a whole new vector, just check the one you have against zero with std::any_of:

template <unsigned N>
Vector<N>::operator bool() const {
    return std::any_of(lhs.x.begin(), lhs.x.end(), [](double v){
        return !double_equals(v, 0.0);
    });
}

Initialization from std::initializer_list<double> can use std::copy:

 std::copy(x.begin(), x.end(), li.begin());

Normalize

I suspect the typical case here might be for valid vectors, so prefer to take the length first and compare that against zero, rather than invoking operator!:

double magnitude = length();
if (magnitude > 0) {
    (*this) /= magnitude;
}
return *this;

I also question the name length(). That seems to suggest that this is a container more along the lines of std::vector, but really you're talking about the norm of the vector. So a better name, please.

Use range-based for

Whenever you loop in a way that doesn't actually involve the index, prefer to use range-based for. It's just easier to write. For instance, your operator*=:

for (double& val : x) {
    val *= d;
}
return *this;

Friend Operators

Rather than forward declaring template functions, then friending them, prefer to write these operators as non-member non-template friends. So operator+ would be:

template <unsigned N>
class Vector {
public:
    friend Vector operator+(Vector const& lhs, Vector const& rhs) {
        Vector sum = lhs;
        sum += rhs;
        return sum;
    }
};

This will be found by ADL and nothing else, and then you also don't have to worry about some of the other issues that writing function templates leads you to. It's just simpler.

rotateCoordinates

This seems like it should be a non-member to me.

Answer 2

I have a similar class for physical vectors. In addition to the number of elements, it alse parameterizes element type, which must be an arithmetic type (controlled with static_assert). In addition to the members/friends you defined, I have the following (all implemented such that they are found only by ADL)

element_type operator*(vector const&, vector const&);  // dot(x,y)
element_type norm(vector const&);                      // dot(x,x)
element_type abs(vector const&);                       // sqrt(norm(x))
element_type dist_sq(vector const&, vector const&);    // norm(x-y)
element_type distance(vector const&, vector const&);   // sqrt(dist_sq(x,y))
vector operator^(vector<3> const&, vector<3> const&);  // cross_product
element_type operator^(vector<2> const&, vector<2> const&); // z-component of cross product

and several more non-arithmetic operations, as well as an iterator interface a la std::array. I find overloading the * and ^ operators as vector dot and cross product useful (though the latter has lower precedence than all the arithmetic operators +, -, *, /, so some care is required). One can go one and define a corresponding matrix (or tensor) class and its interoperability with vectors.