diff --git a/mini-odeint/mini-odeint-tests.cpp b/mini-odeint/mini-odeint-tests.cpp
index 3d37ccd..40481d7 100644
--- a/mini-odeint/mini-odeint-tests.cpp
+++ b/mini-odeint/mini-odeint-tests.cpp
@@ -4,8 +4,6 @@
 
 #include "mini-odeint.hpp"
 
-#include <ranges>
-
 unsigned int Factorial(unsigned int number) {
   return number <= 1 ? number : Factorial(number - 1) * number;
 }
diff --git a/mini-odeint/mini-odeint.hpp b/mini-odeint/mini-odeint.hpp
index 91d84d5..a218a6b 100644
--- a/mini-odeint/mini-odeint.hpp
+++ b/mini-odeint/mini-odeint.hpp
@@ -158,29 +158,47 @@ template <typename E> struct Vec2 {
   E x, y;
 
   constexpr Vec2() = default;
-  constexpr Vec2(E x, E y) : x(x), y(y) {}
+  constexpr Vec2(E x, E y) : x(std::move(x)), y(std::move(y)) {}
 
   constexpr explicit Vec2(std::array<E, 3> v) : x(v[0]), y(v[1]) {}
 
   friend constexpr Vec2 operator+(Vec2 lhs, const Vec2 &rhs) {
     return lhs += rhs;
   }
+  friend constexpr Vec2 operator-(Vec2 lhs, const Vec2 &rhs) {
+    return lhs -= rhs;
+  }
   friend constexpr Vec2 operator*(Vec2 lhs, const value_type &rhs) {
     return lhs *= rhs;
   }
   friend constexpr Vec2 operator*(const value_type &lhs, Vec2 rhs) {
     return rhs *= lhs;
   }
+  friend constexpr Vec2 operator/(Vec2 lhs, const value_type &rhs) {
+    return lhs /= rhs;
+  }
+  constexpr Vec2 operator-() const { return Vec2{-x, -y}; }
+
   constexpr Vec2 &operator+=(const Vec2 &rhs) {
     x += rhs.x;
     y += rhs.y;
     return *this;
   }
+  constexpr Vec2 &operator-=(const Vec2 &rhs) {
+    x -= rhs.x;
+    y -= rhs.y;
+    return *this;
+  }
   constexpr Vec2 &operator*=(const value_type &rhs) {
     x *= rhs;
     y *= rhs;
     return *this;
   }
+  constexpr Vec2 &operator/=(const value_type &rhs) {
+    x /= rhs;
+    y /= rhs;
+    return *this;
+  }
 };
 
 template <typename E> struct Vec3 {
@@ -189,31 +207,56 @@ template <typename E> struct Vec3 {
   E x, y, z;
 
   constexpr Vec3() = default;
-  constexpr Vec3(E x, E y, E z) : x(x), y(y), z(z) {}
+  constexpr Vec3(E x, E y, E z)
+      : x(std::move(x)), y(std::move(y)), z(std::move(z)) {}
 
-  constexpr explicit Vec3(std::array<E, 3> v) : x(v[0]), y(v[1]), z(v[2]) {}
+  explicit constexpr Vec3(const std::array<E, 3> &v)
+      : x(v[0]), y(v[1]), z(v[2]) {}
+
+  explicit constexpr Vec3(std::span<const E, 3> v)
+      : x(v[0]), y(v[1]), z(v[2]) {}
 
   friend constexpr Vec3 operator+(Vec3 lhs, const Vec3 &rhs) {
     return lhs += rhs;
   }
+  friend constexpr Vec3 operator-(Vec3 lhs, const Vec3 &rhs) {
+    return lhs -= rhs;
+  }
   friend constexpr Vec3 operator*(Vec3 lhs, const value_type &rhs) {
     return lhs *= rhs;
   }
   friend constexpr Vec3 operator*(const value_type &lhs, Vec3 rhs) {
     return rhs *= lhs;
   }
+  friend constexpr Vec3 operator/(Vec3 lhs, const value_type &rhs) {
+    return lhs /= rhs;
+  }
+  constexpr Vec3 operator-() const { return Vec3{-x, -y, -z}; }
+
   constexpr Vec3 &operator+=(const Vec3 &rhs) {
     x += rhs.x;
     y += rhs.y;
     z += rhs.z;
     return *this;
   }
+  constexpr Vec3 &operator-=(const Vec3 &rhs) {
+    x -= rhs.x;
+    y -= rhs.y;
+    z -= rhs.z;
+    return *this;
+  }
   constexpr Vec3 &operator*=(const value_type &rhs) {
     x *= rhs;
     y *= rhs;
     z *= rhs;
     return *this;
   }
+  constexpr Vec3 &operator/=(const value_type &rhs) {
+    x /= rhs;
+    y /= rhs;
+    z /= rhs;
+    return *this;
+  }
 };
 
 namespace func {