An N-dimensional mathematical vector template class deriving from std::array (C++-17)

Question

I've coded up a Vector class for use in my simulation support code, which I offer up for review.

I decided to extend std::array. Now, I know that std::array is an aggregate type, but with C++-17 we can now extend an aggregate base type and still use 'curly brace' initialization. Extending std::array means that each Vector consumes the same, fixed amount of memory as a regular std::array. I've added methods for arithmetic and various Vector operations, such as the cross-product (for N=3 dimensions only), the dot product, and length.

I'd welcome opinions over whether this is a good idea in principle, and also specific suggestions about the details of the code.

The full source code follows. Note that it makes use of a randomization class which is outside the scope of this question, but if you want to see it it is here: https://github.com/ABRG-Models/morphologica/blob/master/morph/Random.h

/*!
 * \file
 * \brief An N dimensional vector class template which derives from std::array.
 *
 * \author Seb James
 * \date April 2020
 */
#pragma once

#include <cmath>
using std::abs;
using std::sqrt;

#include <array>
using std::array;

#include <iostream>
using std::cout;
using std::endl;
using std::ostream;

#include <string>
using std::string;
#include <sstream>
using std::stringstream;

#include <type_traits>
using std::enable_if;
using std::enable_if_t;
using std::is_integral;
using std::is_scalar;
using std::decay_t;

#include "Random.h"
using morph::RandUniformReal;
using morph::RandUniformInt;

namespace morph {

    /*!
     * \brief N-D vector class
     *
     * An N dimensional vector class template which derives from std::array. Vector
     * components are of scalar type S. It is anticipated that S will be set either to
     * floating point scalar types such as float or double, or to integer scalar types
     * such as int, long long int and so on. Thus, a typical (and in fact, the default)
     * signature would be:
     *
     * Vector<float, 3> v;
     *
     * The class inherits std:array's fixed-size array of memory for storing the
     * components of the vector. It adds numerous methods which allow objects of type
     * Vector to have arithmetic operations applied to them, either scalar (add a scalar
     * to all elements; divide all elements by a scalar, etc) or vector (including dot
     * and cross products, normalization and so on.
     *
     * Because morph::Vector extends std::array, it works best when compiled with a
     * c++-17 compiler (although it can be compiled with a c++-11 compiler). This is
     * because std::array is an 'aggregate class' with no user-provided constructors,
     * and morph::Vector does not add any of its own constructors. Prior to c++-17,
     * aggregate classes were not permitted to have base classes. So, if you want to do:
     *
     * Vector<float, 3> v = { 1.0f , 1.0f, 1.0f };
     *
     * You need c++-17. Otherwise, restrict your client code to doing:
     *
     * Vector<float, 3> v;
     * v[0] = 1.0f; v[1] = 1.0f; v[2] = 1.0f;
     */
    template <typename S, size_t N> struct Vector;

    /*!
     * Template friendly mechanism to overload the stream operator.
     *
     * Note forward declaration of the Vector template class and this template for
     * stream operator overloading. Example adapted from
     * https://stackoverflow.com/questions/4660123
     */
    template <typename S, size_t N> ostream& operator<< (ostream&, const Vector<S, N>&);

    template <typename S=float, size_t N=3>
    struct Vector : public array<S, N>
    {
        //! \return the first component of the vector
        template <size_t _N = N, enable_if_t<(_N>0), int> = 0>
        S x (void) const {
            return (*this)[0];
        }
        //! \return the second component of the vector
        template <size_t _N = N, enable_if_t<(_N>1), int> = 0>
        S y (void) const {
            return (*this)[1];
        }
        //! \return the third component of the vector
        template <size_t _N = N, enable_if_t<(_N>2), int> = 0>
        S z (void) const {
            return (*this)[2];
        }
        //! \return the fourth component of the vector
        template <size_t _N = N, enable_if_t<(_N>3), int> = 0>
        S w (void) const {
            return (*this)[3];
        }

        /*!
         * \brief Unit vector threshold
         *
         * The threshold outside of which the vector is no longer considered to be a
         * unit vector. Note this is hard coded as a constexpr, to avoid messing with
         * the initialization of the Vector with curly brace initialization.
         *
         * Clearly, this will be the wrong threshold for some cases. Possibly, a
         * template parameter could set this; so size_t U could indicate the threshold;
         * 0.001 could be U=-3 (10^-3).
         *
         * Another idea would be to change unitThresh based on the type S. Or use
         * numeric_limits<S>::epsilon and find out what multiple of epsilon would make
         * sense.
         */
        static constexpr S unitThresh = 0.001;

        /*!
         * Set data members from an array the of same size and type.
         */
        void setFrom (const array<S, N> v) {
            for (size_t i = 0; i < N; ++i) {
                (*this)[i] = v[i];
            }
        }

        /*!
         * Set the data members of this Vector from the passed in, larger vector, v,
         * ignoring the last element of v. Used when working with 4D vectors in graphics
         * applications involving 4x4 transform matrices.
         */
        void setFrom (const array<S, (N+1)> v) {
            for (size_t i = 0; i < N; ++i) {
                (*this)[i] = v[i];
            }
        }

        /*!
         * Output the vector to stdout
         */
        void output (void) const {
            cout << "Vector" << this->asString();
        }

