diff --git a/src/include/dim.hxx b/src/include/dim.hxx
index 089097a27..4e4453728 100644
--- a/src/include/dim.hxx
+++ b/src/include/dim.hxx
@@ -22,7 +22,8 @@ struct Invar
     return {!InvarX::value, !InvarY::value, !InvarZ::value};
   }
 
-  static bool is_invar(int dim)
+  // TODO c++20 make consteval?
+  static constexpr bool is_invar(int dim)
   {
     switch (dim) {
       case 0: return InvarX::value;
diff --git a/src/include/rng.hxx b/src/include/rng.hxx
index 9b8c7236e..2842c5710 100644
--- a/src/include/rng.hxx
+++ b/src/include/rng.hxx
@@ -59,6 +59,12 @@ struct Uniform
 
   Real get() { return dist(*gen); }
 
+  Real get(Real min_override, Real max_override)
+  {
+    return std::uniform_real_distribution<Real>(min_override,
+                                                max_override)(*gen);
+  }
+
 private:
   detail::GeneratorPtr gen;
   std::uniform_real_distribution<Real> dist;
@@ -81,11 +87,11 @@ struct Normal
   Normal() : Normal(0, 1) {}
 
   Real get() { return dist(*gen); }
-  // FIXME remove me, or make standalone func
-  Real get(Real mean, Real stdev)
+
+  Real get(Real mean_override, Real stdev_override)
   {
     // should be possible to pass params to existing dist
-    return std::normal_distribution<Real>(mean, stdev)(*gen);
+    return std::normal_distribution<Real>(mean_override, stdev_override)(*gen);
   }
 
 private:
diff --git a/src/kg/include/kg/Vec3.h b/src/kg/include/kg/Vec3.h
index aff5ba926..b77395e3b 100644
--- a/src/kg/include/kg/Vec3.h
+++ b/src/kg/include/kg/Vec3.h
@@ -9,6 +9,7 @@
 #include "cuda_compat.h"
 #include <kg/Macros.h>
 #include <psc/gtensor.h>
+#include "psc_bits.h"
 
 namespace kg
 {
@@ -70,6 +71,14 @@ struct Vec : gt::sarray<T, N>
     return *this;
   }
 
+  KG_INLINE Vec& operator+=(T s)
+  {
+    for (size_t i = 0; i < N; i++) {
+      (*this)[i] += s;
+    }
+    return *this;
+  }
+
   KG_INLINE Vec& operator-=(const Vec& w)
   {
     for (size_t i = 0; i < N; i++) {
@@ -78,6 +87,14 @@ struct Vec : gt::sarray<T, N>
     return *this;
   }
 
+  KG_INLINE Vec& operator-=(T s)
+  {
+    for (size_t i = 0; i < N; i++) {
+      (*this)[i] -= s;
+    }
+    return *this;
+  }
+
   KG_INLINE Vec& operator*=(const Vec& w)
   {
     for (size_t i = 0; i < N; i++) {
@@ -102,6 +119,148 @@ struct Vec : gt::sarray<T, N>
     return *this;
   }
 
+  KG_INLINE Vec& operator/=(T s)
+  {
+    for (size_t i = 0; i < N; i++) {
+      (*this)[i] /= s;
+    }
+    return *this;
+  }
+
+  /**
+   * @brief Calculates the sum of this vec's elements.
+   * @return the sum
+   */
+  KG_INLINE T sum() const
+  {
+    T sum{0};
+    for (size_t i = 0; i < N; i++) {
+      sum += (*this)[i];
+    }
+    return sum;
+  }
+
+  /**
+   * @brief Calculates the product of this vec's elements.
+   * @return the product
+   */
+  KG_INLINE T prod() const
+  {
+    T prod{1};
+    for (size_t i = 0; i < N; i++) {
+      prod *= (*this)[i];
+    }
+    return prod;
+  }
+
+  /**
+   * @brief Calculates the dot product of this and another vec.
+   * @param w the other vec
+   * @return the dot product
+   */
+  KG_INLINE T dot(const Vec& w) const
+  {
+    T sum{0};
+    for (size_t i = 0; i < N; i++) {
+      sum += (*this)[i] * w[i];
+    }
+    return sum;
+  }
+
+  /**
+   * @brief Calculates the cross product of this and another vec. This method
+   * assumes `N` is 3.
+   * @param w the right-hand-side vec
+   * @return the cross product
+   */
+  KG_INLINE Vec cross(const Vec& w) const
+  {
+    return {(*this)[1] * w[2] - (*this)[2] * w[1],
+            (*this)[2] * w[0] - (*this)[0] * w[2],
+            (*this)[0] * w[1] - (*this)[1] * w[0]};
+  }
+
+  /**
+   * @brief Calculates the square of the magnitude of this vec.
+   * @return the squared magnitude
+   */
+  KG_INLINE T mag2() const { return this->dot(*this); }
+
+  /**
+   * @brief Calculates the magnitude of this vec by taking the square root of
+   * the squared magnitude.
+   * @return the magnitude
+   */
+  KG_INLINE T mag() const { return sqrt(this->mag2()); }
+
+  /**
+   * @brief Determines the maximal value of any of this vec's elements.
+   * @return the maximal value
+   */
+  KG_INLINE T max() const
+  {
+    return *std::max_element(this->begin(), this->end());
+  }
+
+  /**
+   * @brief Determines the minimal value of any of this vec's elements.
+   * @return the minimal value
+   */
+  KG_INLINE T min() const
+  {
+    return *std::min_element(this->begin(), this->end());
+  }
+
+  /**
+   * @brief Calculates the elementwise maximum between two vecs.
+   * @param v the first vec
+   * @param w the second vec
+   * @return a vec containing the maximal values
+   */
+  static KG_INLINE Vec max(const Vec& v, const Vec& w)
+  {
+    Vec res;
+    for (int i = 0; i < N; i++) {
+      res[i] = std::max(v[i], w[i]);
+    }
+    return res;
+  }
+
+  /**
+   * @brief Calculates the elementwise minimum between two vecs.
+   * @param v the first vec
+   * @param w the second vec
+   * @return a vec containing the minimal values
+   */
+  static KG_INLINE Vec min(const Vec& v, const Vec& w)
+  {
+    Vec res;
+    for (int i = 0; i < N; i++) {
+      res[i] = std::min(v[i], w[i]);
+    }
+    return res;
+  }
+
+  /**
+   * @brief Calculates the (multiplicative) inverse of each of this vec's
+   * elements.
+   * @return a vec of the inverses
+   */
+  KG_INLINE Vec inv() const { return T(1) / *this; }
+
+  /**
+   * @brief Calculates the floor of each of this vec's elements.
+   * @return an integer vec of the floors
+   */
+  KG_INLINE Vec<int, N> fint() const
+  {
+    Vec<int, N> res;
+    for (size_t i = 0; i < N; i++) {
+      res[i] = ::fint((*this)[i]);
+    }
+    return res;
+  }
+
   // conversion to pointer
 
   KG_INLINE operator const T*() const { return this->data(); }
@@ -129,6 +288,14 @@ KG_INLINE Vec<T, N> operator+(const Vec<T, N>& v, const Vec<T, N>& w)
   return res;
 }
 
+template <typename T, std::size_t N>
+KG_INLINE Vec<T, N> operator+(T s, const Vec<T, N>& v)
+{
+  Vec<T, N> res = v;
+  res += s;
+  return res;
+}
+
 template <typename T, std::size_t N>
 KG_INLINE Vec<T, N> operator-(const Vec<T, N>& v, const Vec<T, N>& w)
 {
@@ -137,6 +304,23 @@ KG_INLINE Vec<T, N> operator-(const Vec<T, N>& v, const Vec<T, N>& w)
   return res;
 }
 
+template <typename T, std::size_t N>
+KG_INLINE Vec<T, N> operator-(T s, const Vec<T, N>& v)
+{
+  Vec<T, N> res;
+  for (size_t i = 0; i < N; i++) {
+    res[i] = s - v[i];
+  }
+  return res;
+}
+template <typename T, std::size_t N>
+KG_INLINE Vec<T, N> operator-(const Vec<T, N>& v, T s)
+{
+  Vec<T, N> res = v;
+  res -= s;
+  return res;
+}
+
 template <typename T, std::size_t N>
 KG_INLINE Vec<T, N> operator*(const Vec<T, N>& v, const Vec<T, N>& w)
 {
@@ -153,6 +337,24 @@ KG_INLINE Vec<T, N> operator*(T s, const Vec<T, N>& v)
   return res;
 }
 
+template <typename T, std::size_t N>
+KG_INLINE Vec<T, N> operator/(T s, const Vec<T, N>& v)
+{
+  Vec<T, N> res;
+  for (size_t i = 0; i < N; i++) {
+    res[i] = s / v[i];
+  }
+  return res;
+}
+
+template <typename T, std::size_t N>
+KG_INLINE Vec<T, N> operator/(const Vec<T, N>& v, T s)
+{
+  Vec<T, N> res = v;
+  res /= s;
+  return res;
+}
+
 template <typename T, std::size_t N>
 KG_INLINE Vec<T, N> operator/(const Vec<T, N>& v, const Vec<T, N>& w)
 {
diff --git a/src/kg/include/kg/VecRange.hxx b/src/kg/include/kg/VecRange.hxx
new file mode 100644
index 000000000..e8590ea55
--- /dev/null
+++ b/src/kg/include/kg/VecRange.hxx
@@ -0,0 +1,80 @@
+#pragma once
+
+#include "Vec3.h"
+
+/**
+ * @brief Facilitates iteration over `N`-dimensional spaces in row-major order.
+ *
+ * For example,
+ * ```c++
+ * for (Int3 i3 : VecRange(mins, maxs)) {
+ *     // ...
+ * }
+ * ```
+ * is equivalent to
+ * ```c++
+ * for (int ix = mins[0]; ix < maxs[0]; ix++) {
+ *     for (int iy = mins[1]; iy < maxs[1]; iy++) {
+ *         for (int iz = mins[2]; iz < maxs[2]; iz++) {
+ *             Int3 i3{ix, iy, iz};
+ *             // ...
+ *         }
+ *     }
+ * }
+ * ```
+ * @tparam T scalar type (e.g. `int`)
+ * @tparam N number of dimensions
+ */
+template <typename T, std::size_t N>
+class VecRange
+{
+public:
+  using Vec = kg::Vec<T, N>;
+
+  class Iterator
+  {
+  public:
+    Iterator(const Vec& starts, const Vec& stops, bool end = false)
+      : starts_(starts), stops_(stops), current_(starts), end_(end)
+    {}
+
+    Iterator& operator++()
+    {
+      for (int d = N - 1; d >= 0; d--) {
+        ++current_[d];
+
+        if (current_[d] < stops_[d]) {
+          break;
+        } else {
+          current_[d] = starts_[d];
+
+          if (d == 0) {
+            end_ = true;
+          }
+        }
+      }
+      return *this;
+    }
+
+    Vec operator*() { return current_; }
+
+    bool operator!=(const Iterator& other) const { return end_ != other.end_; }
+
+  private:
+    const Vec& starts_;
+    const Vec& stops_;
+    Vec current_;
+    bool end_;
+  };
+
+  VecRange(const Vec& starts, const Vec& stops) : starts_(starts), stops_(stops)
+  {}
+
+  Iterator begin() const { return Iterator(starts_, stops_, false); }
+
+  Iterator end() const { return Iterator(starts_, stops_, true); }
+
+private:
+  const Vec starts_;
+  const Vec stops_;
+};
\ No newline at end of file
diff --git a/src/libpsc/psc_push_particles/push_particles_1vb.hxx b/src/libpsc/psc_push_particles/push_particles_1vb.hxx
index 912ad48ad..18c8ff180 100644
--- a/src/libpsc/psc_push_particles/push_particles_1vb.hxx
+++ b/src/libpsc/psc_push_particles/push_particles_1vb.hxx
@@ -24,7 +24,7 @@ struct PushParticlesVb
   static void push_mprts(Mparticles& mprts, MfieldsState& mflds)
   {
     const auto& grid = mprts.grid();
-    Real3 dxi = Real3{1., 1., 1.} / Real3(grid.domain.dx);
+    Real3 dxi = Real3(grid.domain.dx).inv();
     real_t dq_kind[MAX_NR_KINDS];
     auto& kinds = grid.kinds;
     assert(kinds.size() <= MAX_NR_KINDS);
@@ -60,10 +60,8 @@ struct PushParticlesVb
         auto v = advance.calc_v(prt.u());
         advance.push_x(prt.x(), v);
 
-        Int3 final_index;
-        Real3 final_pos_normalized;
-        find_idx_pos_1st_rel(prt.x(), dxi, final_index, final_pos_normalized,
-                             real_t(0.));
+        Real3 final_pos_normalized = prt.x() * dxi;
+        Int3 final_index = final_pos_normalized.fint();
 
         // CURRENT DENSITY BETWEEN (n+.5)*dt and (n+1.5)*dt
         int lg[3];