# Custom mathematical vector class

My first bigger C++ project: a vector class for personal use and statistical computation.

#include <iostream>
#include <vector>
#include <algorithm>
#include <cmath>

template<class T>
class vect
{
std::vector<T> m;
size_t s;

public:
// default constructor
vect(): m(0), s(0){}

vect(size_t n) :m(n), s(n) {}

vect(std::vector<T> v) :m(v), s(v.size()) {}

// copy constructor
vect(const vect<T>&v): m(v.getData()), s(v.size()){}

// destructor
~vect(){}

std::vector<T> getData() const{
return m;
}

size_t size() const {
return s;
}

m.push_back(value);
s++;
}

s++;
m.emplace(m.begin()+loc, value);
}

void rmFrom(){
m.pop_back();
s--;
}

void rmAt(size_t loc){
m.erase(m.begin()+loc);
s--;
}

// returns a sorted copy of the original vector.
vect<T> sorted(){
std::vector<T> temp = this->getData();
std::sort(temp.begin(), temp.end());
vect<T> v(temp);
return v;
}

double median(){
// linearly interpolated
if (s%2==0)
{
return ((this->sorted()[(s/2)-1]+this->sorted()[(s/2)])/2);
}
else
{
return (this->sorted()[(s)/2]);
}
}

double percentile(double p){
// @p - value between 0 and 1: the percentile
// rounds to the nearest index and return the corresponding
return (this->sorted()[round(s*p)]);
}

double sum(){
// computes the sum of the vectors elements
double sum = 0;
size_t l = 0;
while (l<s) {
sum+=m[l];
l++;
}
return sum;
}

double mean(){
// computes the mean of the vectors elements
return this->sum()/s;
}

double dot(){
// computes the inner product of this  vector
vect<T> temp = *this;
double dot = 0;
size_t l = 0;
while(l<s){
dot+=pow(temp[l],2);
l++;
}
return dot;
}

double dot(vect<T> other){
// computes the inner product of this and another vector
vect<T> temp = *this;
double dot = 0;
size_t l = 0;
while(l<s){
dot+=temp[l]*other[l];
l++;
}
return dot;
}

double magnitude(){
return sqrt(this->dot());
}

double manhattan_norm(){
// computes the manhattan norm
double n = 0;
size_t l = 0;
while(l<s){
n+=abs(m[l]);
l++;
}
return sqrt(n);
}

double p_norm(unsigned int p){
return pow(this->dot(),1/p);
}

vect<T> normalized(){
// devides the vector by it's eucledian distance
// which results in a unity vector, meaning
// the sum the elements of a unity vector add up to one;
T n = 0;
vect<T> unity = *this;
n = unity.magnitude();
return unity/n;
}

vect<T> diff(vect<T> other){
vect<T> diff;
return (*this-other);
}

double cosine(vect<T> other){
// returns the cosine of theta of this vector and another
return this->dot(other)/(this->magnitude()*other.magnitude());
}

double angle(vect<T> other){
return acos(this->cosine(other))*(180/3.14159265);
}

bool is_perpendicular(vect<T> other){
return (this->dot(other)==0);
}

vect<T> direction_cosine(){
vect<T> direction;
double norm = this->norm();
double cosine = 0;
double theta = 0;
for(size_t i = 0; i < s; i++){
cosine = m[i]/norm;
theta = acos(cosine);
}
return direction;
}

vect<T> parallel_comp(vect<T> other){
// projection of other onto *this
vect<T> temp = *this;
return (temp.dot(other)/temp.dot(temp))*temp;
}

vect<T> perpendicular_comp(vect<T> other){
return other-parallel_comp(other);
}

double parallel_magnitude(vect<T> other){
// returns the the projections size
return this->parallel_comp(other).magnitude();
}

double error_magnitude(vect<T> other){
// returns the the parallel components size relative to it's own
return this->perpendicular_comp(other).magnitude();
}