        /*!
         * Create a string representation of the vector
         *
         * \return A 'coordinate format' string such as "(1,1,2)", "(0.2,0.4)" or
         * "(5,4,5,5,40)".
         */
        string asString (void) const {
            stringstream ss;
            auto i = this->begin();
            ss << "(";
            bool first = true;
            while (i != this->end()) {
                if (first) {
                    ss << *i++;
                    first = false;
                } else {
                    ss << "," << *i++;
                }
            }
            ss << ")";
            return ss.str();
        }

        /*!
         * Renormalize the vector to length 1.
         */
        void renormalize (void) {
            S denom = static_cast<S>(0);
            auto i = this->begin();
            while (i != this->end()) {
                denom += ((*i) * (*i));
                ++i;
            }
            denom = sqrt(denom);
            if (denom != static_cast<S>(0.0)) {
                S oneovermag = static_cast<S>(1.0) / denom;
                i = this->begin();
                while (i != this->end()) {
                    *i++ *= oneovermag;
                }
            }
        }

        /*!
         * Randomize the vector
         *
         * Randomly set the elements of the vector consisting of floating point
         * coordinates. Coordinates are set to random numbers drawn from a uniform
         * distribution between 0 and 1 (See morph::RandUniformReal for details).
         *
         * Note that I need a real or int implementation here, depending on the type of
         * S. This allows me to use the correct type of randomizer.
         *
         * Note, if you omit the second template arg from enable_if_t (or enable_if)
         * then the type defaults to void.
         *
         * \tparam F A floating point scalar type
         */
        template <typename F=S, enable_if_t<!is_integral<decay_t<F>>::value, int> = 0 >
        void randomize (void) {
            RandUniformReal<F> ruf (static_cast<F>(0), static_cast<F>(1));
            auto i = this->begin();
            while (i != this->end()) {
                *i++ = ruf.get();
            }
        }

        /*!
         * Randomize the vector
         *
         * Randomly set the elements of the vector consisting of integer
         * coordinates. Coordinates are set to random numbers drawn from a uniform
         * distribution between 0 and 255 (See morph::RandUniformInt for details).
         *
         * Note on the template syntax: Here, if I is integral, then enable_if_t's type
         * is '0' and the function is defined (I think).
         *
         * \tparam I An integer scalar type
         */
        template <typename I=S, enable_if_t<is_integral<decay_t<I>>::value, int> = 0 >
        void randomize (void) {
            RandUniformInt<I> rui (static_cast<I>(0), static_cast<I>(255));
            auto i = this->begin();
            while (i != this->end()) {
                *i++ = rui.get();
            }
        }

        /*!
         * Test to see if this vector is a unit vector (it doesn't *have* to be).
         *
         * \return true if the length of the vector is 1.
         */
        bool checkunit (void) const {
            bool rtn = true;
            S metric = 1.0;
            auto i = this->begin();
            while (i != this->end()) {
                metric -= ((*i) * (*i));
                ++i;
            }
            if (abs(metric) > morph::Vector<S, N>::unitThresh) {
                rtn = false;
            }
            return rtn;
        }

        /*!
         * Find the length of the vector.
         *
         * \return the length
         */
        S length (void) const {
            S sos = static_cast<S>(0);
            auto i = this->begin();
            while (i != this->end()) {
                sos += ((*i) * (*i));
                ++i;
            }
            return sqrt(sos);
        }

        /*!
         * Unary negate operator
         *
         * \return a Vector whose elements have been negated.
         */
        Vector<S, N> operator- (void) const {
            Vector<S, N> rtn;
            auto i = this->begin();
            auto j = rtn.begin();
            while (i != this->end()) {
                *j++ = -(*i++);
            }
            return rtn;
        }

        /*!
         * Unary not operator.
         *
         * \return true if the vector length is 0, otherwise it returns false.
         */
        bool operator! (void) const {
            return (this->length() == static_cast<S>(0.0)) ? true : false;
        }

        /*!
         * Vector multiply * operator.
         *
         * Cross product of this with another vector v2 (if N==3). In
         * higher dimensions, its more complicated to define what the cross product is,
         * and I'm unlikely to need anything other than the plain old 3D cross product.
         */
        template <size_t _N = N, enable_if_t<(_N==3), int> = 0>
        Vector<S, N> operator* (const Vector<S, _N>& v2) const {
            Vector<S, _N> v;
            v[0] = (*this)[1] * v2.z() - (*this)[2] * v2.y();
            v[1] = (*this)[2] * v2.x() - (*this)[0] * v2.z();
            v[2] = (*this)[0] * v2.y() - (*this)[1] * v2.x();
            return v;
        }

        /*!
         * Vector multiply *= operator.
         *
         * Cross product of this with another vector v2 (if N==3). Result written into
         * this.
         */
        template <size_t _N = N, enable_if_t<(_N==3), int> = 0>
        void operator*= (const Vector<S, _N>& v2) {
            Vector<S, _N> v;
            v[0] = (*this)[1] * v2.z() - (*this)[2] * v2.y();
            v[1] = (*this)[2] * v2.x() - (*this)[0] * v2.z();
            v[2] = (*this)[0] * v2.y() - (*this)[1] * v2.x();
            (*this)[0] = v[0];
            (*this)[1] = v[1];
            (*this)[2] = v[2];
        }

