diff --git a/Makefile.am b/Makefile.am
index 51504aac6..aa2b81be0 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -248,7 +248,6 @@ core_SOURCES = \
 	src/core/array.h \
 	src/core/list.h \
 	src/core/vector.h \
-	src/core/vector.imp \
 	src/core/recarray.h \
 	src/core/matrix.h \
 	src/core/matrix.imp \
@@ -258,7 +257,6 @@ core_SOURCES = \
 	src/core/integer.h \
 	src/core/rational.cc \
 	src/core/rational.h \
-	src/core/vector.cc \
 	src/core/matrix.cc \
 	src/core/sqmatrix.cc \
 	src/core/function.cc \
diff --git a/src/core/array.h b/src/core/array.h
index be4923880..223666947 100644
--- a/src/core/array.h
+++ b/src/core/array.h
@@ -81,6 +81,12 @@ template <class T> class Array {
       m_index++;
       return *this;
     }
+    iterator operator++(int)
+    {
+      auto ret = *this;
+      m_index++;
+      return ret;
+    }
     bool operator==(const iterator &it) const
     {
       return (m_array == it.m_array) && (m_index == it.m_index);
@@ -108,6 +114,12 @@ template <class T> class Array {
       m_index++;
       return *this;
     }
+    const_iterator operator++(int)
+    {
+      auto ret = *this;
+      m_index++;
+      return ret;
+    }
     bool operator==(const const_iterator &it) const
     {
       return (m_array == it.m_array) && (m_index == it.m_index);
diff --git a/src/core/matrix.imp b/src/core/matrix.imp
index 9356c5ef2..04d7625a7 100644
--- a/src/core/matrix.imp
+++ b/src/core/matrix.imp
@@ -169,7 +169,7 @@ template <class T> void Matrix<T>::CMultiply(const Vector<T> &in, Vector<T> &out
     T sum = (T)0;
 
     T *src1 = this->data[i] + this->mincol;
-    T *src2 = in.data + this->mincol;
+    auto src2 = in.begin();
     int j = this->maxcol - this->mincol + 1;
     while (j--) {
       sum += *(src1++) * *(src2++);
@@ -217,7 +217,7 @@ template <class T> void Matrix<T>::RMultiply(const Vector<T> &in, Vector<T> &out
     T k = in[i];
 
     T *src = this->data[i] + this->mincol;
-    T *dst = out.data + this->mincol;
+    auto dst = out.begin();
     int j = this->maxcol - this->mincol + 1;
     while (j--) {
       *(dst++) += *(src++) * k;
diff --git a/src/core/vector.cc b/src/core/vector.cc
deleted file mode 100644
index 412dd5975..000000000
--- a/src/core/vector.cc
+++ /dev/null
@@ -1,30 +0,0 @@
-//
-// This file is part of Gambit
-// Copyright (c) 1994-2024, The Gambit Project (http://www.gambit-project.org)
-//
-// FILE: src/libgambit/vector.cc
-// Instantiation of vector types
-//
-// This program is free software; you can redistribute it and/or modify
-// it under the terms of the GNU General Public License as published by
-// the Free Software Foundation; either version 2 of the License, or
-// (at your option) any later version.
-//
-// This program is distributed in the hope that it will be useful,
-// but WITHOUT ANY WARRANTY; without even the implied warranty of
-// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-// GNU General Public License for more details.
-//
-// You should have received a copy of the GNU General Public License
-// along with this program; if not, write to the Free Software
-// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
-//
-
-#include "gambit.h"
-#include "vector.imp"
-
-template class Gambit::Vector<int>;
-template class Gambit::Vector<long>;
-template class Gambit::Vector<double>;
-template class Gambit::Vector<Gambit::Integer>;
-template class Gambit::Vector<Gambit::Rational>;
diff --git a/src/core/vector.h b/src/core/vector.h
index a2cd0521b..8b6ffb185 100644
--- a/src/core/vector.h
+++ b/src/core/vector.h
@@ -23,13 +23,20 @@
 #ifndef LIBGAMBIT_VECTOR_H
 #define LIBGAMBIT_VECTOR_H
 
+#include <numeric>
+
 namespace Gambit {
 
 template <class T> class Matrix;
 
 /// A mathematical vector: a list of numbers with the standard math operators
 template <class T> class Vector : public Array<T> {
-  friend class Matrix<T>;
+private:
+  // check vector for identical boundaries
+  bool Check(const Vector<T> &v) const
+  {
+    return (v.mindex == this->mindex && v.maxdex == this->maxdex);
+  }
 
 public:
   /** Create a vector of length len, starting at 1 */
@@ -42,35 +49,119 @@ template <class T> class Vector : public Array<T> {
   ~Vector() override = default;
 
   /** Assignment operator: requires vectors to be of same length */
-  Vector<T> &operator=(const Vector<T> &V);
+  Vector<T> &operator=(const Vector<T> &V)
+  {
+    if (!Check(V)) {
+      throw DimensionException();
+    }
+    Array<T>::operator=(V);
+    return *this;
+  }
   /** Assigns the value c to all components of the vector */
-  Vector<T> &operator=(T c);
-
-  Vector<T> operator+(const Vector<T> &V) const;
-  Vector<T> &operator+=(const Vector<T> &V);
-
-  Vector<T> operator-();
-  Vector<T> operator-(const Vector<T> &V) const;
-  Vector<T> &operator-=(const Vector<T> &V);
-
-  Vector<T> operator*(T c) const;
-  Vector<T> &operator*=(T c);
-  T operator*(const Vector<T> &V) const;
-
-  Vector<T> operator/(T c) const;
-
-  bool operator==(const Vector<T> &V) const;
+  Vector<T> &operator=(const T &c)
+  {
+    std::fill(this->begin(), this->end(), c);
+    return *this;
+  }
+
+  Vector<T> operator+(const Vector<T> &V) const
+  {
+    if (!Check(V)) {
+      throw DimensionException();
+    }
+    Vector<T> tmp(this->mindex, this->maxdex);
+    std::transform(this->cbegin(), this->cend(), V.cbegin(), tmp.begin(), std::plus<>());
+    return tmp;
+  }
+
+  Vector<T> &operator+=(const Vector<T> &V)
+  {
+    if (!Check(V)) {
+      throw DimensionException();
+    }
+    std::transform(this->cbegin(), this->cend(), V.cbegin(), this->begin(), std::plus<>());
+    return *this;
+  }
+
+  Vector<T> operator-() const
+  {
+    Vector<T> tmp(this->mindex, this->maxdex);
+    std::transform(tmp.cbegin(), tmp.cend(), tmp.begin(), std::negate<>());
+    return tmp;
+  }
+
+  Vector<T> operator-(const Vector<T> &V) const
+  {
+    if (!Check(V)) {
+      throw DimensionException();
+    }
+    Vector<T> tmp(this->mindex, this->maxdex);
+    std::transform(this->cbegin(), this->cend(), V.cbegin(), tmp.begin(), std::minus<>());
+    return tmp;
+  }
+
+  Vector<T> &operator-=(const Vector<T> &V)
+  {
+    if (!Check(V)) {
+      throw DimensionException();
+    }
+    std::transform(this->cbegin(), this->cend(), V.cbegin(), this->begin(), std::minus<>());
+    return *this;
+  }
+
+  Vector<T> operator*(const T &c) const
+  {
+    Vector<T> tmp(this->mindex, this->maxdex);
+    std::transform(this->cbegin(), this->cend(), tmp.begin(), [&](const T &v) { return v * c; });
+    return tmp;
+  }
+
+  Vector<T> &operator*=(const T &c)
+  {
+    std::transform(this->cbegin(), this->cend(), this->begin(), [&](const T &v) { return v * c; });
+    return *this;
+  }
+
+  T operator*(const Vector<T> &V) const
+  {
+    if (!Check(V)) {
+      throw DimensionException();
+    }
+    return std::inner_product(this->begin(), this->end(), V.begin(), static_cast<T>(0));
+  }
+
+  Vector<T> operator/(const T &c) const
+  {
+    Vector<T> tmp(this->mindex, this->maxdex);
+    std::transform(this->cbegin(), this->cend(), tmp.begin(), [&](const T &v) { return v / c; });
+    return tmp;
+  }
+
+  bool operator==(const Vector<T> &V) const
+  {
+    if (!Check(V)) {
+      throw DimensionException();
+    }
+    return Array<T>::operator==(V);
+  }
   bool operator!=(const Vector<T> &V) const { return !(*this == V); }
 
   /** Tests if all components of the vector are equal to a constant c */
-  bool operator==(T c) const;
-  bool operator!=(T c) const;
+  bool operator==(const T &c) const
+  {
+    return std::all_of(this->begin(), this->end(), [&](const T &v) { return v == c; });
+  }
+  bool operator!=(const T &c) const
+  {
+    return std::any_of(this->begin(), this->end(), [&](const T &v) { return v != c; });
+  }
 
   // square of length
-  T NormSquared() const;
-
-  // check vector for identical boundaries
-  bool Check(const Vector<T> &v) const;
+  T NormSquared() const
+  {
+    return std::accumulate(this->begin(), this->end(), static_cast<T>(0),
+                           [](const T &t, const T &v) { return t + v * v; });
+  }
 };
 
 } // end namespace Gambit
diff --git a/src/core/vector.imp b/src/core/vector.imp
deleted file mode 100644
index 1c0a1cd5b..000000000
--- a/src/core/vector.imp
+++ /dev/null
@@ -1,199 +0,0 @@
-//
-// This file is part of Gambit
-// Copyright (c) 1994-2024, The Gambit Project (http://www.gambit-project.org)
-//
-// FILE: src/libgambit/vector.imp
-// Implementation of vector class
-//
-// This program is free software; you can redistribute it and/or modify
-// it under the terms of the GNU General Public License as published by
-// the Free Software Foundation; either version 2 of the License, or
-// (at your option) any later version.
-//
-// This program is distributed in the hope that it will be useful,
-// but WITHOUT ANY WARRANTY; without even the implied warranty of
-// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-// GNU General Public License for more details.
-//
-// You should have received a copy of the GNU General Public License
-// along with this program; if not, write to the Free Software
-// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
-//
-
-#include "vector.h"
-
-namespace Gambit {
-
-//------------------------------------------------------------------------
-//     Vector<T>: Constructors, destructor, constructive operators
-//------------------------------------------------------------------------
-
-template <class T> Vector<T> &Vector<T>::operator=(const Vector<T> &V)
-{
-  if (!Check(V)) {
-    throw DimensionException();
-  }
-  Array<T>::operator=(V);
-  return *this;
-}
-
-//------------------------------------------------------------------------
-// inline arithmetic operators
-//------------------------------------------------------------------------
-
-template <class T> bool Vector<T>::operator!=(T c) const { return !(*this == c); }
-
-//------------------------------------------------------------------------
-// inline internal functions
-//------------------------------------------------------------------------
-
-template <class T> bool Vector<T>::Check(const Vector<T> &v) const
-{
-  return (v.mindex == this->mindex && v.maxdex == this->maxdex);
-}
-
-//------------------------------------------------------------------------
-//	Vector: arithmetic operators
-//------------------------------------------------------------------------
-
-template <class T> Vector<T> &Vector<T>::operator=(T c)
-{
-  for (int i = this->mindex; i <= this->maxdex; i++) {
-    (*this)[i] = c;
-  }
-  return (*this);
-}
-
-template <class T> bool Vector<T>::operator==(const Vector<T> &V) const
-{
-  if (!Check(V)) {
-    throw DimensionException();
-  }
-
-  for (int i = this->mindex; i <= this->maxdex; i++) {
-    if ((*this)[i] != V[i]) {
-      return false;
-    }
-  }
-  return true;
-}
-
-// arithmetic operators
-template <class T> Vector<T> Vector<T>::operator+(const Vector<T> &V) const
-{
-  if (!Check(V)) {
-    throw DimensionException();
-  }
-
-  Vector<T> tmp(this->mindex, this->maxdex);
-  for (int i = this->mindex; i <= this->maxdex; i++) {
-    tmp[i] = (*this)[i] + V[i];
-  }
-  return tmp;
-}
-
-template <class T> Vector<T> Vector<T>::operator-(const Vector<T> &V) const
-{
-  if (!Check(V)) {
-    throw DimensionException();
-  }
-
-  Vector<T> tmp(this->mindex, this->maxdex);
-  for (int i = this->mindex; i <= this->maxdex; i++) {
-    tmp[i] = (*this)[i] - V[i];
-  }
-  return tmp;
-}
-
-template <class T> Vector<T> &Vector<T>::operator+=(const Vector<T> &V)
-{
-  if (!Check(V)) {
-    throw DimensionException();
-  }
-
-  for (int i = this->mindex; i <= this->maxdex; i++) {
-    (*this)[i] += V[i];
-  }
-  return (*this);
-}
-
-template <class T> Vector<T> &Vector<T>::operator-=(const Vector<T> &V)
-{
-  if (!Check(V)) {
-    throw DimensionException();
-  }
-
-  for (int i = this->mindex; i <= this->maxdex; i++) {
-    (*this)[i] -= V[i];
-  }
-  return (*this);
-}
-
-template <class T> Vector<T> Vector<T>::operator-()
-{
-  Vector<T> tmp(this->mindex, this->maxdex);
-  for (int i = this->mindex; i <= this->maxdex; i++) {
-    tmp[i] = -(*this)[i];
-  }
-  return tmp;
-}
-
-template <class T> Vector<T> Vector<T>::operator*(T c) const
-{
-  Vector<T> tmp(this->mindex, this->maxdex);
-  for (int i = this->mindex; i <= this->maxdex; i++) {
-    tmp[i] = (*this)[i] * c;
-  }
-  return tmp;
-}
-
-template <class T> Vector<T> &Vector<T>::operator*=(T c)
-{
-  for (int i = this->mindex; i <= this->maxdex; i++) {
-    (*this)[i] *= c;
-  }
-  return (*this);
-}
-
-template <class T> T Vector<T>::operator*(const Vector<T> &V) const
-{
-  if (!Check(V)) {
-    throw DimensionException();
-  }
-
-  T sum = (T)0;
-  for (int i = this->mindex; i <= this->maxdex; i++) {
-    sum += (*this)[i] * V[i];
-  }
-  return sum;
-}
-
-template <class T> Vector<T> Vector<T>::operator/(T c) const
-{
-  Vector<T> tmp(this->mindex, this->maxdex);
-  for (int i = this->mindex; i <= this->maxdex; i++) {
-    tmp[i] = (*this)[i] / c;
-  }
-  return tmp;
-}
-
-template <class T> bool Vector<T>::operator==(T c) const
-{
-  for (int i = this->mindex; i <= this->maxdex; i++) {
-    if ((*this)[i] != c) {
-      return false;
-    }
-  }
-  return true;
-}
-
-template <class T> T Vector<T>::NormSquared() const
-{
-  T answer = (T)0;
-  for (int i = 1; i <= this->Length(); i++) {
-    answer += (*this)[i] * (*this)[i];
-  }
-  return answer;
-}
-
-} // end namespace Gambit
diff --git a/src/solvers/enumpoly/gcomplex.cc b/src/solvers/enumpoly/gcomplex.cc
index 063de3fe6..016fda390 100644
--- a/src/solvers/enumpoly/gcomplex.cc
+++ b/src/solvers/enumpoly/gcomplex.cc
@@ -77,7 +77,3 @@ gComplex pow(const gComplex &x, long y)
     return answer;
   }
 }
-
-#include "core/vector.imp"
-
-template class Gambit::Vector<gComplex>;
diff --git a/src/solvers/enumpoly/linrcomb.cc b/src/solvers/enumpoly/linrcomb.cc
index dc625101a..ef56e1719 100644
--- a/src/solvers/enumpoly/linrcomb.cc
+++ b/src/solvers/enumpoly/linrcomb.cc
@@ -21,8 +21,5 @@
 //
 
 #include "linrcomb.imp"
-#include "core/vector.imp"
-#include "core/matrix.imp"
 
 template class LinearCombination<double>;
-// template class LinearCombination<Gambit::Rational>;