vect<T>& operator=(const vect<T>& other){
if (this!=&other)
{
m = other.getData();
s = other.size();
}
return *this;
}
T& operator[] (size_t i) {
return m[i];
}
const T& operator[] (size_t i) const {
return m[i];
}

};

template<class T>
std::ostream& operator<<(std::ostream &os, const vect<T>&v){
for(int i = 0; i < v.size(); i++){
os << v.getData()[i] << " ";
}
return os;
}

template<class T>
std::istream& operator>>(std::istream& is, vect<T>& v){
T value;
is >> value;
return is;
}

template<class T>
vect<T>& operator+=(vect<T>& a, const vect<T>& b)
{
for(size_t i=0; i<a.size(); ++i)
a[i]+=b[i];
return a;
}

template<class T>
vect<T> operator+(const vect<T>& a, const vect<T>& b)
{
vect<T> z(a.size());
for(size_t i=0; i<a.size(); ++i)
z[i] = a[i]+b[i];
return z;
}


Plus many more overloaded arithmetic operators which are all like the last two ones.

I know it's nothing fancy but I am always looking for ways to make things better.

• What you can't do after receiving answers Sep 29, 2014 at 15:49
• What is the purpose of the s member? It looks like it's just the same as m.size(). Jul 16, 2021 at 15:24

### C++14

Its 2014 most modern compilers now support C++14 so you should use it. This code is still very C++03. For this class this simply means adding move semantics (and nothrow on swap).

To add move semantics you need to add a move constructor and move assignment operator.

vect(vect&& other);
vect& operator=(vect&& other);


Potentially there are a couple of places where you can use range based for operator.

### General

This seems redundant (and thus dangerous).

size_t s;


The size is already stored as part of the the other member m;

### Const Correctness

None of these functions

double sum();
double mean();
double dot();


modify the state of your vector. So you should mark them as const members.

double sum()  const;
double mean() const;
double dot()  const;


This allows you pass your vect to a function as a const reference and still call these non mutating functions.

### Efficiency

You may want to cache the sum in a mutable member. Invalidate it if the vector is mutated. There is no point re-computing this value if the vector has not changed.

double sum(){
// computes the sum of the vectors elements
double sum = 0;
size_t l = 0;
while (l<s) {
sum+=m[l];
l++;
}
return sum;
}


Also prefer to use some of the algorithms that are built for you when you can.

 sum = std::accumulate(m.begin(), m.end());


### Copy and Swap

This is an old way of doing this:

// overloaded operators
vect<T>& operator=(const vect<T>& other){
if (this!=&other)
{
m = other.getData();
s = other.size();
}
return *this;
}


The more modern way is to use the copy and swap idium.

// overloaded operators
vect<T>& operator=(vect<T> other){  // Notice the pass by value to get the copy.

other.swap(*this);    // Swap the content of this object.
// With the copy you just created. This updates it
// in an exception safe way,
return *this;         // return yourself.
}                        // Let the destructor of other cleanup your old state.

void swap(vect<T>& other) throws() { // use C++11 nothrow if you upgrade to C++11
std::swap(this->m, other.m);
std::swap(this->s, other.s);
}


### The X operators can be written in terms of the X= operator.

Example:
Define the + operator in terms of the += operator.

template<class T>
vect<T> operator+(const vect<T>& a, const vect<T>& b)
{
vect<T>  result(a);  // Make a copy.
result += b;
return result;
}


Once you studied that and see it does the same. There is a small optimization.

template<class T>
vect<T> operator+(const vect<T> copy, const vect<T>& b)
{
copy += b;       // Notice the copy is passed by value.
// So there is an implicit copy
// So you don't need the manual copy inside the code.
return copy;
}


### Your streaming operations should be symmetrical and distinguishable.

The output operator dumps the whole vector:

template<class T>
std::ostream& operator<<(std::ostream &os, const vect<T>&v){
for(int i = 0; i < v.size(); i++){
os << v.getData()[i] << " ";
}
return os;
}