        /*!
         * \brief Scalar (dot) product
         *
         * Compute the scalar product of this Vector and the Vector, v2.
         *
         * \return scalar product
         */
        S dot (const Vector<S, N>& v2) const {
            S rtn = static_cast<S>(0);
            auto i = this->begin();
            auto j = v2.begin();
            while (i != this->end()) {
                rtn += ((*i++) * (*j++));
            }
            return rtn;
        }

        /*!
         * Scalar multiply * operator
         *
         * This function will only be defined if typename _S is a
         * scalar type. Multiplies this Vector<S, N> by s, element-wise.
         */
        template <typename _S=S, enable_if_t<is_scalar<decay_t<_S>>::value, int> = 0 >
        Vector<S, N> operator* (const _S& s) const {
            Vector<S, N> rtn;
            auto val = this->begin();
            auto rval = rtn.begin();
            // Here's a way to iterate through which the compiler should be able to
            // autovectorise; it knows what i is on each loop:
            for (size_t i = 0; i < N; ++i) {
                *(rval+i) = *(val+i) * static_cast<S>(s);
            }
            return rtn;
        }

        /*!
         * Scalar multiply *= operator
         *
         * This function will only be defined if typename _S is a
         * scalar type. Multiplies this Vector<S, N> by s, element-wise.
         */
        template <typename _S=S, enable_if_t<is_scalar<decay_t<_S>>::value, int> = 0 >
        void operator*= (const _S& s) {
            auto val = this->begin();
            for (size_t i = 0; i < N; ++i) {
                *(val+i) *= static_cast<S>(s);
            }
        }

        /*!
         * Scalar division * operator
         */
        template <typename _S=S, enable_if_t<is_scalar<decay_t<_S>>::value, int> = 0 >
        Vector<S, N> operator/ (const _S& s) const {
            Vector<S, N> rtn;
            auto val = this->begin();
            auto rval = rtn.begin();
            for (size_t i = 0; i < N; ++i) {
                *(rval+i) = *(val+i) / static_cast<S>(s);
            }
            return rtn;
        }

        /*!
         * Scalar division *= operator
         */
        template <typename _S=S, enable_if_t<is_scalar<decay_t<_S>>::value, int> = 0 >
        void operator/= (const _S& s) {
            auto val = this->begin();
            for (size_t i = 0; i < N; ++i) {
                *(val+i) /= static_cast<S>(s);
            }
        }

        /*!
         * Vector addition operator
         */
        Vector<S, N> operator+ (const Vector<S, N>& v2) const {
            Vector<S, N> v;
            auto val = this->begin();
            auto val2 = v2.begin();
            for (size_t i = 0; i < N; ++i) {
                v[i] = *(val+i) + *(val2+i);
            }
            return v;
        }

        /*!
         * Vector addition operator
         */
        void operator+= (const Vector<S, N>& v2) {
            auto val = this->begin();
            auto val2 = v2.begin();
            for (size_t i = 0; i < N; ++i) {
                *(val+i) += *(val2+i);
            }
        }

        /*!
         * Vector subtraction
         */
        Vector<S, N> operator- (const Vector<S, N>& v2) const {
            Vector<S, N> v;
            auto val = this->begin();
            auto val2 = v2.begin();
            for (size_t i = 0; i < N; ++i) {
                v[i] = *(val+i) - *(val2+i);
            }
            return v;
        }

        /*!
         * Vector subtraction
         */
        void operator-= (const Vector<S, N>& v2) {
            auto val = this->begin();
            auto val2 = v2.begin();
            for (size_t i = 0; i < N; ++i) {
                *(val+i) -= *(val2+i);
            }
        }

        /*!
         * Scalar addition
         */
        template <typename _S=S, enable_if_t<is_scalar<decay_t<_S>>::value, int> = 0 >
        Vector<S, N> operator+ (const _S& s) const {
            Vector<S, N> rtn;
            auto val = this->begin();
            auto rval = rtn.begin();
            for (size_t i = 0; i < N; ++i) {
                *(rval+i) = *(val+i) + static_cast<S>(s);
            }
            return rtn;
        }

        /*!
         * Scalar addition
         */
        template <typename _S=S, enable_if_t<is_scalar<decay_t<_S>>::value, int> = 0 >
        void operator+= (const _S& s) {
            auto val = this->begin();
            for (size_t i = 0; i < N; ++i) {
                *(val+i) += static_cast<S>(s);
            }
        }

        /*!
         * Scalar subtraction
         */
        template <typename _S=S, enable_if_t<is_scalar<decay_t<_S>>::value, int> = 0 >
        Vector<S, N> operator- (const _S& s) const {
            Vector<S, N> rtn;
            auto val = this->begin();
            auto rval = rtn.begin();
            for (size_t i = 0; i < N; ++i) {
                *(rval+i) = *(val+i) - static_cast<S>(s);
            }
            return rtn;
        }

        /*!
         * Scalar subtraction
         */
        template <typename _S=S, enable_if_t<is_scalar<decay_t<_S>>::value, int> = 0 >
        void operator-= (const _S& s) {
            auto val = this->begin();
            for (size_t i = 0; i < N; ++i) {
                *(val+i) -= static_cast<S>(s);
            }
        }

        /*!
         * Overload the stream output operator
         */
        friend ostream& operator<< <S, N> (ostream& os, const Vector<S, N>& v);
    };

