#
Implementing Vector types in C++

Vector types are used in almost any kind of game code. Position, velocity, acceleration, forces, points in polygons and polyhedra etc. are all using some sort of vector type. Writing a vector class isn't particularly hard, but since it's such a central part of the code it might be worth it to have a closer look. Over the last couple years I implemented and customized several versions of vector classes and will now try to present my collected findings and thoughts here.

[subheading]Two types of Vectors[/subheading]

While mathematically the same, "small" and "big" vectors are not identical from a programmers perspective. Especially in the case of game programming we typically have the small vectors that are two or three dimensional and are used for positions, velocities and so forth. Then there are the vectors we use for example when solving linear equation systems. These are often way larger and their dimension isn't necessarily known beforehand. Simply put, there are the vectors that scale in dimension with the problem size and those that don't. For example the vectors you need to represent positions of units in a game have fixed size. No matter how many units there are in the game, their positions always are two- or three-dimensional. On the other hand if you are calculating splines (which requires you to solve a linear equation system) the vector size in your linear system depends on the amount of points you want to calculate the spline for. This concept of knowing the size beforehand is relevant to programming C++ since you have to decide whether to have the size defined statically at compile time or dynamically at run time. This is also a performance concern since dynamic sizes always require memory allocations while static sizes allow you to place the vector data directly on the call stack and compilers usually have an easier time optimizing for the static case.

In this entry I'll focus on the static/small vector case.

[subheading]So how should my class look like?[/subheading]

Coming from mathematics a lot of people like calling the components of their vectors x, y and z. This leads to the obvious implementation of a vector struct as:

[source lang="cpp"]

struct Vector3d {

Vector3d() : x(0), y(0), z(0) {};

Vector3d(double xx, double yy, double zz) : x(xx), y(yy), z(zz) {};

double x, y, z;

};

[/source]

Then we implement a bunch of overloaded operators and functions like

[source lang="cpp"]

inline Vector3d operator+(const Vector3d &a, const Vector3d &b)

{

return Vector3d(a.x+b.x, a.y+b.y, a.z+b.z);

}

inline double dot(const Vector3d &a, const Vector3d &b)

{

return a.x*b.x + a.y*b.y + a.z*b.z;

}

//...

[/source]

But in C++ there isn't really a reason why we should restrict ourselves to a type and dimension at this point. So we probably should make type and dimension template arguments. For the type this is trivial, we just replace every "double" with "T". But the dimension is a problem with the way we named the components. A obvious solution to allow arbitrary dimensions is to have the components in an array.

[source lang="cpp"]

template

class Vector {

public:

Vector3d() { std::fill(data, data+D, T()); }

inline T& operator[](unsigned i) { return data; }

inline const T& operator[](unsigned i) const { return data; }

private:

T data[D];

};

[/source]

Sadly we lose the ability to access the components via x, y and z. If we really want that we can explicitly specialize the template and add references. For the 3D case this would look like:

[source lang="cpp"]

template

class Vector {

public:

Vector3d() : x(data[0]), y(data[1]), z(data[2]) { std::fill(data, data+D, T()); }

inline T& operator[](unsigned i) { return data; }

inline const T& operator[](unsigned i) const { return data; }

T &x, &y, &z;

private:

T data[D];

};

[/source]

I personally don't like doing this since writing myvec[0] instead of myvec.x is just as intuitive to me and it automatically reduces the chances that I'm writing copy/paste code since having the index operator in the sources makes it more obvious where I should use loops. Also as soon as you write a function using x, y and z, chances are you're restricting yourself to writing that function for the 2D or 3D case when you could just as well write it for the general case.

The operators and functions obviously have to be rewritten:

[source lang="cpp"]

template

inline Vector operator+(const Vector &a, const Vector &b)