But the read operator only reads a single value into the vector (so its not symmetrical).

template<class T>
std::istream& operator>>(std::istream& is, vect<T>& v){
T value;
is >> value;    // Also note. That if the read fails.
// You should probably test to make sure the read works.
// if (is >> value) {v.addTo(value);}
return is;
}


Personally I would make it dump the data in a way that makes it obvious of the start/end point (or) you can prefix the dump with a count and then the values.

template<class T>
std::ostream& operator<<(std::ostream &os, const vect<T>&v){
os << v.size() << ": ";
for(int i = 0; i < v.size(); i++){
os << v.getData()[i] << " ";
}
return os;
}

// Note: This does not work for T == std::string as it does not
//       have symmetric input output operators. But all other normal types do.
template<class T>
std::istream& operator>>(std::istream& is, vect<T>& v){
std::size  s       = 0;
char       marker  = 'B';

if (is >> s >> marker)
{
if (marker != ':')
{   // The mark is not what we expect.
// So mark the stream as bad so that processing stops.
is.setstate(std::ios::failbit);
}
else
{
T  value;
for(;s > 0 && (is >> value);--s)
{
// Have values and successful read.
}
}
}
return is;
}

• Some minor remarks: "range based for operator" looks like a typo (?). "You may want to cache the sum in a mutable member." Might violate Single Responsibility Principle. "Also prefer to use some of the algorithms" Somehow was swallowed by a code block. accumulate requires a third argument. "The more modern way is to use the copy and swap idium." Typo, and don't let Howard Hinnant hear this ;) "There is a small optimization" First parameter must not be const.
– dyp
Sep 27, 2014 at 0:00

### Naming

Your "vector" is not really like std::vector, but vect sounds like it. It would be more intuitive to rename it some thing different, for example vectorstats, or just vstats, or something.

The m and s variables are terrible.

### size

I don't really see the point of the s variable. Why do you want to count the size yourself? Why not use m.size()? That would be a lot less error-prone, because as a general rule, the less lines of code you write yourself the better.

### sum

It would be more natural to use the iterator pattern, shorter and better:

double sum() {
double sum = 0;
for (std::vector<T>::const_iterator it = m.begin(); it != m.end(); ++it) {
sum += *it;
}
return sum;
}


### dot

The temp variable is completely pointless in this method:

double dot(){
// computes the inner product of this  vector
vect<T> temp = *this;
double dot = 0;
size_t l = 0;
while(l<s){
dot+=pow(temp[l],2);
l++;
}
return dot;
}


You could rewrite as:

double dot() {
double dot = 0;
for (std::vector<T>::const_iterator it = m.begin(); it != m.end(); ++it) {
dot += pow(*it, 2);
}
return dot;
}


The same is true for the overloaded version of this method.

### Excessive parentheses

The too many parantheses are seriously hurting readability here:

return ((this->sorted()[(s/2)-1]+this->sorted()[(s/2)])/2);


And the lack of spaces around operators don't help!

This should have been:

return (this->sorted()[s / 2 - 1] + this->sorted()[s / 2]) / 2;


Or actually, it's wasteful to sort the vector twice, better sort once and save in a temp variable:

std::vector<T> temp = this->sorted();
return (temp[s / 2 - 1] + temp[s / 2]) / 2;


Now that I look at this cleaner version, it's clear that this won't work if the vector is empty, because if s == 0 this will end up referencing temp[-1]. This was not so easy to see in the original version, now it's obvious.

This was just the worst example I found, but in many many places you use far more parentheses than you need. I suggest to review the entire code and trim a little bit.

Also try to put spaces around operators like I did in this example.

### Placement of curly braces

You're placing curly braces inconsistently. Sometimes you use put the opening brace on the same line as the statement, like this:

double sum(){
// ...
while (l<s) {
// ...
}
}


Other times you place opening brace on the next line, like this:

