Add dynamically allocated versions of matrix and vector

Sometimes it's not possible to know the dimensions at compile time, for
example when calculating regressions from user-inputted data.

diff --git a/source/linal/dummy.cpp b/source/linal/dummy.cpp
index 05e6ff2c39de45d549cd6346b34af29b8caa9568..84d7395c2478da62f682b9f1839c897a073609f1 100644 (file)
--- a/source/linal/dummy.cpp
+++ b/source/linal/dummy.cpp
@@ -1,6 +1,8 @@
 #include "vector.h"
 #include "matrix.h"
 #include "squarematrix.h"
+#include "dynamicvector.h"
+#include "dynamicmatrix.h"
 
 /* This dummy file is needed to make Builder happy, and also serves as a syntax
 check by pulling in all the headers. */

diff --git a/source/linal/dynamicmatrix.h b/source/linal/dynamicmatrix.h
new file mode 100644 (file)
index 0000000..3852424
--- /dev/null
+++ b/source/linal/dynamicmatrix.h
@@ -0,0 +1,323 @@
+#ifndef MSP_LINAL_DYNAMICMATRIX_H_
+#define MSP_LINAL_DYNAMICMATRIX_H_
+
+#include <algorithm>
+#include "dynamicvector.h"
+#include "matrixops.h"
+
+namespace Msp {
+namespace LinAl {
+
+/**
+A general mathematical matrix.  Size is specified at runtime.  May be slower
+than a fixed-size matrix.  There are no compile-time diagnostics for mismatched
+matrix sizes.
+*/
+template<typename T>
+class DynamicMatrix
+{
+public:
+       typedef T ElementType;
+
+private:
+       unsigned rows_;
+       unsigned columns_;
+       T *data;
+
+public:
+       DynamicMatrix(unsigned, unsigned);
+       DynamicMatrix(unsigned, unsigned, const T *);
+       DynamicMatrix(const DynamicMatrix &);
+       DynamicMatrix &operator=(const DynamicMatrix &);
+       ~DynamicMatrix();
+
+       unsigned rows() const { return rows_; }
+       unsigned columns() const { return columns_; }
+
+       T &element(unsigned, unsigned);
+       const T &element(unsigned, unsigned) const;
+       T &operator()(unsigned i, unsigned j) { return element(i, j); }
+       const T &operator()(unsigned i, unsigned j) const { return element(i, j); }
+
+       DynamicMatrix &operator*=(T);
+       DynamicMatrix &operator/=(T);
+       DynamicMatrix &operator+=(const DynamicMatrix &);
+       DynamicMatrix &operator-=(const DynamicMatrix &);
+
+       DynamicMatrix &exchange_rows(unsigned, unsigned);
+       DynamicMatrix &multiply_row(unsigned, T);
+       DynamicMatrix &add_row(unsigned, unsigned, T);
+
+       DynamicMatrix &invert();
+};
+
+
+template<typename T>
+inline DynamicMatrix<T>::DynamicMatrix(unsigned r, unsigned c):
+       rows_(r),
+       columns_(c),
+       data(new T[rows_*columns_])
+{
+       std::fill(data, data+rows_*columns_, T());
+}
+
+template<typename T>
+inline DynamicMatrix<T>::DynamicMatrix(unsigned r, unsigned c, const T *d):
+       rows_(r),
+       columns_(c),
+       data(new T[rows_*columns_])
+{
+       std::copy(d, d+rows_*columns_, data);
+}
+
+template<typename T>
+inline DynamicMatrix<T>::DynamicMatrix(const DynamicMatrix &other):
+       rows_(other.rows()),
+       columns_(other.columns()),
+       data(new T[rows_*columns_])
+{
+       std::copy(other.data, other.data+rows_*columns_, data);
+}
+
+template<typename T>
+inline DynamicMatrix<T> &DynamicMatrix<T>::operator=(const DynamicMatrix &other)
+{
+       if(rows_!=other.rows() || columns_!=other.columns())
+       {
+               delete[] data;
+               rows_ = other.rows();
+               columns_ = other.columns();
+               data = new T[rows_*columns_];
+       }
+
+       std::copy(other.data, other.data+rows_*columns_, data);
+
+       return *this;
+}
+
+template<typename T>
+inline DynamicMatrix<T>::~DynamicMatrix()
+{
+       delete[] data;
+}
+
+template<typename T>
+inline T &DynamicMatrix<T>::element(unsigned i, unsigned j)
+{
+       if(i>=rows_ || j>=columns_)
+               throw std::out_of_range("DynamicMatrix::element");
+
+       return data[i+rows_*j];
+}
+
+template<typename T>
+inline const T &DynamicMatrix<T>::element(unsigned i, unsigned j) const
+{
+       if(i>=rows_ || j>=columns_)
+               throw std::out_of_range("DynamicMatrix::element");
+
+       return data[i+rows_*j];
+}
+
+template<typename T>
+inline DynamicMatrix<T> &DynamicMatrix<T>::operator*=(T s)
+{
+       for(unsigned i=0; i<rows_*columns_; ++i)
+               data[i] *= s;
+       return *this;
+}
+
+template<typename T>
+inline DynamicMatrix<T> operator*(const DynamicMatrix<T> &m, T s)
+{
+       DynamicMatrix<T> r(m);
+       return r *= s;
+}
+
+template<typename T>
+inline DynamicMatrix<T> operator*(T s, const DynamicMatrix<T> &m)
+{
+       return m*s;
+}
+
+template<typename T>
+inline DynamicMatrix<T> operator*(const DynamicMatrix<T> &m1, const DynamicMatrix<T> &m2)
+{
+       if(m1.columns()!=m2.rows())
+               throw size_mismatch("matrix*matrix");
+
+       DynamicMatrix<T> r(m1.rows(), m2.columns());
+       for(unsigned i=0; i<m1.rows(); ++i)
+               for(unsigned j=0; j<m2.columns(); ++j)
+                       for(unsigned k=0; k<m1.columns(); ++k)
+                               r(i, j) += m1(i, k)*m2(k, j);
+
+       return r;
+}
+
+template<typename T>
+inline DynamicVector<T> operator*(const DynamicMatrix<T> &m, const DynamicVector<T> &v)
+{
+       if(m.columns()!=v.size())
+               throw size_mismatch("matrix*vector");
+
+       DynamicVector<T> r(m.rows());
+       for(unsigned i=0; i<m.rows(); ++i)
+               for(unsigned j=0; j<m.columns(); ++j)
+                       r[i] += m(i, j)*v[j];
+
+       return r;
+}
+
+template<typename T>
+inline DynamicVector<T> operator*(const DynamicVector<T> &v, const DynamicMatrix<T> &m)
+{
+       if(v.size()!=m.rows())
+               throw size_mismatch("vector*matrix");
+
+       DynamicVector<T> r(m.columns());
+       for(unsigned j=0; j<m.columns(); ++j)
+               for(unsigned i=0; i<m.rows(); ++i)
+                       r[i] += v[i]*m(i, j);
+
+       return r;
+}
+
+template<typename T>
+inline DynamicMatrix<T> &DynamicMatrix<T>::operator/=(T s)
+{
+       for(unsigned i=0; i<rows_*columns_; ++i)
+               data[i] /= s;
+       return *this;
+}
+
+template<typename T>
+inline DynamicMatrix<T> operator/(const DynamicMatrix<T> &m, T s)
+{
+       DynamicMatrix<T> r(m);
+       return r /= s;
+}
+
+template<typename T>
+inline DynamicMatrix<T> &DynamicMatrix<T>::operator+=(const DynamicMatrix<T> &m)
+{
+       if(rows_!=m.rows_ || columns_!=m.columns_)
+               throw size_mismatch("matrix+matrix");
+
+       for(unsigned i=0; i<rows_*columns_; ++i)
+               data[i] += m.data[i];
+
+       return *this;
+}
+
+template<typename T>
+inline DynamicMatrix<T> operator+(const DynamicMatrix<T> &m1, const DynamicMatrix<T> &m2)
+{
+       DynamicMatrix<T> r(m1);
+       return r += m2;
+}
+
+template<typename T>
+inline DynamicMatrix<T> &DynamicMatrix<T>::operator-=(const DynamicMatrix<T> &m)
+{
+       if(rows_!=m.rows_ || columns_!=m.columns_)
+               throw size_mismatch("matrix-matrix");
+
+       for(unsigned i=0; i<rows_*columns_; ++i)
+               data[i] += m.data[i];
+
+       return *this;
+}
+
+template<typename T>
+inline DynamicMatrix<T> operator-(const DynamicMatrix<T> &m1, const DynamicMatrix<T> &m2)
+{
+       DynamicMatrix<T> r(m1);
+       return r += m2;
+}
+
+template<typename T>
+inline bool operator==(const DynamicMatrix<T> &m1, const DynamicMatrix<T> &m2)
+{
+       if(m1.rows()!=m2.rows() || m1.columns()!=m2.columns())
+               throw size_mismatch("matrix==matrix");
+
+       for(unsigned i=0; i<m1.rows(); ++i)
+               for(unsigned j=0; j<m1.columns(); ++j)
+                       if(m1(i, j)!=m2(i, j))
+                               return false;
+
+       return true;
+}
+
+template<typename T>
+inline DynamicMatrix<T> &DynamicMatrix<T>::exchange_rows(unsigned i, unsigned j)
+{
+       if(i>=rows_ || j>=rows_)
+               throw std::out_of_range("DynamicMatrix::exchange_rows");
+
+       using std::swap;
+       for(unsigned k=0; k<columns_; ++k)
+               swap(element(i, k), element(j, k));
+
+       return *this;
+}
+
+template<typename T>
+inline DynamicMatrix<T> &DynamicMatrix<T>::multiply_row(unsigned i, T s)
+{
+       if(i>=rows_)
+               throw std::out_of_range("DynamicMatrix::multiply_row");
+
+       for(unsigned k=0; k<columns_; ++k)
+               element(i, k) *= s;
+
+       return *this;
+}
+
+template<typename T>
+inline DynamicMatrix<T> &DynamicMatrix<T>::add_row(unsigned i, unsigned j, T s)
+{
+       if(i>=rows_ || j>=rows_)
+               throw std::out_of_range("DynamicMatrix::exchange_rows");
+
+       for(unsigned k=0; k<columns_; ++k)
+               element(j, k) += element(i, k)*s;
+
+       return *this;
+}
+
+template<typename T>
+inline DynamicMatrix<T> transpose(const DynamicMatrix<T> &m)
+{
+       DynamicMatrix<T> r(m.columns(), m.rows());
+       for(unsigned i=0; i<m.rows(); ++i)
+               for(unsigned j=0; j<m.columns(); ++j)
+                       r(j, i) = m(i, j);
+       return r;
+}
+
+template<typename T>
+inline DynamicMatrix<T> &DynamicMatrix<T>::invert()
+{
+       if(rows_!=columns_)
+               throw size_mismatch("matrix invert");
+
+       DynamicMatrix<T> r(rows_, columns_);
+       for(unsigned i=0; i<rows_; ++i)
+               r(i, i) = T(1);
+
+       return invert_matrix(*this, r);
+}
+
+template<typename T>
+inline DynamicMatrix<T> invert(const DynamicMatrix<T> &m)
+{
+       DynamicMatrix<T> r = m;
+       return r.invert();
+}
+
+} // namespace LinAl
+} // namespace Msp
+
+#endif