    template <typename S=float, size_t N=3>
    ostream& operator<< (ostream& os, const Vector<S, N>& v)
    {
        os << v.asString();
        return os;
    }

} // namespace morph
```

Answer 1

There’s a lot of going on here. At a cursory glance:

• You must not use using in headers at file scope.
• It’s not clear why your template parameter is called S instead of T, and the more I read your source code the more this confuses me. What does S stand for? — The use of T for such member types is near-universal.
• I’m not sure public inheritance is a good idea here: why are you modelling an is-a relationship? Why not private inheritance or a std::array member?
• Don’t declare functions with a (void) parameter list. In C this is necessary to create the correct prototype. In C++ it has no purpose — () does the same, and is conventional.
• Your setFrom member functions should be constructors and assignment operators instead.
• Don’t pass std::array by value, pass it by const& — otherwise a very expensive copy might be created. At the very least make this depend on N as well as sizeof(S) so you can optimise for arrays small enough to be passed inside a single register.
• Use algorithms (std::copy, constructors, assignment) instead of copying your arrays in a for loop.
• output is redundant if you define a suitable formatted output stream operator.
• If you want to follow C++ convention, asString should be called str. That’s not necessary of course.
• S denom = static_cast<S>(0); can usually be written as auto denom = S{0};, and those cases where this fails because no suitable constructor exists are probably cases where you want this to fail.
• Don’t use while loops to iterate over ranges, it’s unidiomatic and therefore confusing: either use for loops or, better yet, ranged-for loops where possible. And once again use appropriate algorithms. The loop that calculates denom can be replaced by a call to std::accumulate, for instance. That way you can also declare denom as const and initialise it directly.
• randomize guards against S being an integral type; renomalize does not, but also needs this constraint.

Answer 2

Many thanks to Konrad Rudolph and Miguel Avila for time spent reading the code and responding.

I'm responding in an answer, because I want to post an update to the code in which I've followed most of Konrad's suggestions.

In answer to Konrad's points:

• Thank you for pointing out the error with using in headers. Fixed.
• I named the template parameter 'S' for Scalar. The idea was to hint that, although this derives from std::array, it is not to be treated as a generic container of any other type. I realise this might look confusing, when T is conventional, but I have left the naming scheme in my updated code as I think this is perhaps a matter of preference rather than correctness.
• I did want an is-a relationship. This relates to Miguel's answer too. I will discuss after finishing this list.
• fn(void) all changed to fn(); I'm happy to follow this convention; it looks better, too.
• setFrom(): First, I changed to set_from() in the update because lower case looks more in-keeping with STL derived code. Second it can't be a constructor, because I want to keep the Vector as an aggregate type. I don't anticipate using set_from() often, but there is some of my code which needs the feature.
• Good point. I've updated and made use of algorithms; this was a good exercise for me.
• Good point; I've removed output()
• Good point; asString() is now str()
• I like this alternative to static_casts, so I've updated the code accordingly.
• Many of the while's disappeared with std::aggregate and std::transform. Others have been replaced with fors
• I've added a guard for renormalise; I think that occurred to me while I was writing it, but it got forgotten!

So, the last point to discuss is whether deriving from std::array is a good idea. You guys are suspicious, I think, but I still like it.

What does it do for the Vector class?

• No need to write constructors - they're banned anyway
• compiler knows exactly the size of the object
• I can use STL iterators on my Vector objects in client code. Consider that these vectors may be used for computations where they really do have high dimensionality (neural nets are one application)

Downsides

• If you really need to assign from another type, such as std::array, you have to write an ugly method like set_from().
• Because it derives from an STL container, a coder using Vector might think they can set S=array_type or S=SomeClass - but the type_trait tests should address that.

If you can think of any other concrete downsides, do please list them!

Here is the updated class template. Note that I've included the Random number generator code so that this listing is compilable, even though that code is out of scope. I'll also list a test program to compile against it.

Thanks again!

/*!
 * \file
 * \brief An N dimensional vector class template which derives from std::array.
 *
 * \author Seb James (with thanks to Konrad Rudolph  and Miguel Avila for code review)
 * \date April 2020
 */
#pragma once

#include <cmath>
#include <array>
#include <iostream>
#include <string>
#include <sstream>
#include <type_traits>
#include <numeric>
#include <algorithm>
#include <functional>
// For random number generator only (would normally be a separate <Random.h>)
#include <random>
#include <vector>
#include <limits>

namespace morph {

    /*!
     * \brief N-D vector class
     *
     * An N dimensional vector class template which derives from std::array. Vector
     * components are of scalar type S. It is anticipated that S will be set either to
     * floating point scalar types such as float or double, or to integer scalar types
     * such as int, long long int and so on. Thus, a typical (and in fact, the default)
     * signature would be:
     *
     * Vector<float, 3> v;
     *
     * The class inherits std:array's fixed-size array of memory for storing the
     * components of the vector. It adds numerous methods which allow objects of type
     * Vector to have arithmetic operations applied to them, either scalar (add a scalar
     * to all elements; divide all elements by a scalar, etc) or vector (including dot
     * and cross products, normalization and so on.
     *
     * Because morph::Vector extends std::array, it works best when compiled with a
     * c++-17 compiler (although it can be compiled with a c++-11 compiler). This is
     * because std::array is an 'aggregate class' with no user-provided constructors,
     * and morph::Vector does not add any of its own constructors. Prior to c++-17,
     * aggregate classes were not permitted to have base classes. So, if you want to do:
     *
     * Vector<float, 3> v = { 1.0f , 1.0f, 1.0f };
     *
     * You need c++-17. Otherwise, restrict your client code to doing:
     *
     * Vector<float, 3> v;
     * v[0] = 1.0f; v[1] = 1.0f; v[2] = 1.0f;
     */
    template <typename S, size_t N> struct Vector;

    /*!
     * Template friendly mechanism to overload the stream operator.
     *
     * Note forward declaration of the Vector template class and this template for
     * stream operator overloading. Example adapted from
     * https://stackoverflow.com/questions/4660123
     */
    template <typename S, size_t N> std::ostream& operator<< (std::ostream&, const Vector<S, N>&);

    //! forward declaration of random number generator classes
    template <typename T> class RandUniformReal;
    template <typename T> class RandUniformInt;

    template <typename S=float, size_t N=3>
    struct Vector : public std::array<S, N>
    {
        //! \return the first component of the vector
        template <size_t _N = N, std::enable_if_t<(_N>0), int> = 0>
        S x() const {
            return (*this)[0];
        }
        //! \return the second component of the vector
        template <size_t _N = N, std::enable_if_t<(_N>1), int> = 0>
        S y() const {
            return (*this)[1];
        }
        //! \return the third component of the vector
        template <size_t _N = N, std::enable_if_t<(_N>2), int> = 0>
        S z() const {
            return (*this)[2];
        }
        //! \return the fourth component of the vector
        template <size_t _N = N, std::enable_if_t<(_N>3), int> = 0>
        S w() const {
            return (*this)[3];
        }

        /*!
         * \brief Unit vector threshold
         *
         * The threshold outside of which the vector is no longer considered to be a
         * unit vector. Note this is hard coded as a constexpr, to avoid messing with
         * the initialization of the Vector with curly brace initialization.
         *
         * Clearly, this will be the wrong threshold for some cases. Possibly, a
         * template parameter could set this; so size_t U could indicate the threshold;
         * 0.001 could be U=-3 (10^-3).
         *
         * Another idea would be to change unitThresh based on the type S. Or use
         * numeric_limits<S>::epsilon and find out what multiple of epsilon would make
         * sense.
         */
        static constexpr S unitThresh = 0.001;

        /*!
         * Set data members from an array the of same size and type.
         */
        void set_from (const std::array<S, N>& ar) {
            std::copy (ar.begin(), ar.end(), this->begin());
        }

        /*!
         * Set the data members of this Vector from the passed in, larger array, \a ar,
         * ignoring the last element of \a ar. Used when working with 4D vectors in
         * graphics applications involving 4x4 transform matrices.
         */
        void set_from (const std::array<S, (N+1)>& ar) {
            // Don't use std::copy here, because ar has more elements than *this.
            for (size_t i = 0; i < N; ++i) {
                (*this)[i] = ar[i];
            }
        }

        /*!
         * Set an N-D Vector from an N+1 D Vector. Intended to convert 4D vectors (that
         * have been operated on by 4x4 matrices) into 3D vectors.
         */
        void set_from (const Vector<S, (N+1)>& v) {
            for (size_t i = 0; i < N; ++i) {
                (*this)[i] = v[i];
            }
        }

        /*!
         * Create a string representation of the vector
         *
         * \return A 'coordinate format' string such as "(1,1,2)", "(0.2,0.4)" or
         * "(5,4,5,5,40)".
         */
        std::string str() const {
            std::stringstream ss;
            ss << "(";
            bool first = true;
            for (auto i : *this) {
                if (first) {
                    ss << i;
                    first = false;
                } else {
                    ss << "," << i;
                }
            }
            ss << ")";
            return ss.str();
        }

        /*!
         * Renormalize the vector to length 1.0. Only for S types that are floating point.
         */
        template <typename F=S, std::enable_if_t<!std::is_integral<std::decay_t<F>>::value, int> = 0 >
        void renormalize() {
            auto add_squared = [](F a, F b) { return a + b * b; };
            const F denom = std::sqrt (std::accumulate (this->begin(), this->end(), F{0}, add_squared));
            if (denom != F{0}) {
                F oneovermag = F{1} / denom;
                auto x_oneovermag = [oneovermag](F f) { return f * oneovermag; };
                std::transform (this->begin(), this->end(), this->begin(), x_oneovermag);
            }
        }

        /*!
         * Randomize the vector
         *
         * Randomly set the elements of the vector consisting of floating point
         * coordinates. Coordinates are set to random numbers drawn from a uniform
         * distribution between 0 and 1 (See morph::RandUniformReal for details).
         *
         * Note that I need a real or int implementation here, depending on the type of
         * S. This allows me to use the correct type of randomizer.
         *
         * Note, if you omit the second template arg from enable_if_t (or enable_if)
         * then the type defaults to void.
         *
         * \tparam F A floating point scalar type
         */
        template <typename F=S, std::enable_if_t<!std::is_integral<std::decay_t<F>>::value, int> = 0 >
        void randomize() {
            RandUniformReal<F> ruf (static_cast<F>(0), static_cast<F>(1));
            for (auto& i : *this) {
                i = ruf.get();
            }
        }

        /*!
         * Randomize the vector
         *
         * Randomly set the elements of the vector consisting of integer
         * coordinates. Coordinates are set to random numbers drawn from a uniform
         * distribution between 0 and 255 (See morph::RandUniformInt for details).
         *
         * Note on the template syntax: Here, if I is integral, then enable_if_t's type
         * is '0' and the function is defined (I think).
         *
         * \tparam I An integer scalar type
         */
        template <typename I=S, std::enable_if_t<std::is_integral<std::decay_t<I>>::value, int> = 0 >
        void randomize() {
            RandUniformInt<I> rui (static_cast<I>(0), static_cast<I>(255));
            for (auto& i : *this) {
                i = rui.get();
            }
        }

        /*!
         * Test to see if this vector is a unit vector (it doesn't *have* to be).
         *
         * \return true if the length of the vector is 1.
         */
        bool checkunit() const {
            auto subtract_squared = [](S a, S b) { return a - b * b; };
            const S metric = std::accumulate (this->begin(), this->end(), S{1}, subtract_squared);
            if (std::abs(metric) > Vector<S, N>::unitThresh) {
                return false;
            }
            return true;
        }

        /*!
         * Find the length of the vector.
         *
         * \return the length
         */
        S length() const {
            auto add_squared = [](S a, S b) { return a + b * b; };
            const S len = std::sqrt (std::accumulate (this->begin(), this->end(), S{0}, add_squared));
            return len;
        }

        /*!
         * Unary negate operator
         *
         * \return a Vector whose elements have been negated.
         */
        Vector<S, N> operator-() const {
            Vector<S, N> rtn;
            std::transform (this->begin(), this->end(), rtn.begin(), std::negate<S>());
            return rtn;
        }

        /*!
         * Unary not operator.
         *
         * \return true if the vector length is 0, otherwise it returns false.
         */
        bool operator!() const {
            return (this->length() == S{0}) ? true : false;
        }

        /*!
         * Vector multiply * operator.
         *
         * Cross product of this with another vector \a v (if N==3). In
         * higher dimensions, its more complicated to define what the cross product is,
         * and I'm unlikely to need anything other than the plain old 3D cross product.
         */
        template <size_t _N = N, std::enable_if_t<(_N==3), int> = 0>
        Vector<S, N> operator* (const Vector<S, _N>& v) const {
            Vector<S, _N> vrtn;
            vrtn[0] = (*this)[1] * v.z() - (*this)[2] * v.y();
            vrtn[1] = (*this)[2] * v.x() - (*this)[0] * v.z();
            vrtn[2] = (*this)[0] * v.y() - (*this)[1] * v.x();
            return vrtn;
        }

        /*!
         * Vector multiply *= operator.
         *
         * Cross product of this with another vector v (if N==3). Result written into
         * this.
         */
        template <size_t _N = N, std::enable_if_t<(_N==3), int> = 0>
        void operator*= (const Vector<S, _N>& v) {
            Vector<S, _N> vtmp;
            vtmp[0] = (*this)[1] * v.z() - (*this)[2] * v.y();
            vtmp[1] = (*this)[2] * v.x() - (*this)[0] * v.z();
            vtmp[2] = (*this)[0] * v.y() - (*this)[1] * v.x();
            std::copy (vtmp.begin(), vtmp.end(), this->begin());
        }

        /*!
         * \brief Scalar (dot) product
         *
         * Compute the scalar product of this Vector and the Vector, v.
         *
         * \return scalar product
         */
        S dot (const Vector<S, N>& v) const {
            auto vi = v.begin();
            auto dot_product = [vi](S a, S b) mutable { return a + b * (*vi++); };
            const S rtn = std::accumulate (this->begin(), this->end(), S{0}, dot_product);
            return rtn;
        }

        /*!
         * Scalar multiply * operator
         *
         * This function will only be defined if typename _S is a
         * scalar type. Multiplies this Vector<S, N> by s, element-wise.
         */
        template <typename _S=S, std::enable_if_t<std::is_scalar<std::decay_t<_S>>::value, int> = 0 >
        Vector<_S, N> operator* (const _S& s) const {
            Vector<_S, N> rtn;
            auto mult_by_s = [s](_S coord) { return coord * s; };
            std::transform (this->begin(), this->end(), rtn.begin(), mult_by_s);
            return rtn;
        }

        /*!
         * Scalar multiply *= operator
         *
         * This function will only be defined if typename _S is a
         * scalar type. Multiplies this Vector<S, N> by s, element-wise.
         */
        template <typename _S=S, std::enable_if_t<std::is_scalar<std::decay_t<_S>>::value, int> = 0 >
        void operator*= (const _S& s) {
            auto mult_by_s = [s](_S coord) { return coord * s; };
            std::transform (this->begin(), this->end(), this->begin(), mult_by_s);
        }

        /*!
         * Scalar divide by s
         */
        template <typename _S=S, std::enable_if_t<std::is_scalar<std::decay_t<_S>>::value, int> = 0 >
        Vector<_S, N> operator/ (const _S& s) const {
            Vector<_S, N> rtn;
            auto div_by_s = [s](_S coord) { return coord / s; };
            std::transform (this->begin(), this->end(), rtn.begin(), div_by_s);
            return rtn;
        }

        /*!
         * Scalar divide by s
         */
        template <typename _S=S, std::enable_if_t<std::is_scalar<std::decay_t<_S>>::value, int> = 0 >
        void operator/= (const _S& s) {
            auto div_by_s = [s](_S coord) { return coord / s; };
            std::transform (this->begin(), this->end(), this->begin(), div_by_s);
        }

        /*!
         * Vector addition operator
         */
        Vector<S, N> operator+ (const Vector<S, N>& v) const {
            Vector<S, N> vrtn;
            auto vi = v.begin();
            auto add_v = [vi](S a) mutable { return a + (*vi++); };
            std::transform (this->begin(), this->end(), vrtn.begin(), add_v);
            return vrtn;
        }

        /*!
         * Vector addition operator
         */
        void operator+= (const Vector<S, N>& v) {
            auto vi = v.begin();
            auto add_v = [vi](S a) mutable { return a + (*vi++); };
            std::transform (this->begin(), this->end(), this->begin(), add_v);
        }

        /*!
         * Vector subtraction operator
         */
        Vector<S, N> operator- (const Vector<S, N>& v) const {
            Vector<S, N> vrtn;
            auto vi = v.begin();
            auto subtract_v = [vi](S a) mutable { return a - (*vi++); };
            std::transform (this->begin(), this->end(), vrtn.begin(), subtract_v);
            return vrtn;
        }

        /*!
         * Vector subtraction operator
         */
        void operator-= (const Vector<S, N>& v) {
            auto vi = v.begin();
            auto subtract_v = [vi](S a) mutable { return a - (*vi++); };
            std::transform (this->begin(), this->end(), this->begin(), subtract_v);
        }

        /*!
         * Scalar addition
         */
        template <typename _S=S, std::enable_if_t<std::is_scalar<std::decay_t<_S>>::value, int> = 0 >
        Vector<_S, N> operator+ (const _S& s) const {
            Vector<_S, N> rtn;
            auto add_s = [s](_S coord) { return coord + s; };
            std::transform (this->begin(), this->end(), rtn.begin(), add_s);
            return rtn;
        }

        /*!
         * Scalar addition
         */
        template <typename _S=S, std::enable_if_t<std::is_scalar<std::decay_t<_S>>::value, int> = 0 >
        void operator+= (const _S& s) {
            auto add_s = [s](_S coord) { return coord + s; };
            std::transform (this->begin(), this->end(), this->begin(), add_s);
        }

        /*!
         * Scalar subtraction
         */
        template <typename _S=S, std::enable_if_t<std::is_scalar<std::decay_t<_S>>::value, int> = 0 >
        Vector<_S, N> operator- (const _S& s) const {
            Vector<_S, N> rtn;
            auto subtract_s = [s](_S coord) { return coord - s; };
            std::transform (this->begin(), this->end(), rtn.begin(), subtract_s);
            return rtn;
        }

        /*!
         * Scalar subtraction
         */
        template <typename _S=S, std::enable_if_t<std::is_scalar<std::decay_t<_S>>::value, int> = 0 >
        void operator-= (const _S& s) {
            auto subtract_s = [s](_S coord) { return coord - s; };
            std::transform (this->begin(), this->end(), this->begin(), subtract_s);
        }

        /*!
         * Overload the stream output operator
         */
        friend std::ostream& operator<< <S, N> (std::ostream& os, const Vector<S, N>& v);
    };

    template <typename S=float, size_t N=3>
    std::ostream& operator<< (std::ostream& os, const Vector<S, N>& v)
    {
        os << v.str();
        return os;
    }

    /*
     * Random number generator classes outside scope of code review, but included so
     * that file compiles:
     */

    //! Generate uniform random numbers in a floating point format.
    template <typename T = double>
    class RandUniformReal
    {
    private:
        std::random_device rd{};
        std::mt19937_64 generator{rd()};
        std::uniform_real_distribution<T> dist;
    public:
        //! Default constructor gives RN generator which works in range [0,1)
        RandUniformReal (void) {
            typename std::uniform_real_distribution<T>::param_type prms (T{0}, T{1});
            this->dist.param (prms);
        }
        //! This constructor gives RN generator which works in range [a,b)
        RandUniformReal (T a, T b) {
            typename std::uniform_real_distribution<T>::param_type prms (a, b);
            this->dist.param (prms);
        }
        //! Get 1 random number from the generator
        T get (void) { return this->dist (this->generator); }
        //! Get n random numbers from the generator
        std::vector<T> get (size_t n) {
            std::vector<T> rtn (n, T{0});
            for (size_t i = 0; i < n; ++i) {
                rtn[i] = this->dist (this->generator);
            }
            return rtn;
        }
        T min (void) { return this->dist.min(); }
        T max (void) { return this->dist.max(); }
    };

    //! Generate uniform random numbers in a integer format
    template <typename T = unsigned int>
    class RandUniformInt
    {
    private:
        std::random_device rd{};
        std::mt19937_64 generator{rd()};
        std::uniform_int_distribution<T> dist;
    public:
        //! Default constructor gives an integer RNG which works in range [0,(type max))
        RandUniformInt (void) {
            typename std::uniform_int_distribution<T>::param_type prms (std::numeric_limits<T>::min(),
                                                                        std::numeric_limits<T>::max());
            this->dist.param (prms);
        }
        //! This constructor gives RN generator which works in range [a,b)
        RandUniformInt (T a, T b) {
            typename std::uniform_int_distribution<T>::param_type prms (a, b);
            this->dist.param (prms);
        }
        //! Get 1 random number from the generator
        T get (void) { return this->dist (this->generator); }
        //! Get n random numbers from the generator
        std::vector<T> get (size_t n) {
            std::vector<T> rtn (n, T{0});
            for (size_t i = 0; i < n; ++i) {
                rtn[i] = this->dist (this->generator);
            }
            return rtn;
        }
        T min (void) { return this->dist.min(); }
        T max (void) { return this->dist.max(); }
    };

} // namespace morph

The test program:

#include "Vector.h"
using morph::Vector;
using std::cout;
using std::endl;
using std::array;

int main() {
    int rtn = 0;
    Vector<float, 4> v = {1,2,3};
    cout << "x component of v: " << v.x() << endl;
    v.renormalize();
    cout << "After renormalize: " << v << "; unit vector? " << (v.checkunit() ? "yes" : "no") << endl;
    v.randomize();
    cout << "After randomize: " << v << endl;
    Vector<int, 5> vi;
    vi.randomize();
    cout << "After randomize of 5D int vector: " << vi << endl;
    cout << "Length: " << vi.length() << endl;
    // Test assignment
    Vector<int, 5> vi2 = vi;
    cout << "Copy of int vector: " << vi2 << endl;
    // Test comparison
    cout << "vi == vi2? " << (vi == vi2 ? "yes" : "no") << endl;
    // Test negate
    Vector<int, 5> vi3 = -vi;
    cout << "-ve Copy of int vector: " << vi3 << endl;
    // Test comparison
    cout << "vi == vi3? " << (vi == vi3 ? "yes" : "no") << endl;
    // Test cross product (3D only
    Vector<double, 3> a = {1,0,0};
    Vector<double, 3> b = {0,1,0};
    Vector<double, 3> c = a * b;
    cout << a << "*" << b << "=" << c << endl;
    // dot product
    Vector<int, 2> vv1 = {5,8};
    Vector<int, 2> vv2 = {2,5};
    int dp = vv1.dot (vv2);
    cout << vv1 << "." << vv2 << " = " << dp << endl;
#if 0 // No good:
    // Test init from array
    array<float, 3> arr = { 2,3,4 };
    Vector<float, 3> varr = arr; // Tried overloading operator= to no avail.
    cout << "Vector from array: " << varr << endl;
#endif
    // Scalar multiply
    cout << vv2;
    vv2 *= 2UL;
    cout << " * 2 = " << vv2 << endl;
    Vector<int, 2> vv4 = vv1 * 3;
    cout << "vv1 * 3 = " << vv4 << endl;
#if 0 // No good, as expected:
    // Scalar multiply with different type
    double dbl = 3.0;
    Vector<int, 2> vv5 = vv1 * dbl;
    cout << "vv1 * 3.0 = " << vv4 << endl;
#endif
    // Scalar division
    Vector<double, 3> d = a/3.0;
    cout << "a / 3.0 = " << d << endl;
    // Vector addition
    Vector<double, 3> e = a+b;
    cout << "a + b = " << e << endl;
    // Vector subtraction
    Vector<double, 3> f = a-b;
    cout << "a - b = " << f << endl;
    // Test default template args
    Vector<double> vd_def;
    vd_def.randomize();
    cout << "Vector<double> and randomize: " << vd_def << endl;
    Vector<> v_def;
    v_def.randomize();
    cout << "Vector<> and randomize: " << v_def << endl;

    // dot product of large vector
    const size_t n = 1000000; // Approx limit on my laptop: 1045000
    Vector<float, n> big1;
    Vector<float, n> big2;
    big1.randomize();
    big2.randomize();
    float bdp = big1.dot(big2);
    cout << "(N="<<n<<") big1.big2 = " << bdp << endl;

    // Test set_from
    Vector<double, 3> d1;
    array<double, 3> a1 = { 5,6,7 };
    d1.set_from (a1);
    cout << "After set_from(), d1 should be (5,6,7): " << d1 << endl;
    array<double, 4> a2 = { 5,6,8,8 };
    d1.set_from (a2);
    cout << "After set_from(), d1 should be (5,6,8): " << d1 << endl;

    return rtn;
}
```

Answer 3

I think, I can comprehend your main idea, "to create a specialized class for mathematical N-dimensional vectors" but there are few problems by doing it using inheritance in your case. You should use a truly specialized class for that, think of it, std::vector is a stl container, it is supposed to handle data as a data structure, any class which extends from it, should handle data.

From a design perspective it would be better to create a class which either adapts std::vector as a member (is composed by) or one which has a vector of objects (commonly numerical types, but in extension mathematical objects too).

/*using google c++ style guide (only for naming purposes)*/
template <typename MathObject>
class Vector
{
private:
    //allows to handle a runtime determined size vector
    //could be implemented using std::unique_ptr
    //(define a default constructor in the case that MathObject be a class)
    MathObject *vector_;

    size_t dimension_; //allows you to check the dimension of a vector
public:
    Vector(MathObject* const& vector, size_t dimension) :
        vector_(vector), dimension_(dimension)
    {
    }

    ~Vector()
    {
        delete[] vector_;
    }
    //...
};

Many thanks to Konrad Rudolph who has mentioned the mistakes on your implementation.