if (s%2==0)
{
// ...
}
else
{
// ...
}


I suggest to pick either of these styles and stick to it. (I prefer the first.)

• The class name should be capitalized so that it's easily distinguished from functions and variables.

• m is a very undescriptive name as it's only one letter. Write out the entire name instead of giving a random letter or some abbreviation. It's especially important to do this so that others can know what this vector is storing.

• You don't need s; the size is already known to the vector structure. Just access its size via the size() member function.

• If you're just summing the values in a vector, you can use std::accumulate():

return std::accumulate(m.cbegin(), m.cend(), 0);


(That would be the entire body of sum().)

1. You don't need an empty destructor. Since it is a no-op. You can leave it out and the compiler will provide a default for you.

2. The default constructor, vect() should not init m with 0.

3. The copy constructor, vect(const vect<T>&v) is not needed either. The compiler will generate the exact same code for you if you don't provide it.

4. I don't really see the need for the s variable here. It is just keeping a copy of m.size(). Why not just use std::vector::size() for that? It is a bunch of extra work to keep that var up-to-date.

5. std::vector<T> getData() const might be more efficient if returning by reference. If the intended use for this method is to allow read-only access to the internal vector, then there is no reason to return by copy. Return by const reference: const std::vector<T> & getData() const.

6. rmFrom() is a very unclear name. rm is short for remove? Then just call it remove.

7. The methods that add data to the vector, such as addTo()/addAt() should be provided in two flavors: one taking a const T& and one taking a T&&. If you refer to std::vector:push_back() you will see that it provides these two overloads to optimize for move semantics (move semantics will require a C++11 capable compiler).

8. You should probably also provide a move constructor and a move assignment operator for your vect. See the rule of three/five/zero.

9. Unneeded temps inside dot(). Why are you creating temporary copies of the vector inside both dot() methods? That is nonsense since both function only read from the vector.

10. Your code is not fully const correct. There are several methods that need to be made const.

11. In vect<T> diff(vect<T> other) there is a local variable inside (also named diff) that is never used.

12. You are using the <cmath> functions. So functions like sin, cos, acos, etc should be prefixed with std::. E.g.: std::acos(). That would be the correct form.

13. The array access operators [] would benefit from some runtime bounds checking. I suggest adding asserts to those. E.g.: assert(i < m.size());

Few missing poins.

• STL provides convenient algorithms for pretty much all loops in your code. Besides std::accumulate, take a look at std::inner_product and std::transform algorithms.

• I question a mathematical validity of p_norm. It calculates $(\sum {m_i}^2)^\frac{1}{p}$. Shouldn't it be $(\sum {m_i}^p)^\frac{1}{p}$ instead?

• Certain operations (such as operator+= and dot) only make sense for vectors of the same dimension. Depending on the size of your vectors you'd be getting either a not very meaningful result or an exception. Doesn't seem right. I don't know your use case; I'd recommend to address the problem anyway.

I'd probably implement this as at least two distinct classes: one class to implement statistical functions, and one or more classes to implement linear algebra.

The statistical classes could have additional functions such as sample variance, estimated population variance, and covariances. This might lead to a bit of duplication of code (since some of the operations here look a lot like inner product, in fact are inner product if you take a vector-space approach to your statistics) but I think this might be worth the benefit of making it clear what functions are best suited to whatever domain in which you're currently working. (A particular side benefit is that in a development environment that offers autocompletion of member function names, it would offer the ones that you might actually want to use and not a lot of functions that make no sense at all.)

