/*========================================================================
    NVector.h : Vector class for general constant size vectors
  ------------------------------------------------------------------------
    Copyright (C) 2003  James R. Bruce
    School of Computer Science, Carnegie Mellon University
  ------------------------------------------------------------------------
    This software is distributed under the GNU General Public License,
    version 2.  If you do not have a copy of this licence, visit
    www.gnu.org, or write: Free Software Foundation, 59 Temple Place,
    Suite 330 Boston, MA 02111-1307 USA.  This program is distributed
    in the hope that it will be useful, but WITHOUT ANY WARRANTY,
    including MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
  ========================================================================*/

#ifndef __NVECTOR_H__
#define __NVECTOR_H__

#include <math.h>
#include "util.h"

namespace Vec {

template <class num,const int dim>
struct nvector{
  num v[dim];

  num &operator[](int i)
    {return(v[i]);}
};

template <class num,class num2,const int dim>
void set(nvector<num,dim> &r,num2 val)
{
  switch(dim){
    case 4: r.v[3] = val;
    case 3: r.v[2] = val;
    case 2: r.v[1] = val;
    case 1: r.v[0] = val;
      break;
    default:
      for(int i=0; i<dim; i++){
        r.v[i] = val;
      }
  }
}

template <class num,const int dim>
void copy(nvector<num,dim> &r,const nvector<num,dim> &a)
{
  switch(dim){
    case 4: r.v[3] = a.v[3];
    case 3: r.v[2] = a.v[2];
    case 2: r.v[1] = a.v[1];
    case 1: r.v[0] = a.v[0];
      break;
    default:
      for(int i=0; i<dim; i++){
        r.v[i] = a.v[i];
      }
  }
}

template <class num,const int dim>
void zero(nvector<num,dim> &r)
{
  set(r,0);
}

template <class num,const int dim>
void unit(nvector<num,dim> &r,int dir)
{
  switch(dim){
    case 4: r.v[3] = (dir == 3);
    case 3: r.v[2] = (dir == 2);
    case 2: r.v[1] = (dir == 1);
    case 1: r.v[0] = (dir == 0);
      break;
    default:
      for(int i=0; i<dim; i++){
        r.v[i] = (dir == i);
      }
  }
}

template <class num,const int dim>
num sqlength(const nvector<num,dim> &a)
{
  num sum = 0;
  switch(dim){
    case 4: sum += a.v[3] * a.v[3];
    case 3: sum += a.v[2] * a.v[2];
    case 2: sum += a.v[1] * a.v[1];
    case 1: sum += a.v[0] * a.v[0];
      break;
    default:
      for(int i=0; i<dim; i++){
        sum += a.v[i] * a.v[i];
      }
  }
  return(sum);
}

template <class num,const int dim>
num length(const nvector<num,dim> &a)
{
  return(sqrt(sqlength(a)));
}

template <class num,const int dim>
num sqdist(const nvector<num,dim> &a,const nvector<num,dim> &b)
{
  num sum = 0;
  switch(dim){
    case 4: sum += sq(a.v[3] - b.v[3]);
    case 3: sum += sq(a.v[2] - b.v[2]);
    case 2: sum += sq(a.v[1] - b.v[1]);
    case 1: sum += sq(a.v[0] - b.v[0]);
      break;
    default:
      for(int i=0; i<dim; i++){
        sum += sq(a.v[i] - b.v[i]);
      }
  }
  return(sum);
}

template <class num,const int dim>
num dist(const nvector<num,dim> &a,const nvector<num,dim> &b)
{
  return(sqrt(sqdist(a,b)));
}

template <class num,const int dim>
num dot(const nvector<num,dim> &a,const nvector<num,dim> &b)
{
  num sum = 0;
  switch(dim){
    case 4: sum += a.v[3] * b.v[3];
    case 3: sum += a.v[2] * b.v[2];
    case 2: sum += a.v[1] * b.v[1];
    case 1: sum += a.v[0] * b.v[0];
      break;
    default:
      for(int i=0; i<dim; i++){
        sum += a.v[i] * b.v[i];
      }
  }
  return(sum);
}


#define NVECTOR_EQUAL_BINARY_OPERATOR(fname,opr) \
  template <class num,const int dim>    \
  void fname(nvector<num,dim> &r,       \
             const nvector<num,dim> &a) \
  {                               \
    switch(dim){                  \
      case 4: r.v[3] opr a.v[3];  \
      case 3: r.v[2] opr a.v[2];  \
      case 2: r.v[1] opr a.v[1];  \
      case 1: r.v[0] opr a.v[0];  \
        break;                    \
      default:                    \
        for(int i=0; i<dim; i++){ \
          r.v[i] opr a.v[i];      \
        }                         \
    }                             \
  }

NVECTOR_EQUAL_BINARY_OPERATOR(add,+=)
NVECTOR_EQUAL_BINARY_OPERATOR(sub,-=)
NVECTOR_EQUAL_BINARY_OPERATOR(mul,*=)
NVECTOR_EQUAL_BINARY_OPERATOR(div,/=)


#define NVECTOR_BINARY_OPERATOR(fname,opr) \
  template <class num,const int dim>       \
  void fname(nvector<num,dim> &r,          \
             const nvector<num,dim> &a,    \
             const nvector<num,dim> &b)    \
  {                                        \
    switch(dim){                           \
      case 4: r.v[3] = a.v[3] opr b.v[3];  \
      case 3: r.v[2] = a.v[2] opr b.v[2];  \
      case 2: r.v[1] = a.v[1] opr b.v[1];  \
      case 1: r.v[0] = a.v[0] opr b.v[0];  \
        break;                             \
      default:                             \
        for(int i=0; i<dim; i++){          \
          r.v[i] = a.v[i] opr b.v[i];      \
        }                                  \
    }                                      \
  }

NVECTOR_BINARY_OPERATOR(add,+)
NVECTOR_BINARY_OPERATOR(sub,-)
NVECTOR_BINARY_OPERATOR(mul,*)
NVECTOR_BINARY_OPERATOR(div,/)


#define NVECTOR_SCALAR_BINARY_OPERATOR(fname,opr) \
  template <class num,class num2,const int dim>   \
  void fname(nvector<num,dim> &r,       \
             const nvector<num,dim> &a, \
             num2 s)                    \
  {                                     \
    switch(dim){                        \
      case 4: r.v[3] = a.v[3] opr s;    \
      case 3: r.v[2] = a.v[2] opr s;    \
      case 2: r.v[1] = a.v[1] opr s;    \
      case 1: r.v[0] = a.v[0] opr s;    \
        break;                          \
      default:                          \
        for(int i=0; i<dim; i++){       \
          r.v[i] = a.v[i] opr s;        \
        }                               \
    }                                   \
  }

NVECTOR_SCALAR_BINARY_OPERATOR(mul,*)
NVECTOR_SCALAR_BINARY_OPERATOR(div,/)


#define NVECTOR_SCALAR_EQUAL_BINARY_OPERATOR(fname,opr) \
  template <class num,class num2,const int dim>         \
  void fname(nvector<num,dim> &r, num2 s) \
  {                                       \
    switch(dim){                          \
      case 4: r.v[3] opr s;               \
      case 3: r.v[2] opr s;               \
      case 2: r.v[1] opr s;               \
      case 1: r.v[0] opr s;               \
        break;                            \
      default:                            \
        for(int i=0; i<dim; i++){         \
          r.v[i] opr s;                   \
        }                                 \
    }                                     \
  }

NVECTOR_SCALAR_EQUAL_BINARY_OPERATOR(mul,*=)
NVECTOR_SCALAR_EQUAL_BINARY_OPERATOR(div,/=)


template <class num,const int dim>
void neg(nvector<num,dim> &r,const nvector<num,dim> &a)
{
  mul(r,a,-1);
}

template <class num,const int dim>
void neg(nvector<num,dim> &r)
{
  mul(r,-1);
}

template <class num,const int dim>
void norm(nvector<num,dim> &r,const nvector<num,dim> &a)
{
  num l = length(a);
  mul(r,a,1/l);
}

template <class num,const int dim>
void norm(nvector<num,dim> &r)
{
  num l = length(r);
  mul(r,1/l);
}

}; // namespace Vec

#endif