diff --git a/source/linal/dynamicvector.h b/source/linal/dynamicvector.h
new file mode 100644 (file)
index 0000000..8839c33
--- /dev/null
+++ b/source/linal/dynamicvector.h
@@ -0,0 +1,249 @@
+#ifndef MSP_LINAL_DYNAMICVECTOR_H_
+#define MSP_LINAL_DYNAMICVECTOR_H_
+
+#include <algorithm>
+#include <cmath>
+#include <stdexcept>
+
+namespace Msp {
+namespace LinAl {
+
+class size_mismatch: public std::logic_error
+{
+public:
+       size_mismatch(const std::string &w): std::logic_error(w) { }
+       virtual ~size_mismatch() throw() { }
+};
+
+
+/**
+A general mathematical vector.  Size is specified at runtime.  May be slower
+than a fixed-size vector.  There are no compile-time diagnostics for mismatched
+vector sizes.
+*/
+template<typename T>
+class DynamicVector
+{
+public:
+       typedef T ElementType;
+
+private:
+       unsigned size_;
+       T *data;
+
+public:
+       DynamicVector(unsigned);
+       DynamicVector(unsigned, const T *);
+       DynamicVector(const DynamicVector &);
+       DynamicVector &operator=(const DynamicVector &);
+       ~DynamicVector();
+
+       unsigned size() const { return size_; }
+
+       T &operator[](unsigned);
+       const T &operator[](unsigned) const;
+
+       DynamicVector &operator*=(T);
+       DynamicVector &operator/=(T);
+       DynamicVector &operator+=(const DynamicVector &);
+       DynamicVector &operator-=(const DynamicVector &);
+
+       T norm() const;
+       DynamicVector &normalize();
+};
+
+template<typename T>
+inline DynamicVector<T>::DynamicVector(unsigned s):
+       size_(s),
+       data(new T[size_])
+{
+       std::fill(data, data+size_, T());
+}
+
+template<typename T>
+inline DynamicVector<T>::DynamicVector(unsigned s, const T *d):
+       size_(s),
+       data(new T[size_])
+{
+       std::copy(d, d+size_, data);
+}
+
+template<typename T>
+inline DynamicVector<T>::DynamicVector(const DynamicVector &other):
+       size_(other.size()),
+       data(new T[size_])
+{
+       std::copy(other.data, other.data+size_, data);
+}
+
+template<typename T>
+inline DynamicVector<T> &DynamicVector<T>::operator=(const DynamicVector<T> &other)
+{
+       if(size_!=other.size())
+       {
+               delete[] data;
+               size_ = other.size();
+               data = new T[size_];
+       }
+
+       std::copy(other.data, other.data+size_, data);
+
+       return *this;
+}
+
+template<typename T>
+inline DynamicVector<T>::~DynamicVector()
+{
+       delete[] data;
+}
+
+template<typename T>
+inline T &DynamicVector<T>::operator[](unsigned i)
+{
+       if(i>=size_)
+               throw std::out_of_range("DynamicVector::operator[]");
+
+       return data[i];
+}
+
+template<typename T>
+inline const T &DynamicVector<T>::operator[](unsigned i) const
+{
+       if(i>=size_)
+               throw std::out_of_range("DynamicVector::operator[]");
+
+       return data[i];
+}
+
+template<typename T>
+inline DynamicVector<T> &DynamicVector<T>::operator*=(T s)
+{
+       for(unsigned i=0; i<size_; ++i)
+               data[i] *= s;
+       return *this;
+}
+
+template<typename T>
+inline DynamicVector<T> operator*(const DynamicVector<T> &v, T s)
+{
+       DynamicVector<T> r(v);
+       return r *= s;
+}
+
+template<typename T>
+inline DynamicVector<T> operator*(T s, const DynamicVector<T> &v)
+{
+       return v*s;
+}
+
+template<typename T>
+inline DynamicVector<T> &DynamicVector<T>::operator/=(T s)
+{
+       for(unsigned i=0; i<size_; ++i)
+               data[i] /= s;
+       return *this;
+}
+
+template<typename T>
+inline DynamicVector<T> operator/(const DynamicVector<T> &v, T s)
+{
+       DynamicVector<T> r(v);
+       return r /= s;
+}
+
+template<typename T>
+inline DynamicVector<T> &DynamicVector<T>::operator+=(const DynamicVector<T> &v)
+{
+       if(size_!=v.size())
+               throw size_mismatch("vector+vector");
+
+       for(unsigned i=0; i<size_; ++i)
+               data[i] += v.data[i];
+
+       return *this;
+}
+
+template<typename T>
+inline DynamicVector<T> operator+(const DynamicVector<T> &v1, const DynamicVector<T> &v2)
+{
+       DynamicVector<T> r(v1);
+       return r += v2;
+}
+
+template<typename T>
+inline DynamicVector<T> &DynamicVector<T>::operator-=(const DynamicVector<T> &v)
+{
+       if(size_!=v.size())
+               throw size_mismatch("vector-vector");
+
+       for(unsigned i=0; i<size_; ++i)
+               data[i] -= v.data[i];
+
+       return *this;
+}
+
+template<typename T>
+inline DynamicVector<T> operator-(const DynamicVector<T> &v1, const DynamicVector<T> &v2)
+{
+       DynamicVector<T> r(v1);
+       return r -= v2;
+}
+
+template<typename T>
+inline DynamicVector<T> operator-(const DynamicVector<T> &v)
+{
+       DynamicVector<T> r(v);
+       for(unsigned i=0; i<r.size(); ++i)
+               r[i] = -r[i];
+       return r;
+}
+
+template<typename T>
+inline bool operator==(const DynamicVector<T> &v1, const DynamicVector<T> &v2)
+{
+       if(v1.size()!=v2.size())
+               throw size_mismatch("vector==vector");
+
+       for(unsigned i=0; i<v1.size(); ++i)
+               if(v1[i]!=v2[i])
+                       return false;
+
+       return true;
+}
+
+template<typename T>
+inline T inner_product(const DynamicVector<T> &v1, const DynamicVector<T> &v2)
+{
+       if(v1.size()!=v2.size())
+               throw size_mismatch("inner_product");
+
+       T r = T();
+       for(unsigned i=0; i<v1.size(); ++i)
+               r += v1[i]*v2[i];
+       return r;
+}
+
+template<typename T>
+inline T DynamicVector<T>::norm() const
+{
+       using std::sqrt;
+       return sqrt(inner_product(*this, *this));
+}
+
+template<typename T>
+inline DynamicVector<T> &DynamicVector<T>::normalize()
+{
+       return *this /= norm();
+}
+
+template<typename T>
+inline DynamicVector<T> normalize(const DynamicVector<T> &v)
+{
+       DynamicVector<T> r(v);
+       return r.normalize();
+}
+
+} // namespace LinAL
+} // namespace Msp
+
+#endif