For a linear-algebra vector, you could make the dimension of the vector space be one of the template parameters. It should still be possible to construct a higher-dimension vector by adjoining a new component to a lower-dimension vector. (If you do a lot of work in two-, three-, or four-dimensional vector spaces, you might want vector templates for each of those specific dimensions in which the iterative functions are completely unrolled, though it's possible that the compiler will do this anyway when the dimension is explicit in the class definition; it would be interesting to benchmark two templates against each other, one designed exclusively for two-dimensional vectors and the other an N-dimensional vector with the template parameter N set to 2.)

I would provide a better approximation of pi (more digits). I probably would define pi as a constant (possibly a static const member variable, although possibly just a static variable in the compilation unit).

Arc cosine does not accurately measure the angle between two vectors when that angle is small. The angle between two vectors is twice the arc sine of half the magnitude of the difference between two unit vectors parallel to the original two vectors; this works well for small angles, reasonably well for other acute angles (though it's probably more computation than the arc cosine method), but is inaccurate for angles near pi, where it would be better to take the sum of the unit vectors rather than their difference. I'd suggest comparing the dot product to the product of the magnitudes of the two vectors at runtime, and using that result to decide whether to use the cosine method or one of the other methods. There are some tweaks you can do to try to minimize the number of square roots you have to compute (such as, initially compare the squares of two quantities rather than first computing the quantities themselves).

The percentile function looks dangerous. What happens if you ask for the 99th percentile of a vector with only forty elements? Since 0.99*40 = 39.6, which rounds to 40, you would be trying to access m[40], which is beyond the end of your data set. I think you need to decide what you want percentiles to do, and throw an exception when the function call cannot produce a correct result. It seems to me you might want d.percentile(0.5) to do exactly the same thing as d.median(), which is not true in your current implementation when d has an even number of elements.

• Thank you for your comments. I did not split the vector in two different ones, rather took the methods (linalg) and (stats) into two different namespaces outside of the class. I agree that the class was ambiguous in it's methods. Oct 7, 2014 at 6:07

Class Design Changes Most notably is that I decided follow the advice of @glampert and decided to go with the rule of zero, and drop any other destructors/copy/move constructors. I played around with different version and came to the conclusion that at this points no additional constructors are required. So I dropped them.

Use of STD Library Another recurring topic was to make use of the and libraries. This was a great advice. I incorporated them quite a view times.

Other changes Changes of class and function names. Adding to additional constructor: initialiser list and "parameter pack". Dropping the unnecessary size_t in the private data section. Dropping some of the functions.

 #include <iostream>          // std::cout
#include <initializer_list>  // allows to write data onto in this form: MVector v = {...}
#include <vector>            // container for private data of the class
#include <algorithm>         // swap...
#include <numeric>           // accumulate, innerproduct
#include <functional>        // function as parameter
#include <random>            // random numer generatore and distributions
#include <cmath>             // trigon functions
#include <assert.h>

// ************ Description (short) ***************
// Custom vector class for personal usage.

// ************** Definition Vector ***************
// Definition http://en.wikipedia.org/wiki/Euclidean_vector:
// In mathematics, physics, and engineering, a Euclidean vector
// (sometimes called a geometric or spatial vector or—as here—simply a vector)
// is a geometric object that has magnitude (or length) and direction
// and can be added to other vectors according to vector algebra.

// *************** Purpose of Class ***************
// Primary purpose: of this class is primarily to store numeric data.
// Secondary purpose: simple construction of a vector object (ease of use)
// Quaternary purpose: and inclusion of vector operations and algorithms
// commonly used in linear algebra and statistics.

// ******************** Data **********************
// The data is stored in a std::vector, which provides random access to it's elements.

#ifndef MVECTOR_H
#define MVECTOR_H

template<class T>
class MVector
{
std::vector<T> data;

public:

// create a vector of size n either with xvalue or lvalue
// requires an object of std::size_t to construct a vector
// of n zeros.
explicit MVector(const std::size_t& n) :data(n) {}
explicit MVector(std::size_t&& n) :data(n) {}

// construct a vector by initialising an unspecified amount of values
template<class ... Data>
explicit MVector(T first, Data&&... values): data{first,
std::forward<T>(static_cast<T>(values))... }{}

// construct a vector by initializer list
MVector(std::initializer_list<T> il) :data(il){}

// construct with a std::vector by initializer list
explicit MVector(const std::vector<T>& v) :data(v){}
explicit MVector(std::vector<T>&& v) :data(v){}

// returns the data
const std::vector<T>& getData() const{
return data;
}

// returns the amount of elements stored in the vector
std::size_t size() const{
return data.size();
}

// add an elments at the end of the vector
// @ param "value" - object of class T
data.push_back(value);
}

// add an elments at the end of the vector
// @ param "value" - object of class T
data.push_back(value);
}

// add an elments at specified location in the vector
// @ param "value" - object of class T
// @ param "loc" - location of the input
void addAt(const T& value, const std::size_t& loc=0){
data.emplace(data.begin()+loc, value);
}

// add an elments at specified location in the vector
// @ param "value" - object of class T
// @ param "loc" - location of the input
data.emplace(data.begin()+loc, value);
}

