3D mathematical vector class I've been working on a 3D mathematical vector class which should be as streamlined as possible for use in numerical simulations. It will be used to model 3D-physical vectors.

Here, 3D-vector should be taken in mathematical sense, meaning a tuple (a,b,c).

I hoped to design it in a modern and fast way - but one is never perfect. So, I would be interested in some input from your side.  Any tips for making this faster?

//threevector.h
#ifndef threevector_h_
#define threevector_h_

#include &lt;fstream&gt;
#include &lt;cmath&gt;
#include &lt;array&gt;


template &lt;class T&gt;
class threevector
{
private:
  static const int dim = 3; //dimension of vector
  std::array&lt;T,dim&gt; container;  

public:

  //constructors and assignment
  threevector(const double a = 0, const double b = 0, const double c = 0):
  container({{a,b,c}}) {}; //standard constructor
  threevector(const threevector&amp; a): container(a.container) {}; 
  //copy constructor

// add once gcc 4.7 is used
//   threevector(threevector&amp;&amp; other): threevector() {swap(*this, other);} 
  // move constructor

  threevector&amp; operator=(threevector rhs) //assignment
  {
    swap(*this, rhs);
    return *this;
  }

  void swap(threevector&amp; first, threevector&amp; second) 
  {first.container.swap(second.container);}


  //operators
  threevector&amp; operator+=(const threevector&amp; rhs) 
  {
    container[0] += rhs.container[0];
    container[1] += rhs.container[1];
    container[2] += rhs.container[2];
    return *this;
  }

  threevector&amp; operator-=(const threevector&amp; rhs) 
  {
    *this += -rhs;
    return *this;
  }

  threevector&amp; operator*=(const double rhs) //scalar multiplication assignment
  {
    container[0] *= rhs;
    container[1] *= rhs;
    container[2] *= rhs;
    return *this;
  }

  threevector&amp; operator/=(const double rhs) //scalar division assignment
  {
    *this *= 1./rhs;
    return *this;
  }


  threevector operator+() const //unary plus
  {
    threevector a(*this);
    return a;
  }

  threevector operator-() const //unary minus
  {
    threevector a(*this);
    a *= -1;
    return a;
  }

  T&amp; operator[](const int input) {return container[input];} //access operator

  const T&amp; operator[](const int input) const {return container[input];} 
  // const access operator

  //utility functions
  double abs() const 
  {return sqrt(container[0]*container[0]+container[1]*container[1]+
   container[2]*container[2]);}
  double abs_sq() const {return pow(abs(),2);}
  void reset() { container[0] = 0; container[1] = 0; container[2] = 0;}
};

//output operator
template&lt;class T&gt;
std::ostream&amp; operator&lt;&lt;(std::ostream&amp; os, const threevector&lt;T&gt;&amp; obj)
{
  os &lt;&lt; std::fixed &lt;&lt; "(" &lt;&lt; obj[0] &lt;&lt; "," &lt;&lt; obj[1] &lt;&lt; "," &lt;&lt; obj[2] &lt;&lt; ")";
  return os;
}

//addition operator
template&lt;class T&gt;
inline threevector&lt;T&gt; operator+(threevector&lt;T&gt; lhs, const threevector&lt;T&gt;&amp; rhs)
{
  lhs += rhs;
  return lhs;
}

//subtraction operator
template&lt;class T&gt;
inline threevector&lt;T&gt; operator-(threevector&lt;T&gt; lhs, const threevector&lt;T&gt;&amp; rhs)
{
  lhs -= rhs;
  return lhs;
}

//scalar product
template&lt;class T&gt;
inline double operator*(const threevector&lt;T&gt;&amp; lhs, const threevector&lt;T&gt;&amp; rhs) 
{return lhs[0] * rhs[0] + lhs[1] * rhs[1] + lhs[2] * rhs[2];}


//product with scalar
template&lt;class T&gt;
inline threevector&lt;T&gt; operator*(const double lhs, threevector&lt;T&gt; rhs) 
{
  rhs[0] *= lhs;
  rhs[1] *= lhs;
  rhs[2] *= lhs; 
  return rhs;
}

//product with scalar
template&lt;class T&gt;
inline threevector&lt;T&gt; operator*(threevector&lt;T&gt; lhs, const double rhs){return rhs*lhs;}

// scalar division
template&lt;class T&gt;
inline threevector&lt;T&gt; operator/(threevector&lt;T&gt; lhs, const double rhs){return lhs*(1./rhs);};


#endif
Your design is a little bit confusing. You define your class as a template class, but just about everything takes and returns a double. Generally, you should try and be consistent: will this class work for any numeric-like T? Then make sure it has constructors that take a T instead of a double:

threevector(T a = T(), T b = T(), T c = T())


Similarly for the overloaded operators that multiply/divide by a scalar, and so on.

The name threevector is awkward. At the least it should be three_vector, but that's still a little awkward. vector_3d is more of the "standard" name for this kind of thing.

Your include guards should generally be in all-caps. This is because they are global in scope, and you want to do as much as possible to minimize the possibility of potential name clashes.

#ifndef threevector_h_ // should be THREE_VECTOR_H_
#define threevector_h_ // or VECTOR_3D_H_


You have some semi-colons where they aren't needed:

threevector(const double a = 0, const double b = 0, const double c = 0):
container({{a,b,c}}) {}; //standard constructor
threevector(const threevector&amp; a): container(a.container) {}; // &lt;--- Not needed


Comments like //standard constructor are obvious and don't really need to be there. Further, your copy constructor and copy assignment operator are redundant; the compiler generated ones are sufficient. 

It's often easiest to write operator+ (and -, * etc) in terms of operator+=:

vector_3d operator+(const vector_3d&amp; a, const vector_3d&amp; b)
{
     vector_3d result(a);
     result += b;
     return result;
}


For something this simple, it's unlikely you'll be able to speed it up much. The only potential that I can see is replacing the call to pow() in:

double abs_sq() const {return pow(abs(),2);}


with something like:

T abs_sq() const
{
    const T ab = abs();  // abs() should also return a T
    return ab * ab;
}


which avoids the overhead of a function call to pow. In practice, this is unlikely to make much of a difference, however.