diff --git a/source/linal/matrix.h b/source/linal/matrix.h
index 929c48247e517aa47d837df1438c94465e25064a..a57ce08d68b78f6d2f49cd3bffbc16fbc223847f 100644 (file)
--- a/source/linal/matrix.h
+++ b/source/linal/matrix.h
@@ -13,6 +13,9 @@ A general mathematical matrix with M rows and N columns.
 template<typename T, unsigned M, unsigned N>
 class Matrix
 {
+public:
+       typedef T ElementType;
+
 private:
        T data[M*N];
 
@@ -25,11 +28,17 @@ public:
        static Matrix from_columns(const Vector<T, M> *);
        static Matrix from_rows(const Vector<T, N> *);
 
+       unsigned rows() const { return M; }
+       unsigned columns() const { return N; }
+
        T &element(unsigned i, unsigned j) { return data[i+M*j]; }
        const T &element(unsigned i, unsigned j) const { return data[i+M*j]; }
        T &operator()(unsigned i, unsigned j) { return element(i, j); }
        const T &operator()(unsigned i, unsigned j) const { return element(i, j); }
 
+       template<unsigned P, unsigned Q>
+       Matrix<T, P, Q> select(const Vector<unsigned, P> &, const Vector<unsigned, Q> &) const;
+
        template<unsigned P, unsigned Q>
        Matrix<T, P, Q> block(unsigned, unsigned) const;
 
@@ -84,6 +93,17 @@ inline Matrix<T, M, N> Matrix<T, M, N>::from_rows(const Vector<T, N> *v)
        return m;
 }
 