// removes an elments from the end of the vector
void remove_From(){
data.pop_back();
}

// removes an elments from the specified location of the vector
// @ param "loc" - index which element to delete
void remove_At(const std::size_t& loc){
data.erase(data.begin()+loc);
}

// removes an elments from the specified location of the vector
// @ param "loc" - index which element to delete
void remove_At(std::size_t&& loc){
data.erase(data.begin()+loc);
}

// resizes the vector to the size of to
// @ param "to" - size of which the vector will be changed to
void resize(std::size_t to){
data.resize(to);
}

// returns a sorted COPY of the original vector.
MVector<T> sorted(){
std::vector<T> temp = this->getData();
std::sort(temp.begin(), temp.end());
MVector<T> v(temp);
return v;
}

// convenience functions

// returns the median, linearly interpolated, of the vector
double median() {
MVector<T> temp = this->sorted();
if (data.size()%2==0){
return (temp[data.size() / 2 - 1] + temp[data.size() / 2]) / 2;
}
else{
return temp[data.size() / 2];
}
}

// returns the p'th percentile of the vector
// rounds to the nearest index and return the corresponding
// @ param "p" - value between 0 and 1: the percentile
double percentile(double p) {
assert(p>=0 && p<=1);
return (this->sorted()[round(data.size()*p)]);
}

// returs the sum of all objects in the vector
// @ param "loc" - adds a constant value to all parameters
double sum(T loc=0) const{
return std::accumulate(data.begin(), data.end(), loc);
}

// returs the arithmetic mean of the vector
double mean() const{
return this->sum()/data.size();
}

// returs the inner product of the vector
// @ param "loc" - adds a constant value to all parameters
double dot(T loc=0) const{
return std::inner_product(data.begin(), data.end(),
data.begin(), loc);
}

// returs the inner product of two different vectors
// @ param "loc" - adds a constant value to all parameters
double dot(MVector<T> other, T loc=0) const{
return std::inner_product(data.begin(), data.end(),
other.getData().begin(), loc);
}

// returns the "norm" or "length" of a vector
double magnitude() const{
return sqrt(this->dot());
}

// returns the manhatten "norm" or "length" of a vector
double manhattan_norm(double loc=0) const{
return std::sqrt(
std::accumulate(data.begin(), data.end(), loc,
[](double x, double y){
return std::abs(x)+std::abs(y);}));
}

// devides the vector by it's eucledian distance
// which results in a unity vector, meaning
// the sum the elements of a unity vector add up to one;
MVector<T> normalized() const{
double n = 0;
MVector<T> unity = *this;
n = unity.magnitude();
return unity/n;
}

// returns the cosine of theta of this vector and another
// @ param "other" - other vector
double cosine(MVector<T> &other) const{
return this->dot(other)/(this->magnitude()*other.magnitude());
}

// returns the angel between to vectors (not exact, pi is approximated)
// @ param "other" - other vector
double angle(MVector<T> &other) const{
return std::acos(this->cosine(other))*(180/3.14159265);
}

