Template Specialization

#ifndef VECTOR3_H
#define VECTOR3_H

#include <ostream>
#include <cassert>

/** Class representing a vector of arbitrary size and type.
  * Provides all the functions necessary to operate on Vectors.
  * The * operator is overloaded to provide the dot product. The % operator is 
  * overloaded to provide the cross product.
  * Internally all variables are represented in a statically allocated array.
  */
template<typename T, int N>
class Vector
{
  public:
    /** Default Constructor.
      * Initializes all values to zero.
      */
    Vector()
    {
      for(int i = 0; i < N; i++)
        data = 0;
    }

    // Copy and assignment constructors autogenerated by the compiler.
    
    /// Compound addition operator
    const Vector<T,N>& operator+=(const Vector<T,N> &rhs)
    {
      for(int i = 0; i < N; i++)
        data += rhs.data;
    }

    /// Compound subtraction operator
    const Vector<T,N>& operator-=(const Vector<T,N> &rhs)
    {
      for(int i = 0; i < N; i++)
        data -= rhs.data;
    }
    /// Compound scalar multiplaction operator
    const Vector<T,N>& operator*=(const float scal)
    {
      for(int i = 0; i < N; i++)
        data *= scal;
    }
    /// Compound scalar division operator
    const Vector<T,N>& operator/=(const float scal)
    {
      for(int i = 0; i < N; i++)
        data /= scal;
    }

    /// Addition operator
    inline const Vector<T,N> operator+(const Vector<T,N> &rhs) const
    {
      return Vector<T,N>(*this) += rhs;
    }
    /// Subtraction operator
    inline const Vector<T,N> operator-(const Vector<T,N> &rhs) const
    {
      return Vector<T,N>(*this) -= rhs;
    }
    /// Scalar multiplication operator
    inline const Vector<T,N> operator*(const float scal) const
    {
      return Vector<T,N>(*this) *= scal;
    }
    /// Scalar division operator
    inline const Vector<T,N> operator/(const float scal) const
    {
      return Vector<T,N>(*this) /= scal;
    }

    /// Unary minus operator
    const Vector<T,N> operator-() const
    {
      return Vector<T,N>(*this) *= -1;
    }

    /// Const operator[] for accessing elements in non assignment contexts
    inline const float & operator[](int index) const
    {
      assert(index > -1 && index < N);
      return data[index];
    }
    /// Operator[] for accessing elements in assignment contexts
    inline float & operator[](int index)
    {
      assert(index > -1 && index < N);
      return data[index];
    }
    /// Dot product operator
    const float operator*(const Vector<T,N>& rhs) const
    {
      T result =  0;
      for(int i = 0; i < N; i++)
        result += data * rhs.data;
      return result;
    }
    /** Cross product operator
      * This is not implemented in the template class.
      * This is due to the fact that not all vector spaces allow the cross prod-
      * uct.  The solution to this is to use template specialization to create 
      * methods for spaces that do contain the cross product.
      */
    const Vector<T,N> operator%(const Vector<T,N>& rhs) const
    {
      assert( false && "A non specialization of the cross product has been cal"
                       "led.  This should never happen!" );
    }

    /** The magnitude squared.
      * This function is much faster than magnitude as it avoids a costly call 
      * to sqrt.
      * @return The magnitude^2.
      */
    float magnitude2() const
    {
      // The magnitude^2 is the dot product of a vector with itself.
      return (*this) * (*this);
    }

    /** The magnitude of this vector.
      * @return The magnitude of this vector.
      */
    inline float magnitude() const
    {
      return sqrtf( (*this) * (*this) );
    }

    /** Normalizes this vector.
      * This is done IN PLACE so this vector will be modified.
      * Use the function normal() if you just want the normalized vector and 
      * this vector to remain intact.  This function is faster than normal().
      * @return A const reference to this vector.
      */
    const Vector<T,N> & normalize()
    {
      *this /= magnitude();
      return *this;
    }

    /** Return the normal of this vector.
      * @return The normal of this vector.
      */
    inline const Vector<T,N> normal() const
    {
      return *this / magnitude();
    }


  private:
    T data[N];
};

/// Scalar * vector operator
template<typename T, int N>
const Vector<T,N> operator*(const float &scal, const Vector<T,N> &rhs)
{
  return rhs * scal;
}

/// Insertion operator
template<typename T, int N>
std::ostream& operator<<(std::ostream &os, const Vector<T,N> &vec)
{
  os << '(';
  for(int i = 0; i < N-1; i++)
    os << vec << ' ';
  os << vec[N-1] << ')';
  return os;
}

// Typedefs
typedef Vector<float,3> Vector3F;

#endif

/// Specializations
template<typename T>
class Vector<T,3>
{
  public:
    Vector()
    {
      data[0] = data[1] = data[2] = 0;
    }

    Vector(T x, T y, T z)
    {
      data[0] = x;
      data[1] = y;
      data[2] = z;
    }

  private:
    T data[3];
};

template<typename T>
const Vector<T,3> Vector<T,3>::operator%(const Vector<T,3> &rhs) const
{
  T i = (data[1] * rhs[2]) - (data[2] * rhs[1]);
  T j = (data[0] * rhs[2]) - (data[2] * rhs[0]);
  T k = (data[0] * rhs[1]) - (data[1] * rhs[0]);;

  return Vector<T,3>(i,j,k);
}