+template<typename T, unsigned M, unsigned N>
+template<unsigned P, unsigned Q>
+inline Matrix<T, P, Q> Matrix<T, M, N>::select(const Vector<unsigned, P> &rows, const Vector<unsigned, Q> &cols) const
+{
+       Matrix<T, P, Q> r;
+       for(unsigned j=0; j<P; ++j)
+               for(unsigned i=0; i<Q; ++i)
+                       r(j, i) = element(rows[j], cols[i]);
+       return r;
+}
+
 template<typename T, unsigned M, unsigned N>
 template<unsigned P, unsigned Q>
 inline Matrix<T, P, Q> Matrix<T, M, N>::block(unsigned y, unsigned x) const

diff --git a/source/linal/matrixops.h b/source/linal/matrixops.h
new file mode 100644 (file)
index 0000000..790dfbf
--- /dev/null
+++ b/source/linal/matrixops.h
@@ -0,0 +1,62 @@
+#ifndef MSP_LINAL_MATRIXOPS_H_
+#define MSP_LINAL_MATRIXOPS_H_
+
+#include <cmath>
+#include <stdexcept>
+
+namespace Msp {
+namespace LinAl {
+
+class not_invertible: public std::domain_error
+{
+public:
+       not_invertible(): domain_error(std::string()) { }
+       virtual ~not_invertible() throw() { }
+};
+
+
+template<typename T>
+inline T &invert_matrix(T &m, T &r)
+{
+       typedef typename T::ElementType V;
+       using std::abs;
+
+       for(unsigned i=0; i<m.rows(); ++i)
+       {
+               unsigned pivot = i;
+               for(unsigned j=i+1; j<m.rows(); ++j)
+                       if(abs(m.element(j, i))>abs(m.element(pivot, i)))
+                               pivot = j;
+
+               if(m.element(pivot, i)==V(0))
+                       throw not_invertible();
+
+               if(pivot!=i)
+               {
+                       m.exchange_rows(i, pivot);
+                       r.exchange_rows(i, pivot);
+               }
+
+               for(unsigned j=i+1; j<m.rows(); ++j)
+               {
+                       V a = -m.element(j, i)/m.element(i, i);
+                       m.add_row(i, j, a);
+                       r.add_row(i, j, a);
+               }
+
+               V a = V(1)/m.element(i, i);
+               m.multiply_row(i, a);
+               r.multiply_row(i, a);
+       }
+
+       for(unsigned i=m.rows(); i-->0; )
+               for(unsigned j=i; j-->0; )
+                       r.add_row(i, j, -m.element(j, i));
+
+       return m = r;
+}
+
+} // namespace LinAl
+} // namespace Msp
+
+#endif