// returns true if the angel between to vectors
// is 90 degrees, which is exactly the case when the dot product is zero.
// @ param "other" - other vector
bool is_perpendicular(const MVector<T> &other) const{
return (this->dot(other)==0);
}

// projection of other onto *this
// requires copy of *this because multiplying *this resulted in compilation error
// @ param "other" - other vector
MVector<T> parallel_comp(const MVector<T>& other) const{
MVector<T> temp = *this;
return (temp.dot(other)/temp.dot(temp))*temp;
}

// return the "error component" of two vectors
// @ param "other" - other vector
MVector<T> perpendicular_comp(const MVector<T>& other) const{
return other-parallel_comp(other);
}

T& operator[] (std::size_t i) {
assert(i<data.size());
// make sure that the index is not out of bound
return data[i];
}

const T& operator[] (std::size_t i) const {
assert(i<data.size());
// make sure that the index is not out of bound
return data[i];
}

MVector<T>& operator=(MVector<T> other){
this->data.swap(other.data);
return *this;
}

// end of class
};

// *****************************************************************************
// bool  operators
// *****************************************************************************

template<class T>
bool operator==(const MVector<T> &lhs,
const MVector<T> &rhs) {
return std::equal(lhs.getData().begin(), lhs.getData().end(),
rhs.getData().begin());
}

template<class T>
bool operator!=(const MVector<T> &lhs,
const MVector<T> &rhs) {
return (! (lhs == rhs) );
}

// *****************************************************************************
// input output operators
// *****************************************************************************

template<class T>
std::ostream& operator<<(std::ostream &os, const MVector<T>& v){
for(int i = 0; i < v.size(); i++){
os << v.getData()[i] << " ";
}
return os;
}

template<class T>
std::istream& operator>>(std::istream& is, MVector<T>& v){
T value;
is >> value;
return is;
}

// *****************************************************************************

// // vector arithmetic overloaded operators

// *****************************************************************************

// *****************************************************************************
// *****************************************************************************

template<class T>
MVector<T>& operator+=(MVector<T>& lhs, const MVector<T>& rhs){
assert(lhs.size() == rhs.size());
for(std::size_t i=0; i < lhs.size(); ++i){
lhs[i]+=rhs[i];
}
return lhs;
}

template<class T>
MVector<T> operator+(MVector<T> lhs, const MVector<T>& rhs){
assert(lhs.size()==rhs.size());
lhs+=rhs;
return lhs;
}

template<class T>
MVector<T>& operator+=(MVector<T>& lhs, const T& rhs){
for(std::size_t i=0; i<lhs.size(); ++i)
lhs[i]+=rhs;
return lhs;
}

template<class T>
MVector<T> operator+(MVector<T> lhs, const T& rhs){
lhs+=rhs;
return lhs;
}

template<class T>
MVector<T> operator+(const T& lhs, MVector<T> rhs){
MVector<T> result(rhs.size());
for(std::size_t i=0; i<rhs.size(); ++i)
result[i] = lhs+rhs[i];
return result;
}

// *****************************************************************************
// // vector creation routines: constant/random space
// *****************************************************************************

template <class T>
MVector<T> range(T to, T from = 0, T by = 0){
MVector<T> v;
if (from>to){
for (T i=from; i > to; i-=by) {
}
}
else{
for (T i=from ; i < to; i+=by) {
}
}
return v;
}

template<class T>
MVector<T> range_function(std::size_t to, std::function<T(T)> fun, T from=0, T by=1){
MVector<T> v;
if (from<to){
for (T i = from; i < to; i+=(by)) {
}
}
else{
for (T i = from; i > to; i-=(by)) {
}
}
return v;
}

template <class T>
MVector<T> rgeom(std::size_t length, T probability){
assert(probability>=0 && probability<=1);
std::mt19937 gen;
gen.seed(time(NULL));
std::geometric_distribution<T> geometric(probability);
MVector<T> v;
for (std::size_t i = 0; i < length; i++) {