{

Vector result;

for(unsigned i = 0;i {

result = a + b;

}

return result;

}

template

inline T dot(const Vector &a, const Vector &b)

{

T result = 0;

for(unsigned i = 0;i {

result += a * b;

}

return result;

}

//...

[/source]

Woha... suddenly we have loops and local variables? Isn't that going to affect performance?

Luckily compilers optimize this kind of stuff. Since D is known at compile time and usually small, the compiler will unroll that loop and after inlineing it will probably also optimize away the local variables if possible

[subheading]Getting rid of temporary objects[/subheading]

There is still room for improvement in the way vector operations are handled. The operators produce temporary objects for intermediate results and that can be expensive when the temporary object is of nontrivial size.

The compiler will very likely turn this

[source lang="cpp"]

Vector a,b,c,d;

a = b+c+d;

[/source]

into something like

[source lang="cpp"]

Vector tmp1;

for(unsigned i = 0;i tmp1 = b+c;

Vector tmp2;

for(unsigned i = 0;i tmp2 = tmp1+d;

for(unsigned i = 0;i a = tmp2;

[/source]

when the whole operation could actually be written as

[source lang="cpp"]

for(unsigned i = 0;i a = b + c + d;

[/source]

which is obviously preferable.

The solution to this problem are expression templates. The idea is to have operators and functions return a operation object instead of carrying out the operation themselves. So an expression like the one above will be turned into a tree of operation objects which can then be evaluated by the assignment operator.

To achieve this we first need a base class for everything that has "vectorlike behavior":

[source lang="cpp"]

template

class VectorExpr {

public:

inline operator const A&() const //cast operator for convenience

{

return *static_cast(this);

}

};

[/source]

The new vector class now inherits from this and also gets a modified operator= and "copy constructors"

[source lang="cpp"]

template

class Vector : public VectorExpr > {

public:

static const unsigned Dim = D;

typedef T Type;

Vector() { std::fill(data, data+Dim, T()); }

template

Vector(const VectorExpr& a)

{

const A& ao ( a );

for(unsigned i = 0;i data = ao;

}

template

Vector& operator=(const VectorExpr& a)

{

const A& ao ( a );

for(unsigned i = 0;i data = ao;

return *this;

}

T& operator[] (unsigned i) { return data; }

const T& operator[] (unsigned i) const { return data; }

private:

T data[Dim];

};

[/source]

And finally the operators look like

[source lang="cpp"]

template

class AddOp : public VectorExpr > {

public:

AddOp(const A& pa, const B& pb) : a(pa), b(pb) { }

inline T operator[](unsigned i) const { return a + b; }

private:

const A& a;

const B& b;

};

template

inline AddOp

operator+(const VectorExpr &a, const VectorExpr &b)

{

return AddOp(a, b);

}

[/source]

The whole AddOp object will get optimized away at compile time and the result is the nice loop we wanted without any temporary objects (apart from the operation objects that aren't actually constructed at runtime).

And here is the full source code of the vector class I'm using at the moment:

[source lang="cpp"]

/*

* MathVector.h - Copyright (c) 2011 Jakob Progsch

*

* This software is provided 'as-is', without any express or implied

* warranty. In no event will the authors be held liable for any damages

* arising from the use of this software.

*

* Permission is granted to anyone to use this software for any purpose,

* including commercial applications, and to alter it and redistribute it

* freely, subject to the following restrictions:

*

* 1. The origin of this software must not be misrepresented; you must not

* claim that you wrote the original software. If you use this software

* in a product, an acknowledgment in the product documentation would be

* appreciated but is not required.

*

* 2. Altered source versions must be plainly marked as such, and must not be

* misrepresented as being the original software.

*

* 3. This notice may not be removed or altered from any source

* distribution.

*/

/*

* MathVector.h provides a simple static vector template class with

* a basic expression template ansatz for vector operations.

*/

#ifndef VECTOR_H

#define VECTOR_H

#include

#include

#include

#include

template class Vector;

//base class for all expression templates

template

class VectorExpr {

public:

inline operator const A&() const

{

return *static_cast(this);

}

};

//better use macros instead of copy pasting this stuff all over the place

#define MAKE_VEC_VEC_EXPRESSION(NAME, EXPR, FUNCTION) \

template \

class NAME : public VectorExpr > { \

public: \

NAME(const A& pa, const B& pb) : a(pa), b(pb) { } \

inline T operator[](unsigned i) const { return EXPR; } \

private: \

const A& a; \

const B& b; \

}; \

\

template \

inline NAME \

FUNCTION(const VectorExpr &a, const VectorExpr &b) \

{ \

return NAME(a, b); \

}

#define MAKE_VEC_SCAL_EXPRESSION(NAME, EXPR, FUNCTION) \

template \

class NAME : public VectorExpr > { \

public: \

NAME(const A& pa, const T& pb) : a(pa), b(pb) { } \

inline T operator[](unsigned i) const { return EXPR; } \

private: \

const A& a; \

const T& b; \

}; \

\

template \

inline NAME \

FUNCTION(const VectorExpr &a, const T &b) \

{ \

return NAME(a, b); \

}

#define MAKE_SCAL_VEC_EXPRESSION(NAME, EXPR, FUNCTION) \

template \

class NAME : public VectorExpr > { \

public: \

NAME(const T& pa, const A& pb) : a(pa), b(pb) { } \

inline T operator[](unsigned i) const { return EXPR; } \

private: \

const T& a; \

const A& b; \

}; \

\

template \

inline NAME \

FUNCTION(const T &a, const VectorExpr &b) \

{ \

return NAME(a, b); \

}

#define MAKE_VEC_EXPRESSION(NAME, EXPR, FUNCTION) \

template \

class NAME : public VectorExpr > { \

public: \

NAME(const A& pa) : a(pa) { } \

inline T operator[](unsigned i) const { return EXPR; } \

private: \

const A& a; \

}; \

\

template \

inline NAME \

FUNCTION(const VectorExpr &a) \

{ \

return NAME(a); \

}

//create actual functions and operators

MAKE_VEC_VEC_EXPRESSION (EMulExpr, a * b, multiply_elements)

MAKE_VEC_VEC_EXPRESSION (EDivExpr, a / b, divide_elements)

MAKE_VEC_VEC_EXPRESSION (AddExpr, a + b, operator+)

MAKE_VEC_VEC_EXPRESSION (SubExpr, a - b, operator-)

MAKE_VEC_SCAL_EXPRESSION(DivExpr, a / b, operator/)

MAKE_VEC_SCAL_EXPRESSION(MulExpr1, a * b, operator*)

MAKE_SCAL_VEC_EXPRESSION(MulExpr2, a * b, operator*)

MAKE_VEC_EXPRESSION (NegExpr, -a, operator-)

MAKE_VEC_VEC_EXPRESSION (EMaxExpr, std::max(a,b), max)

MAKE_VEC_VEC_EXPRESSION (EMinExpr, std::min(a,b), min)

//sub vector proxy

template

class SubVectorExpr : public VectorExpr > {

public:

SubVectorExpr(const A& pa, const unsigned& o)

: a(pa), offset(o) { }

inline T operator[](unsigned i) const

{

return a[i+offset];

}

private:

const A& a;

const unsigned offset;

};

template

inline SubVectorExpr

SubVector(const VectorExpr &a, const unsigned &o)

{

return SubVectorExpr(a, o);

}

//actual vector class

template

class Vector : public VectorExpr > {

public:

static const unsigned Dim = D;

typedef T Type;

Vector()

{

std::fill(data, data+Dim, T());

}

Vector(const T &p1)

{

BOOST_STATIC_ASSERT(Dim==1);

data[0] = p1;

}

Vector(const T &p1, const T &p2)

{

BOOST_STATIC_ASSERT(Dim==2);

data[0] = p1;

data[1] = p2;

}

Vector(const T &p1, const T &p2, const T &p3)

{

BOOST_STATIC_ASSERT(Dim==3);

data[0] = p1;

data[1] = p2;

data[2] = p3;

}

Vector(const T &p1, const T &p2, const T &p3, const T &p4)

{

BOOST_STATIC_ASSERT(Dim==4);

data[0] = p1;

data[1] = p2;

data[2] = p3;

data[3] = p4;

}

T* raw() {

return data;

}

template

Vector(const VectorExpr& a)

{

const A& ao ( a );

for(unsigned i = 0; i data = ao;

}

T& operator[] (unsigned i) {

return data;

}

const T& operator[] (unsigned i) const {

return data;

}

const Vector& operator*=(const T &b)

{

for(unsigned i = 0; i data *= b;

return *this;

}

const Vector& operator/=(const T &b)

{

for(unsigned i = 0; i data /= b;

return *this;

}

template

const Vector& operator+=(const VectorExpr& a)

{

const A& ao ( a );

for(unsigned i = 0; i data += ao;

return *this;

}

template

const Vector& operator-=(const VectorExpr& a)

{

const A& ao ( a );

for(unsigned i = 0; i data -= ao;

return *this;

}

template

Vector& operator=(const VectorExpr& a)

{

const A& ao ( a );

for(unsigned i = 0; i data = ao;

return *this;

}

Vector& normalize()

{

*this /= abs(*this);

return *this;

}

private:

T data[Dim];

};

//actual InPlaceVector class

template

class InPlaceVector : public VectorExpr > {

public:

static const unsigned Dim = D;

typedef T Type;

InPlaceVector(T *const d) : data(d)

{ }

T* raw() {

return data;

}

T& operator[] (unsigned i) {

return data;

}

const T& operator[] (unsigned i) const {

return data;

}

const InPlaceVector& operator*=(const T &b)

{

for(unsigned i = 0; i data *= b;

return *this;

}

const InPlaceVector& operator/=(const T &b)

{

for(unsigned i = 0; i data /= b;

return *this;

}

template

const InPlaceVector& operator+=(const VectorExpr& a)

{

const A& ao ( a );

for(unsigned i = 0; i data += ao;

return *this;

}

template

const InPlaceVector& operator-=(const VectorExpr& a)

{

const A& ao ( a );

for(unsigned i = 0; i data -= ao;

return *this;

}

template

InPlaceVector& operator=(const VectorExpr& a)

{

const A& ao ( a );

for(unsigned i = 0; i data = ao;

return *this;

}

InPlaceVector& normalize()

{

*this /= abs(*this);

return *this;

}

private:

T *const data;

};

template

InPlaceVector MakeVector(T *const p)

{

return InPlaceVector(p);

}

template

inline Vector eval(const VectorExpr& a)

{

return Vector(a);

}

//"reduction" functions that don't return expression templates

template

inline T sum(const VectorExpr& a)

{

const A& ao ( a );

T res = 0;

for(unsigned i = 0; i res += ao;

return res;

}

template

inline T dot(const VectorExpr& a, const VectorExpr& b)

{

return sum(multiply_elements(a, b));

}

template

inline T squared_norm(const VectorExpr& a)

{

const A& ao ( a );

T res = 0;

for(unsigned i = 0; i {

T tmp = ao;

res += tmp*tmp;

}

return res;

}

template

inline T abs(const VectorExpr& a)

{

return std::sqrt(squared_norm(a));

}

template

inline T norm(const VectorExpr& a)

{

return std::sqrt(squared_norm(a));

}

template

inline Vector normalize(const VectorExpr& a)

{

return eval(a/norm(a));

}

template

std::ostream& operator<<(std::ostream& out, const VectorExpr& a)

{

const A& ao ( a );

out << '(' << ao[0];

for(unsigned i = 1; i {

out << ", " << ao;

}

out << ')';

return out;

}

#undef MAKE_VEC_VEC_EXPRESSION

#undef MAKE_VEC_SCAL_EXPRESSION

#undef MAKE_SCAL_VEC_EXPRESSION

#undef MAKE_VEC_EXPRESSION

#endif

[/source]

Suggestions and critique appreciated.

[subheading]Two types of Vectors[/subheading]

While mathematically the same, "small" and "big" vectors are not identical from a programmers perspective. Especially in the case of game programming we typically have the small vectors that are two or three dimensional and are used for positions, velocities and so forth. Then there are the vectors we use for example when solving linear equation systems. These are often way larger and their dimension isn't necessarily known beforehand. Simply put, there are the vectors that scale in dimension with the problem size and those that don't. For example the vectors you need to represent positions of units in a game have fixed size. No matter how many units there are in the game, their positions always are two- or three-dimensional. On the other hand if you are calculating splines (which requires you to solve a linear equation system) the vector size in your linear system depends on the amount of points you want to calculate the spline for. This concept of knowing the size beforehand is relevant to programming C++ since you have to decide whether to have the size defined statically at compile time or dynamically at run time. This is also a performance concern since dynamic sizes always require memory allocations while static sizes allow you to place the vector data directly on the call stack and compilers usually have an easier time optimizing for the static case.

In this entry I'll focus on the static/small vector case.

[subheading]So how should my class look like?[/subheading]

Coming from mathematics a lot of people like calling the components of their vectors x, y and z. This leads to the obvious implementation of a vector struct as:

[source lang="cpp"]

struct Vector3d {

Vector3d() : x(0), y(0), z(0) {};

Vector3d(double xx, double yy, double zz) : x(xx), y(yy), z(zz) {};

double x, y, z;

};

[/source]

Then we implement a bunch of overloaded operators and functions like

[source lang="cpp"]

inline Vector3d operator+(const Vector3d &a, const Vector3d &b)

{

return Vector3d(a.x+b.x, a.y+b.y, a.z+b.z);

}

inline double dot(const Vector3d &a, const Vector3d &b)

{

return a.x*b.x + a.y*b.y + a.z*b.z;

}

//...

[/source]

But in C++ there isn't really a reason why we should restrict ourselves to a type and dimension at this point. So we probably should make type and dimension template arguments. For the type this is trivial, we just replace every "double" with "T". But the dimension is a problem with the way we named the components. A obvious solution to allow arbitrary dimensions is to have the components in an array.

[source lang="cpp"]

template

class Vector {

public:

Vector3d() { std::fill(data, data+D, T()); }

inline T& operator[](unsigned i) { return data; }

inline const T& operator[](unsigned i) const { return data; }

private:

T data[D];

};

[/source]

Sadly we lose the ability to access the components via x, y and z. If we really want that we can explicitly specialize the template and add references. For the 3D case this would look like:

[source lang="cpp"]

template

class Vector

public:

Vector3d() : x(data[0]), y(data[1]), z(data[2]) { std::fill(data, data+D, T()); }

inline T& operator[](unsigned i) { return data; }

inline const T& operator[](unsigned i) const { return data; }

T &x, &y, &z;

private:

T data[D];

};

[/source]

I personally don't like doing this since writing myvec[0] instead of myvec.x is just as intuitive to me and it automatically reduces the chances that I'm writing copy/paste code since having the index operator in the sources makes it more obvious where I should use loops. Also as soon as you write a function using x, y and z, chances are you're restricting yourself to writing that function for the 2D or 3D case when you could just as well write it for the general case.

The operators and functions obviously have to be rewritten:

[source lang="cpp"]

template

inline Vector

{

Vector

for(unsigned i = 0;i

result = a + b;

}

return result;

}

template

inline T dot(const Vector

{

T result = 0;

for(unsigned i = 0;i

result += a * b;

}

return result;

}

//...

[/source]

Woha... suddenly we have loops and local variables? Isn't that going to affect performance?

Luckily compilers optimize this kind of stuff. Since D is known at compile time and usually small, the compiler will unroll that loop and after inlineing it will probably also optimize away the local variables if possible

[subheading]Getting rid of temporary objects[/subheading]

There is still room for improvement in the way vector operations are handled. The operators produce temporary objects for intermediate results and that can be expensive when the temporary object is of nontrivial size.

The compiler will very likely turn this

[source lang="cpp"]

Vector

a = b+c+d;

[/source]

into something like

[source lang="cpp"]

Vector

for(unsigned i = 0;i

Vector

for(unsigned i = 0;i

for(unsigned i = 0;i

[/source]

when the whole operation could actually be written as

[source lang="cpp"]

for(unsigned i = 0;i

[/source]

which is obviously preferable.

The solution to this problem are expression templates. The idea is to have operators and functions return a operation object instead of carrying out the operation themselves. So an expression like the one above will be turned into a tree of operation objects which can then be evaluated by the assignment operator.

To achieve this we first need a base class for everything that has "vectorlike behavior":

[source lang="cpp"]

template

class VectorExpr {

public:

inline operator const A&() const //cast operator for convenience

{

return *static_cast

}

};

[/source]

The new vector class now inherits from this and also gets a modified operator= and "copy constructors"

[source lang="cpp"]

template

class Vector : public VectorExpr

public:

static const unsigned Dim = D;

typedef T Type;

Vector() { std::fill(data, data+Dim, T()); }

template

Vector(const VectorExpr

{

const A& ao ( a );

for(unsigned i = 0;i

}

template

Vector& operator=(const VectorExpr

{

const A& ao ( a );

for(unsigned i = 0;i

return *this;

}

T& operator[] (unsigned i) { return data; }

const T& operator[] (unsigned i) const { return data; }

private:

T data[Dim];

};

[/source]

And finally the operators look like

[source lang="cpp"]

template

class AddOp : public VectorExpr

public:

AddOp(const A& pa, const B& pb) : a(pa), b(pb) { }

inline T operator[](unsigned i) const { return a + b; }

private:

const A& a;

const B& b;

};

template

inline AddOp

operator+(const VectorExpr

{

return AddOp

}

[/source]

The whole AddOp object will get optimized away at compile time and the result is the nice loop we wanted without any temporary objects (apart from the operation objects that aren't actually constructed at runtime).

And here is the full source code of the vector class I'm using at the moment:

[source lang="cpp"]

/*

* MathVector.h - Copyright (c) 2011 Jakob Progsch

*

* This software is provided 'as-is', without any express or implied

* warranty. In no event will the authors be held liable for any damages

* arising from the use of this software.

*

* Permission is granted to anyone to use this software for any purpose,

* including commercial applications, and to alter it and redistribute it

* freely, subject to the following restrictions:

*

* 1. The origin of this software must not be misrepresented; you must not

* claim that you wrote the original software. If you use this software

* in a product, an acknowledgment in the product documentation would be

* appreciated but is not required.

*

* 2. Altered source versions must be plainly marked as such, and must not be

* misrepresented as being the original software.

*

* 3. This notice may not be removed or altered from any source

* distribution.

*/

/*

* MathVector.h provides a simple static vector template class with

* a basic expression template ansatz for vector operations.

*/

#ifndef VECTOR_H

#define VECTOR_H

#include

#include

#include

#include

template

//base class for all expression templates

template

class VectorExpr {

public:

inline operator const A&() const

{

return *static_cast

}

};

//better use macros instead of copy pasting this stuff all over the place

#define MAKE_VEC_VEC_EXPRESSION(NAME, EXPR, FUNCTION) \

template

class NAME : public VectorExpr

public: \

NAME(const A& pa, const B& pb) : a(pa), b(pb) { } \

inline T operator[](unsigned i) const { return EXPR; } \

private: \

const A& a; \

const B& b; \

}; \

\

template

inline NAME

FUNCTION(const VectorExpr

{ \

return NAME

}

#define MAKE_VEC_SCAL_EXPRESSION(NAME, EXPR, FUNCTION) \

template

class NAME : public VectorExpr

public: \

NAME(const A& pa, const T& pb) : a(pa), b(pb) { } \

inline T operator[](unsigned i) const { return EXPR; } \

private: \

const A& a; \

const T& b; \

}; \

\

template

inline NAME

FUNCTION(const VectorExpr

{ \

return NAME

}

#define MAKE_SCAL_VEC_EXPRESSION(NAME, EXPR, FUNCTION) \

template

class NAME : public VectorExpr

public: \

NAME(const T& pa, const A& pb) : a(pa), b(pb) { } \

inline T operator[](unsigned i) const { return EXPR; } \

private: \

const T& a; \

const A& b; \

}; \

\

template

inline NAME

FUNCTION(const T &a, const VectorExpr

{ \

return NAME

}

#define MAKE_VEC_EXPRESSION(NAME, EXPR, FUNCTION) \

template

class NAME : public VectorExpr

public: \

NAME(const A& pa) : a(pa) { } \

inline T operator[](unsigned i) const { return EXPR; } \

private: \

const A& a; \

}; \

\

template

inline NAME

FUNCTION(const VectorExpr

{ \

return NAME

}

//create actual functions and operators

MAKE_VEC_VEC_EXPRESSION (EMulExpr, a * b, multiply_elements)

MAKE_VEC_VEC_EXPRESSION (EDivExpr, a / b, divide_elements)

MAKE_VEC_VEC_EXPRESSION (AddExpr, a + b, operator+)

MAKE_VEC_VEC_EXPRESSION (SubExpr, a - b, operator-)

MAKE_VEC_SCAL_EXPRESSION(DivExpr, a / b, operator/)

MAKE_VEC_SCAL_EXPRESSION(MulExpr1, a * b, operator*)

MAKE_SCAL_VEC_EXPRESSION(MulExpr2, a * b, operator*)

MAKE_VEC_EXPRESSION (NegExpr, -a, operator-)

MAKE_VEC_VEC_EXPRESSION (EMaxExpr, std::max(a,b), max)

MAKE_VEC_VEC_EXPRESSION (EMinExpr, std::min(a,b), min)

//sub vector proxy

template

class SubVectorExpr : public VectorExpr

public:

SubVectorExpr(const A& pa, const unsigned& o)

: a(pa), offset(o) { }

inline T operator[](unsigned i) const

{

return a[i+offset];

}

private:

const A& a;

const unsigned offset;

};

template

inline SubVectorExpr

SubVector(const VectorExpr

{

return SubVectorExpr

}

//actual vector class

template

class Vector : public VectorExpr

public:

static const unsigned Dim = D;

typedef T Type;

Vector()

{

std::fill(data, data+Dim, T());

}

Vector(const T &p1)

{

BOOST_STATIC_ASSERT(Dim==1);

data[0] = p1;

}

Vector(const T &p1, const T &p2)

{

BOOST_STATIC_ASSERT(Dim==2);

data[0] = p1;

data[1] = p2;

}

Vector(const T &p1, const T &p2, const T &p3)

{

BOOST_STATIC_ASSERT(Dim==3);

data[0] = p1;

data[1] = p2;

data[2] = p3;

}

Vector(const T &p1, const T &p2, const T &p3, const T &p4)

{

BOOST_STATIC_ASSERT(Dim==4);

data[0] = p1;

data[1] = p2;

data[2] = p3;

data[3] = p4;

}

T* raw() {

return data;

}

template

Vector(const VectorExpr

{

const A& ao ( a );

for(unsigned i = 0; i

}

T& operator[] (unsigned i) {

return data;

}

const T& operator[] (unsigned i) const {

return data;

}

const Vector& operator*=(const T &b)

{

for(unsigned i = 0; i

return *this;

}

const Vector& operator/=(const T &b)

{

for(unsigned i = 0; i

return *this;

}

template

const Vector& operator+=(const VectorExpr

{

const A& ao ( a );

for(unsigned i = 0; i

return *this;

}

template

const Vector& operator-=(const VectorExpr

{

const A& ao ( a );

for(unsigned i = 0; i

return *this;

}

template

Vector& operator=(const VectorExpr

{

const A& ao ( a );

for(unsigned i = 0; i

return *this;

}

Vector& normalize()

{

*this /= abs(*this);

return *this;

}

private:

T data[Dim];

};

//actual InPlaceVector class

template

class InPlaceVector : public VectorExpr

public:

static const unsigned Dim = D;

typedef T Type;

InPlaceVector(T *const d) : data(d)

{ }

T* raw() {

return data;

}

T& operator[] (unsigned i) {

return data;

}

const T& operator[] (unsigned i) const {

return data;

}

const InPlaceVector& operator*=(const T &b)

{

for(unsigned i = 0; i

return *this;

}

const InPlaceVector& operator/=(const T &b)

{

for(unsigned i = 0; i

return *this;

}

template

const InPlaceVector& operator+=(const VectorExpr

{

const A& ao ( a );

for(unsigned i = 0; i

return *this;

}

template

const InPlaceVector& operator-=(const VectorExpr

{

const A& ao ( a );

for(unsigned i = 0; i

return *this;

}

template

InPlaceVector& operator=(const VectorExpr

{

const A& ao ( a );

for(unsigned i = 0; i

return *this;

}

InPlaceVector& normalize()

{

*this /= abs(*this);

return *this;

}

private:

T *const data;

};

template

InPlaceVector

{

return InPlaceVector

}

template

inline Vector

{

return Vector

}

//"reduction" functions that don't return expression templates

template

inline T sum(const VectorExpr

{

const A& ao ( a );

T res = 0;

for(unsigned i = 0; i

return res;

}

template

inline T dot(const VectorExpr

{

return sum(multiply_elements(a, b));

}

template

inline T squared_norm(const VectorExpr

{

const A& ao ( a );

T res = 0;

for(unsigned i = 0; i

T tmp = ao;

res += tmp*tmp;

}

return res;

}

template

inline T abs(const VectorExpr

{

return std::sqrt(squared_norm(a));

}

template

inline T norm(const VectorExpr

{

return std::sqrt(squared_norm(a));

}

template

inline Vector

{

return eval(a/norm(a));

}

template

std::ostream& operator<<(std::ostream& out, const VectorExpr

{

const A& ao ( a );

out << '(' << ao[0];

for(unsigned i = 1; i

out << ", " << ao;

}

out << ')';

return out;

}

#undef MAKE_VEC_VEC_EXPRESSION

#undef MAKE_VEC_SCAL_EXPRESSION

#undef MAKE_SCAL_VEC_EXPRESSION

#undef MAKE_VEC_EXPRESSION

#endif

[/source]

Suggestions and critique appreciated.

0

Sign in to follow this

Followers
0

## 0 Comments

There are no comments to display.

## Create an account or sign in to comment

You need to be a member in order to leave a comment

## Create an account

Sign up for a new account in our community. It's easy!

Register a new account

## Sign in

Already have an account? Sign in here.

Sign In Now