diff --git a/source/linal/squarematrix.h b/source/linal/squarematrix.h
index 5117bb3fa189ee4339d0ce52cdbb4d2ff141b83b..0a97e61172b6489c58e91296421e946cfddfb9d8 100644 (file)
--- a/source/linal/squarematrix.h
+++ b/source/linal/squarematrix.h
@@ -2,19 +2,12 @@
 #define MSP_LINAL_SQUAREMATRIX_H_
 
 #include <cmath>
-#include <stdexcept>
 #include "matrix.h"
+#include "matrixops.h"
 
 namespace Msp {
 namespace LinAl {
 
-class not_invertible: public std::domain_error
-{
-public:
-       not_invertible(): domain_error(std::string()) { }
-       virtual ~not_invertible() throw() { }
-};
-
 /**
 A mathematical matrix with S rows and columns.  Some operations are provided
 here that are only possible for square matrices.
@@ -53,42 +46,8 @@ SquareMatrix<T, S> &SquareMatrix<T, S>::operator*=(const SquareMatrix<T, S> &m)
 template<typename T, unsigned S>
 SquareMatrix<T, S> &SquareMatrix<T, S>::invert()
 {
-       using std::abs;
-
        SquareMatrix<T, S> r = identity();
-       for(unsigned i=0; i<S; ++i)
-       {
-               unsigned pivot = i;
-               for(unsigned j=i+1; j<S; ++j)
-                       if(abs(this->element(j, i))>abs(this->element(pivot, i)))
-                               pivot = j;
-
-               if(this->element(pivot, i)==T(0))
-                       throw not_invertible();
-
-               if(pivot!=i)
-               {
-                       this->exchange_rows(i, pivot);
-                       r.exchange_rows(i, pivot);
-               }
-
-               for(unsigned j=i+1; j<S; ++j)
-               {
-                       T a = -this->element(j, i)/this->element(i, i);
-                       this->add_row(i, j, a);
-                       r.add_row(i, j, a);
-               }
-
-               T a = T(1)/this->element(i, i);
-               this->multiply_row(i, a);
-               r.multiply_row(i, a);
-       }
-
-       for(unsigned i=S; i-->0; )
-               for(unsigned j=i; j-->0; )
-                       r.add_row(i, j, -this->element(j, i));
-
-       return *this = r;
+       return invert_matrix(*this, r);
 }
 
 template<typename T, unsigned S>

diff --git a/source/linal/vector.h b/source/linal/vector.h
index 6dfe6ca0167fb3e1d57e65217003e3d700f6a3b3..152c26e66482e943e5837c5ddecc54fc9e7c2fde 100644 (file)
--- a/source/linal/vector.h
+++ b/source/linal/vector.h
@@ -74,6 +74,8 @@ template<typename T, unsigned N>
 class Vector: public VectorComponents<T, N>
 {
 public:
+       typedef T ElementType;
+
        Vector();
        Vector(const T *);
        Vector(T, T);
@@ -82,6 +84,8 @@ public:
        template<typename U>
        Vector(const Vector<U, N> &);
 
+       unsigned size() const { return N; }
+
        template<unsigned M>
        Vector<T, M> slice(unsigned) const;