diff --git a/src/include/lattice.h b/src/include/lattice.h
index 6b571a3..4544132 100644
--- a/src/include/lattice.h
+++ b/src/include/lattice.h
@@ -13,7 +13,7 @@
 
 namespace pivot {
 
-template <int Dim> struct box;
+template <int Dim, bool Simd> struct box;
 
 /**
  * @brief Represents a closed interval [left_, right_].
@@ -39,7 +39,7 @@ struct interval {
 /**
  * @brief Represents a Dim-dimensional point.
  */
-template <int Dim> class point : boost::multipliable<point<Dim>, int> {
+template <int Dim, bool Simd = false> class point : boost::multipliable<point<Dim, Simd>, int> {
 
 public:
   point() = default;
@@ -77,10 +77,10 @@ template <int Dim> class point : boost::multipliable<point<Dim>, int> {
   std::array<int, Dim> coords_{};
 };
 
-template <typename S, typename T, T Dim> point(std::array<S, Dim>) -> point<Dim>;
+template <typename S, typename T, T Dim> point(std::array<S, Dim>) -> point<Dim, false>;
 
 /** @brief Represents a Dim-dimensional box. */
-template <int Dim> struct box : boost::additive<box<Dim>, point<Dim>> {
+template <int Dim, bool Simd = false> struct box : boost::additive<box<Dim, Simd>, point<Dim, Simd>> {
   std::array<interval, Dim> intervals_;
 
   box() = delete;
@@ -94,7 +94,7 @@ template <int Dim> struct box : boost::additive<box<Dim>, point<Dim>> {
    * in the sense that the input sequences of points (p0, ...) is effectively translated by e0 - p0
    * prior to the box's construction. In other words, the box will always have a vertex at e0.
    */
-  box(std::span<const point<Dim>> points);
+  box(std::span<const point<Dim, Simd>> points);
 
   bool operator==(const box &b) const;
 
@@ -105,9 +105,9 @@ template <int Dim> struct box : boost::additive<box<Dim>, point<Dim>> {
   bool empty() const;
 
   /** @brief Action of a point (understood as a vector) on a box */
-  box<Dim> &operator+=(const point<Dim> &b);
+  box<Dim, Simd> &operator+=(const point<Dim, Simd> &b);
 
-  box<Dim> &operator-=(const point<Dim> &b);
+  box<Dim, Simd> &operator-=(const point<Dim, Simd> &b);
 
   /**
    * @brief Returns the "union" of two boxes.
@@ -129,7 +129,7 @@ template <int Dim> struct box : boost::additive<box<Dim>, point<Dim>> {
   std::string to_string() const;
 };
 
-template <typename S, typename T, S Dim, T N> box(std::array<point<Dim>, N>) -> box<Dim>;
+template <typename S, typename T, S Dim, T N, bool Simd> box(std::array<point<Dim, Simd>, N>) -> box<Dim, Simd>;
 
 struct point_hash {
   int num_steps_;
@@ -138,7 +138,7 @@ struct point_hash {
 
   // This hashing method, which exploits the known range of values that can be taken by the
   // sequence of points in a walk, ppears to result in better performance than other methods tested.
-  template <int Dim> std::size_t operator()(const point<Dim> &p) const;
+  template <int Dim, bool Simd = false> std::size_t operator()(const point<Dim, Simd> &p) const;
 };
 
 /**
@@ -156,7 +156,7 @@ struct point_hash {
  * entries of S by S(i). Note that P e(i) = e(P(i)) and S e(i) = S(i) e(i), where e(i)
  * is the i-th standard unit vector.
  */
-template <int Dim> class transform {
+template <int Dim, bool Simd = false> class transform {
 
 public:
   /** @brief Constructs the identity transformation. */
@@ -174,7 +174,7 @@ template <int Dim> class transform {
    *
    * Note that the resulting pivot transformation is not uniquely defined by the inputs.
    */
-  transform(const point<Dim> &p, const point<Dim> &q);
+  transform(const point<Dim, Simd> &p, const point<Dim, Simd> &q);
 
   /** @brief Produce a uniformly random transfom.*/
   template <typename Gen> static transform rand(Gen &gen) {
@@ -203,7 +203,7 @@ template <int Dim> class transform {
    *
    * which means the P(i)-th component of the result is S(P(i)) p[i].
    */
-  point<Dim> operator*(const point<Dim> &p) const;
+  point<Dim, Simd> operator*(const point<Dim, Simd> &p) const;
 
   /**
    * @brief Composes two transforms.
@@ -238,7 +238,7 @@ template <int Dim> class transform {
    * which means the P(i)-th interval from which SP(B) is constructed has bounds S(P(i)) a[i] and
    * S(P(i)) b[i].
    */
-  box<Dim> operator*(const box<Dim> &b) const;
+  box<Dim, Simd> operator*(const box<Dim, Simd> &b) const;
 
   /** @brief Returns true if the transform is the identity. */
   bool is_identity() const;
diff --git a/src/include/utils.h b/src/include/utils.h
index 1ea524a..b5e6099 100644
--- a/src/include/utils.h
+++ b/src/include/utils.h
@@ -4,10 +4,10 @@
 
 namespace pivot {
 
-template <int Dim> std::vector<point<Dim>> from_csv(const std::string &path);
+template <int Dim, bool Simd> std::vector<point<Dim, Simd>> from_csv(const std::string &path);
 
-template <int Dim> void to_csv(const std::string &path, const std::vector<point<Dim>> &points);
+template <int Dim, bool Simd> void to_csv(const std::string &path, const std::vector<point<Dim, Simd>> &points);
 
-template <int Dim> std::vector<point<Dim>> line(int num_steps);
+template <int Dim, bool Simd> std::vector<point<Dim, Simd>> line(int num_steps);
 
 } // namespace pivot
diff --git a/src/include/walk.h b/src/include/walk.h
index 9be21d7..636cf31 100644
--- a/src/include/walk.h
+++ b/src/include/walk.h
@@ -11,10 +11,10 @@
 
 namespace pivot {
 
-template <int Dim> class walk : public walk_base<Dim> {
+template <int Dim, bool Simd = false> class walk : public walk_base<Dim, Simd> {
 
 public:
-  walk(const std::vector<point<Dim>> &steps, std::optional<unsigned int> seed = std::nullopt);
+  walk(const std::vector<point<Dim, Simd>> &steps, std::optional<unsigned int> seed = std::nullopt);
   walk(int num_steps, std::optional<unsigned int> seed = std::nullopt);
   walk(const std::string &path, std::optional<unsigned int> seed = std::nullopt);
 
@@ -24,15 +24,15 @@ template <int Dim> class walk : public walk_base<Dim> {
 
   ~walk() = default;
 
-  point<Dim> operator[](int i) const { return steps_[i]; }
+  point<Dim, Simd> operator[](int i) const { return steps_[i]; }
 
   int num_steps() const { return steps_.size(); }
 
-  point<Dim> endpoint() const override { return steps_.back(); }
+  point<Dim, Simd> endpoint() const override { return steps_.back(); }
 
-  std::optional<std::vector<point<Dim>>> try_pivot(int step, const transform<Dim> &trans) const;
+  std::optional<std::vector<point<Dim, Simd>>> try_pivot(int step, const transform<Dim, Simd> &trans) const;
 
-  std::pair<int, std::optional<std::vector<point<Dim>>>> try_rand_pivot() const;
+  std::pair<int, std::optional<std::vector<point<Dim, Simd>>>> try_rand_pivot() const;
 
   bool rand_pivot(bool fast = false) override;
 
@@ -43,15 +43,15 @@ template <int Dim> class walk : public walk_base<Dim> {
   void export_csv(const std::string &path) const override;
 
 protected:
-  std::vector<point<Dim>> steps_;
-  boost::unordered_flat_map<point<Dim>, int, point_hash> occupied_;
+  std::vector<point<Dim, Simd>> steps_;
+  boost::unordered_flat_map<point<Dim, Simd>, int, point_hash> occupied_;
 
   mutable std::mt19937 rng_;
   mutable std::uniform_int_distribution<int> dist_;
 
-  void do_pivot(int step, std::vector<point<Dim>> &new_points);
+  void do_pivot(int step, std::vector<point<Dim, Simd>> &new_points);
 
-  point<Dim> pivot_point(int step, int i, const transform<Dim> &trans) const;
+  point<Dim, Simd> pivot_point(int step, int i, const transform<Dim, Simd> &trans) const;
 };
 
 } // namespace pivot
diff --git a/src/include/walk_base.h b/src/include/walk_base.h
index fdb0e87..36ee01e 100644
--- a/src/include/walk_base.h
+++ b/src/include/walk_base.h
@@ -4,7 +4,7 @@
 
 namespace pivot {
 
-template <int Dim> class walk_base {
+template <int Dim, bool Simd = false> class walk_base {
 
 public:
   virtual ~walk_base() = default;
@@ -15,7 +15,7 @@ template <int Dim> class walk_base {
 
   virtual void export_csv(const std::string &path) const = 0;
 
-  virtual point<Dim> endpoint() const = 0;
+  virtual point<Dim, Simd> endpoint() const = 0;
 };
 
 } // namespace pivot
diff --git a/src/include/walk_node.h b/src/include/walk_node.h
index 038c7c2..2f4ec31 100644
--- a/src/include/walk_node.h
+++ b/src/include/walk_node.h
@@ -10,18 +10,19 @@ namespace pivot {
 
 /* FORWARD REFERENCES */
 
-template <int Dim> class walk_node;
+template <int Dim, bool Simd> class walk_node;
 
-template <int Dim>
-bool intersect(const walk_node<Dim> *l_walk, const walk_node<Dim> *r_walk, const point<Dim> &l_anchor,
-               const point<Dim> &r_anchor, const transform<Dim> &l_symm, const transform<Dim> &r_symm);
+template <int Dim, bool Simd = false>
+bool intersect(const walk_node<Dim, Simd> *l_walk, const walk_node<Dim, Simd> *r_walk, const point<Dim, Simd> &l_anchor,
+               const point<Dim, Simd> &r_anchor, const transform<Dim, Simd> &l_symm,
+               const transform<Dim, Simd> &r_symm);
 
-template <int Dim> class walk_tree;
+template <int Dim, bool Simd> class walk_tree;
 
 /* WALK NODE */
 
 /** @brief Represents a node in a walk tree. */
-template <int Dim> class walk_node {
+template <int Dim, bool Simd = false> class walk_node {
 
 public:
   /* CONSTRUCTORS, DESTRUCTOR */
@@ -33,7 +34,7 @@ template <int Dim> class walk_node {
    *
    * @return The root of the walk tree.
    */
-  static walk_node *pivot_rep(const std::vector<point<Dim>> &steps, walk_node *buf = nullptr);
+  static walk_node *pivot_rep(const std::vector<point<Dim, Simd>> &steps, walk_node *buf = nullptr);
 
   /** @brief Returns the root of a walk tree for the balanced representation of a walk given by a sequence of points.
    *
@@ -42,7 +43,7 @@ template <int Dim> class walk_node {
    *
    * @return The root of the walk tree.
    */
-  static walk_node *balanced_rep(const std::vector<point<Dim>> &steps, walk_node *buf = nullptr);
+  static walk_node *balanced_rep(const std::vector<point<Dim, Simd>> &steps, walk_node *buf = nullptr);
 
   /** @brief Copies the given node but none of the nodes it links to. */
   walk_node(const walk_node &w) = default;
@@ -57,11 +58,11 @@ template <int Dim> class walk_node {
 
   int id() const { return id_; }
 
-  const box<Dim> &bbox() const { return bbox_; }
+  const box<Dim, Simd> &bbox() const { return bbox_; }
 
-  const point<Dim> &endpoint() const { return end_; }
+  const point<Dim, Simd> &endpoint() const { return end_; }
 
-  const transform<Dim> &symm() const { return symm_; }
+  const transform<Dim, Simd> &symm() const { return symm_; }
 
   walk_node *left() const { return left_; }
 
@@ -125,7 +126,7 @@ template <int Dim> class walk_node {
    *
    * @return Whether the transform creates an intersection.
    */
-  bool shuffle_intersect(const transform<Dim> &t, std::optional<bool> is_left_child);
+  bool shuffle_intersect(const transform<Dim, Simd> &t, std::optional<bool> is_left_child);
 
   /**
    * @brief Checks if the current walk has an intersection via a top-down algorithm.
@@ -136,7 +137,7 @@ template <int Dim> class walk_node {
 
   /* OTHER FUNCTIONS */
 
-  std::vector<point<Dim>> steps() const;
+  std::vector<point<Dim, Simd>> steps() const;
 
   void todot(const std::string &path) const;
 
@@ -146,15 +147,16 @@ template <int Dim> class walk_node {
   walk_node *parent_{};
   walk_node *left_{};
   walk_node *right_{};
-  transform<Dim> symm_;
-  box<Dim> bbox_;
-  point<Dim> end_;
+  transform<Dim, Simd> symm_;
+  box<Dim, Simd> bbox_;
+  point<Dim, Simd> end_;
 
-  friend class walk_tree<Dim>;
+  friend class walk_tree<Dim, Simd>;
 
   /* CONVENIENCE METHODS */
 
-  walk_node(int id, int num_sites, const transform<Dim> &symm, const box<Dim> &bbox, const point<Dim> &end);
+  walk_node(int id, int num_sites, const transform<Dim, Simd> &symm, const box<Dim, Simd> &bbox,
+            const point<Dim, Simd> &end);
 
   void set_left(walk_node *left) {
     left_ = left;
@@ -185,15 +187,15 @@ template <int Dim> class walk_node {
   Agnode_t *todot(Agraph_t *g, const cgraph_t &cgraph) const;
 
   // recursive helper
-  static walk_node *balanced_rep(std::span<const point<Dim>> steps, int start, const transform<Dim> &glob_symm,
-                                 walk_node *buf);
+  static walk_node *balanced_rep(std::span<const point<Dim, Simd>> steps, int start,
+                                 const transform<Dim, Simd> &glob_symm, walk_node *buf);
 
-  bool shuffle_intersect(const transform<Dim> &t, std::optional<bool> was_left_child,
+  bool shuffle_intersect(const transform<Dim, Simd> &t, std::optional<bool> was_left_child,
                          std::optional<bool> is_left_child);
 
-  template <int D>
-  friend bool intersect(const walk_node<D> *l_walk, const walk_node<D> *r_walk, const point<D> &l_anchor,
-                        const point<D> &r_anchor, const transform<D> &l_symm, const transform<D> &r_symm);
+  template <int D, bool S>
+  friend bool intersect(const walk_node<D, S> *l_walk, const walk_node<D, S> *r_walk, const point<D, S> &l_anchor,
+                        const point<D, S> &r_anchor, const transform<D, S> &l_symm, const transform<D, S> &r_symm);
 };
 
 } // namespace pivot
diff --git a/src/include/walk_tree.h b/src/include/walk_tree.h
index 1afb8ed..63af7c5 100644
--- a/src/include/walk_tree.h
+++ b/src/include/walk_tree.h
@@ -10,10 +10,10 @@
 
 namespace pivot {
 
-template <int Dim> class walk_node;
+template <int Dim, bool Simd> class walk_node;
 
 /** @brief Represents an entire saw-tree (as per Clisby's 2010 paper). */
-template <int Dim> class walk_tree : public walk_base<Dim> {
+template <int Dim, bool Simd = false> class walk_tree : public walk_base<Dim, Simd> {
 
 public:
   /* CONSTRUCTORS, DESTRUCTOR */
@@ -56,7 +56,7 @@ template <int Dim> class walk_tree : public walk_base<Dim> {
    *
    * @warning It is not recommended to set balanced=false.
    */
-  walk_tree(const std::vector<point<Dim>> &steps, std::optional<unsigned int> seed = std::nullopt,
+  walk_tree(const std::vector<point<Dim, Simd>> &steps, std::optional<unsigned int> seed = std::nullopt,
             bool balanced = true);
 
   /**@brief Deallocates the entire tree and every node it contains. */
@@ -64,9 +64,9 @@ template <int Dim> class walk_tree : public walk_base<Dim> {
 
   /* GETTERS, SETTERS, SIMPLE UTILITIES */
 
-  walk_node<Dim> *root() const;
+  walk_node<Dim, Simd> *root() const;
 
-  point<Dim> endpoint() const override;
+  point<Dim, Simd> endpoint() const override;
 
   bool is_leaf() const;
 
@@ -84,7 +84,7 @@ template <int Dim> class walk_tree : public walk_base<Dim> {
    *
    * @return reference to the node
    */
-  walk_node<Dim> &find_node(int n);
+  walk_node<Dim, Simd> &find_node(int n);
 
   /* HIGH-LEVEL FUNCTIONS (see Clisby (2010), Section 2.7) */
 
@@ -99,7 +99,7 @@ template <int Dim> class walk_tree : public walk_base<Dim> {
    *
    * @return Whether the pivot was successful.
    */
-  bool try_pivot(int n, const transform<Dim> &r);
+  bool try_pivot(int n, const transform<Dim, Simd> &r);
 
   /**
    * @brief Attempt to pivot the walk about the given lattice site using Clisby's Attempt_pivot_fast function.
@@ -112,7 +112,7 @@ template <int Dim> class walk_tree : public walk_base<Dim> {
    *
    * @return Whether the pivot was successful.
    */
-  bool try_pivot_fast(int n, const transform<Dim> &r);
+  bool try_pivot_fast(int n, const transform<Dim, Simd> &r);
 
   /**
    * @brief Attempts to pivot the walk about a randomly chosen lattice site with a random transform.
@@ -130,7 +130,7 @@ template <int Dim> class walk_tree : public walk_base<Dim> {
    *
    * @return Sequence of lattice sites.
    */
-  std::vector<point<Dim>> steps() const;
+  std::vector<point<Dim, Simd>> steps() const;
 
   /**
    * @brief Check whether the walk is self-avoiding using a naive algorithm.
@@ -149,10 +149,10 @@ template <int Dim> class walk_tree : public walk_base<Dim> {
   void todot(const std::string &path) const;
 
 private:
-  std::unique_ptr<walk_node<Dim>> root_;
+  std::unique_ptr<walk_node<Dim, Simd>> root_;
   std::mt19937 rng_;
   std::uniform_int_distribution<int> dist_; // distribution for choosing a random lattice site
-  walk_node<Dim> *buf_;                     // buffer into which nodes are allocated (used for fast node lookup by id)
+  walk_node<Dim, Simd> *buf_;               // buffer into which nodes are allocated (used for fast node lookup by id)
 };
 
 } // namespace pivot
diff --git a/src/lattice/box.hpp b/src/lattice/box.hpp
index 2b93179..bf6de71 100644
--- a/src/lattice/box.hpp
+++ b/src/lattice/box.hpp
@@ -24,7 +24,7 @@ std::string interval::to_string() const {
   return result;
 }
 
-template <int Dim> box<Dim> &box<Dim>::operator+=(const point<Dim> &p) {
+template <int Dim, bool Simd> box<Dim, Simd> &box<Dim, Simd>::operator+=(const point<Dim, Simd> &p) {
   for (int i = 0; i < Dim; ++i) {
     intervals_[i].left_ += p[i];
     intervals_[i].right_ += p[i];
@@ -32,7 +32,7 @@ template <int Dim> box<Dim> &box<Dim>::operator+=(const point<Dim> &p) {
   return *this;
 }
 
-template <int Dim> box<Dim> &box<Dim>::operator-=(const point<Dim> &b) {
+template <int Dim, bool Simd> box<Dim, Simd> &box<Dim, Simd>::operator-=(const point<Dim, Simd> &b) {
   for (int i = 0; i < Dim; ++i) {
     intervals_[i].left_ -= b[i];
     intervals_[i].right_ -= b[i];
@@ -40,9 +40,9 @@ template <int Dim> box<Dim> &box<Dim>::operator-=(const point<Dim> &b) {
   return *this;
 }
 
-template <int Dim> box<Dim>::box(const std::array<interval, Dim> &intervals) : intervals_(intervals) {}
+template <int Dim, bool Simd> box<Dim, Simd>::box(const std::array<interval, Dim> &intervals) : intervals_(intervals) {}
 
-template <int Dim> box<Dim>::box(std::span<const point<Dim>> points) {
+template <int Dim, bool Simd> box<Dim, Simd>::box(std::span<const point<Dim, Simd>> points) {
   std::array<int, Dim> min;
   std::array<int, Dim> max;
   min.fill(std::numeric_limits<int>::max());
@@ -60,17 +60,17 @@ template <int Dim> box<Dim>::box(std::span<const point<Dim>> points) {
   }
 }
 
-template <int Dim> bool box<Dim>::operator==(const box &b) const { return intervals_ == b.intervals_; }
+template <int Dim, bool Simd> bool box<Dim, Simd>::operator==(const box &b) const { return intervals_ == b.intervals_; }
 
-template <int Dim> bool box<Dim>::operator!=(const box &b) const { return intervals_ != b.intervals_; }
+template <int Dim, bool Simd> bool box<Dim, Simd>::operator!=(const box &b) const { return intervals_ != b.intervals_; }
 
-template <int Dim> interval box<Dim>::operator[](int i) const { return intervals_[i]; }
+template <int Dim, bool Simd> interval box<Dim, Simd>::operator[](int i) const { return intervals_[i]; }
 
-template <int Dim> bool box<Dim>::empty() const {
+template <int Dim, bool Simd> bool box<Dim, Simd>::empty() const {
   return std::any_of(intervals_.begin(), intervals_.end(), [](const interval &i) { return i.empty(); });
 }
 
-template <int Dim> box<Dim> box<Dim>::operator|(const box<Dim> &b) const {
+template <int Dim, bool Simd> box<Dim, Simd> box<Dim, Simd>::operator|(const box<Dim, Simd> &b) const {
   std::array<interval, Dim> intervals;
   for (int i = 0; i < Dim; ++i) {
     intervals[i] = interval(std::min(intervals_[i].left_, b.intervals_[i].left_),
@@ -79,7 +79,7 @@ template <int Dim> box<Dim> box<Dim>::operator|(const box<Dim> &b) const {
   return box(intervals);
 }
 
-template <int Dim> box<Dim> box<Dim>::operator&(const box<Dim> &b) const {
+template <int Dim, bool Simd> box<Dim, Simd> box<Dim, Simd>::operator&(const box<Dim, Simd> &b) const {
   std::array<interval, Dim> intervals;
   for (int i = 0; i < Dim; ++i) {
     intervals[i] = interval(std::max(intervals_[i].left_, b.intervals_[i].left_),
@@ -88,7 +88,7 @@ template <int Dim> box<Dim> box<Dim>::operator&(const box<Dim> &b) const {
   return box(intervals);
 }
 
-template <int Dim> std::string box<Dim>::to_string() const {
+template <int Dim, bool Simd> std::string box<Dim, Simd>::to_string() const {
   std::string s = "";
   for (int i = 0; i < Dim - 1; ++i) {
     s += intervals_[i].to_string() + " x ";
diff --git a/src/lattice/point.hpp b/src/lattice/point.hpp
index 8264ed1..ddc44ce 100644
--- a/src/lattice/point.hpp
+++ b/src/lattice/point.hpp
@@ -4,17 +4,17 @@
 
 namespace pivot {
 
-template <int Dim> point<Dim>::point(const std::array<int, Dim> &coords) : coords_(coords) {}
+template <int Dim, bool Simd> point<Dim, Simd>::point(const std::array<int, Dim> &coords) : coords_(coords) {}
 
-template <int Dim> point<Dim> point<Dim>::unit(int i) {
-  point<Dim> result;
+template <int Dim, bool Simd> point<Dim, Simd> point<Dim, Simd>::unit(int i) {
+  point<Dim, Simd> result;
   result.coords_[i] = 1;
   return result;
 }
 
-template <int Dim> int point<Dim>::operator[](int i) const { return coords_[i]; }
+template <int Dim, bool Simd> int point<Dim, Simd>::operator[](int i) const { return coords_[i]; }
 
-template <int Dim> bool point<Dim>::operator==(const point &p) const {
+template <int Dim, bool Simd> bool point<Dim, Simd>::operator==(const point &p) const {
   // Faster than using std::equal, at least for small arrays
   // https://stackoverflow.com/questions/39262496/why-is-stdequal-much-slower-than-a-hand-rolled-loop-for-two-small-stdarray
   for (int i = 0; i < Dim; ++i) {
@@ -25,9 +25,9 @@ template <int Dim> bool point<Dim>::operator==(const point &p) const {
   return true;
 }
 
-template <int Dim> bool point<Dim>::operator!=(const point &p) const { return !(*this == p); }
+template <int Dim, bool Simd> bool point<Dim, Simd>::operator!=(const point &p) const { return !(*this == p); }
 
-template <int Dim> point<Dim> point<Dim>::operator+(const point<Dim> &p) const {
+template <int Dim, bool Simd> point<Dim, Simd> point<Dim, Simd>::operator+(const point<Dim, Simd> &p) const {
   std::array<int, Dim> sum;
   for (int i = 0; i < Dim; ++i) {
     sum[i] = coords_[i] + p.coords_[i];
@@ -35,7 +35,7 @@ template <int Dim> point<Dim> point<Dim>::operator+(const point<Dim> &p) const {
   return point(sum);
 }
 
-template <int Dim> point<Dim> point<Dim>::operator-(const point &p) const {
+template <int Dim, bool Simd> point<Dim, Simd> point<Dim, Simd>::operator-(const point &p) const {
   std::array<int, Dim> diff;
   for (int i = 0; i < Dim; ++i) {
     diff[i] = coords_[i] - p.coords_[i];
@@ -43,14 +43,14 @@ template <int Dim> point<Dim> point<Dim>::operator-(const point &p) const {
   return point(diff);
 }
 
-template <int Dim> point<Dim> &point<Dim>::operator*=(int k) {
+template <int Dim, bool Simd> point<Dim, Simd> &point<Dim, Simd>::operator*=(int k) {
   for (int i = 0; i < Dim; ++i) {
     coords_[i] *= k;
   }
   return *this;
 }
 
-template <int Dim> int point<Dim>::norm() const {
+template <int Dim, bool Simd> int point<Dim, Simd>::norm() const {
   int sum = 0;
   for (int i = 0; i < Dim; ++i) {
     sum += coords_[i] * coords_[i];
@@ -58,7 +58,7 @@ template <int Dim> int point<Dim>::norm() const {
   return sum;
 }
 
-template <int Dim> std::string point<Dim>::to_string() const {
+template <int Dim, bool Simd> std::string point<Dim, Simd>::to_string() const {
   std::string s = "(";
   for (int i = 0; i < Dim - 1; ++i) {
     s += std::to_string(coords_[i]) + ", ";
@@ -74,7 +74,7 @@ template <int Dim> std::string point<Dim>::to_string() const {
 
 point_hash::point_hash(int num_steps) : num_steps_(num_steps) {}
 
-template <int Dim> std::size_t point_hash::operator()(const point<Dim> &p) const {
+template <int Dim, bool Simd> std::size_t point_hash::operator()(const point<Dim, Simd> &p) const {
   std::size_t hash = 0;
   for (int i = 0; i < Dim; ++i) {
     hash = num_steps_ + p[i] + 2 * num_steps_ * hash;
diff --git a/src/lattice/transform.hpp b/src/lattice/transform.hpp
index 7b035c0..1fe7fde 100644
--- a/src/lattice/transform.hpp
+++ b/src/lattice/transform.hpp
@@ -4,20 +4,21 @@
 
 namespace pivot {
 
-template <int Dim> transform<Dim>::transform() {
+template <int Dim, bool Simd> transform<Dim, Simd>::transform() {
   for (int i = 0; i < Dim; ++i) {
     perm_[i] = i;  // trivial permutation
     signs_[i] = 1; // standard orientations (no flips)
   }
 }
 
-template <int Dim>
-transform<Dim>::transform(const std::array<int, Dim> &perm, const std::array<int, Dim> &signs)
+template <int Dim, bool Simd>
+transform<Dim, Simd>::transform(const std::array<int, Dim> &perm, const std::array<int, Dim> &signs)
     : perm_(perm), signs_(signs) {}
 
-template <int Dim> transform<Dim>::transform(const point<Dim> &p, const point<Dim> &q) : transform() {
+template <int Dim, bool Simd>
+transform<Dim, Simd>::transform(const point<Dim, Simd> &p, const point<Dim, Simd> &q) : transform() {
   // The input points should differ by 1 in a single coordinate. Start by finding this coordinate or fail.
-  point<Dim> diff = q - p;
+  point<Dim, Simd> diff = q - p;
   int idx = -1;
   for (int i = 0; i < Dim; ++i) {
     if (std::abs(diff[i]) == 1) {
@@ -39,25 +40,26 @@ template <int Dim> transform<Dim>::transform(const point<Dim> &p, const point<Di
   signs_[idx] = diff[idx];
 }
 
-template <int Dim> transform<Dim> transform<Dim>::rand() {
+template <int Dim, bool Simd> transform<Dim, Simd> transform<Dim, Simd>::rand() {
   static std::random_device rd;
   static std::mt19937 gen(rd());
   return rand(gen);
 }
 
-template <int Dim> bool transform<Dim>::operator==(const transform &t) const {
+template <int Dim, bool Simd> bool transform<Dim, Simd>::operator==(const transform &t) const {
   return perm_ == t.perm_ && signs_ == t.signs_;
 }
 
-template <int Dim> point<Dim> transform<Dim>::operator*(const point<Dim> &p) const {
+template <int Dim, bool Simd> point<Dim, Simd> transform<Dim, Simd>::operator*(const point<Dim, Simd> &p) const {
   std::array<int, Dim> coords;
   for (int i = 0; i < Dim; ++i) {
     coords[perm_[i]] = signs_[perm_[i]] * p[i];
   }
-  return point<Dim>(coords);
+  return point<Dim, Simd>(coords);
 }
 
-template <int Dim> transform<Dim> transform<Dim>::operator*(const transform<Dim> &t) const {
+template <int Dim, bool Simd>
+transform<Dim, Simd> transform<Dim, Simd>::operator*(const transform<Dim, Simd> &t) const {
   std::array<int, Dim> perm;
   std::array<int, Dim> signs;
   for (int i = 0; i < Dim; ++i) {
@@ -67,17 +69,17 @@ template <int Dim> transform<Dim> transform<Dim>::operator*(const transform<Dim>
   return transform(perm, signs);
 }
 
-template <int Dim> box<Dim> transform<Dim>::operator*(const box<Dim> &b) const {
+template <int Dim, bool Simd> box<Dim, Simd> transform<Dim, Simd>::operator*(const box<Dim, Simd> &b) const {
   std::array<interval, Dim> intervals;
   for (int i = 0; i < Dim; ++i) {
     int x = signs_[perm_[i]] * b.intervals_[i].left_;
     int y = signs_[perm_[i]] * b.intervals_[i].right_;
     intervals[perm_[i]] = interval(std::min(x, y), std::max(x, y));
   }
-  return box<Dim>(intervals);
+  return box<Dim, Simd>(intervals);
 }
 
-template <int Dim> bool transform<Dim>::is_identity() const {
+template <int Dim, bool Simd> bool transform<Dim, Simd>::is_identity() const {
   for (int i = 0; i < Dim; ++i) {
     if (perm_[i] != i || signs_[i] != 1) {
       return false;
@@ -86,7 +88,7 @@ template <int Dim> bool transform<Dim>::is_identity() const {
   return true;
 }
 
-template <int Dim> transform<Dim> transform<Dim>::inverse() const {
+template <int Dim, bool Simd> transform<Dim, Simd> transform<Dim, Simd>::inverse() const {
   std::array<int, Dim> perm;
   std::array<int, Dim> signs;
   for (int i = 0; i < Dim; ++i) {
@@ -96,7 +98,7 @@ template <int Dim> transform<Dim> transform<Dim>::inverse() const {
   return transform(perm, signs);
 }
 
-template <int Dim> std::array<std::array<int, Dim>, Dim> transform<Dim>::to_matrix() const {
+template <int Dim, bool Simd> std::array<std::array<int, Dim>, Dim> transform<Dim, Simd>::to_matrix() const {
   std::array<std::array<int, Dim>, Dim> matrix = {};
   for (int i = 0; i < Dim; ++i) {
     matrix[perm_[i]][i] = signs_[perm_[i]];
@@ -104,7 +106,7 @@ template <int Dim> std::array<std::array<int, Dim>, Dim> transform<Dim>::to_matr
   return matrix;
 }
 
-template <int Dim> std::string transform<Dim>::to_string() const {
+template <int Dim, bool Simd> std::string transform<Dim, Simd>::to_string() const {
   auto matrix = to_matrix();
   std::string s = "[";
   for (int i = 0; i < Dim; ++i) {
diff --git a/src/loop.h b/src/loop.h
index 8987a79..e37fe3d 100644
--- a/src/loop.h
+++ b/src/loop.h
@@ -6,15 +6,15 @@
 #include "walk.h"
 #include "walk_tree.h"
 
-template <int Dim>
+template <int Dim, bool Simd = false>
 int main_loop(int num_steps, int iters, bool naive, bool fast, int seed, bool require_success, bool verify,
               const std::string &in_path, const std::string &out_dir) {
   std::unique_ptr<pivot::walk_base<Dim>> w;
   if (naive) {
     if (in_path.empty()) {
-      w = std::make_unique<pivot::walk<Dim>>(num_steps, seed);
+      w = std::make_unique<pivot::walk<Dim, Simd>>(num_steps, seed);
     } else {
-      w = std::make_unique<pivot::walk<Dim>>(in_path, seed);
+      w = std::make_unique<pivot::walk<Dim, Simd>>(in_path, seed);
     }
   } else {
     if (in_path.empty()) {
diff --git a/src/utils/utils.cpp b/src/utils/utils.cpp
index 05975f1..cb82d0a 100644
--- a/src/utils/utils.cpp
+++ b/src/utils/utils.cpp
@@ -8,9 +8,9 @@
 
 namespace pivot {
 
-template <int Dim> std::vector<point<Dim>> from_csv(const std::string &path) {
+template <int Dim, bool Simd = false> std::vector<point<Dim, Simd>> from_csv(const std::string &path) {
   std::ifstream file(path);
-  std::vector<point<Dim>> points;
+  std::vector<point<Dim, Simd>> points;
   std::string line;
   while (std::getline(file, line)) {
     std::array<int, Dim> coords;
@@ -23,12 +23,13 @@ template <int Dim> std::vector<point<Dim>> from_csv(const std::string &path) {
       coords[i] = std::stoi(line.substr(start, end - start));
       start = end + 1;
     }
-    points.push_back(point<Dim>(coords));
+    points.push_back(point<Dim, Simd>(coords));
   }
   return points;
 }
 
-template <int Dim> void to_csv(const std::string &path, const std::vector<point<Dim>> &points) {
+template <int Dim, bool Simd = false>
+void to_csv(const std::string &path, const std::vector<point<Dim, Simd>> &points) {
   // TODO: check path exists
   std::ofstream file(path);
   for (const auto &p : points) {
@@ -39,17 +40,18 @@ template <int Dim> void to_csv(const std::string &path, const std::vector<point<
   }
 }
 
-template <int Dim> std::vector<point<Dim>> line(int num_steps) {
-  std::vector<point<Dim>> steps(num_steps);
+template <int Dim, bool Simd = false> std::vector<point<Dim, Simd>> line(int num_steps) {
+  std::vector<point<Dim, Simd>> steps(num_steps);
   for (int i = 0; i < num_steps; ++i) {
-    steps[i] = (i + 1) * point<Dim>::unit(0);
+    steps[i] = (i + 1) * point<Dim, Simd>::unit(0);
   }
   return steps;
 }
 
-#define FROM_CSV_INST(z, n, data) template std::vector<point<n>> from_csv<n>(const std::string &path);
-#define TO_CSV_INST(z, n, data) template void to_csv<n>(const std::string &path, const std::vector<point<n>> &points);
-#define LINE_INST(z, n, data) template std::vector<point<n>> line(int num_steps);
+#define FROM_CSV_INST(z, n, data) template std::vector<point<n>> from_csv<n, false>(const std::string &path);
+#define TO_CSV_INST(z, n, data)                                                                                        \
+  template void to_csv<n, false>(const std::string &path, const std::vector<point<n>> &points);
+#define LINE_INST(z, n, data) template std::vector<point<n>> line<n, false>(int num_steps);
 
 // cppcheck-suppress syntaxError
 BOOST_PP_REPEAT_FROM_TO(1, DIMS_UB, TO_CSV_INST, ~)
diff --git a/src/walks/node/ctors.hpp b/src/walks/node/ctors.hpp
index fd71c4d..8c98b44 100644
--- a/src/walks/node/ctors.hpp
+++ b/src/walks/node/ctors.hpp
@@ -2,27 +2,29 @@
 
 namespace pivot {
 
-template <int Dim>
-walk_node<Dim>::walk_node(int id, int num_sites, const transform<Dim> &symm, const box<Dim> &bbox,
-                          const point<Dim> &end)
+template <int Dim, bool Simd>
+walk_node<Dim, Simd>::walk_node(int id, int num_sites, const transform<Dim, Simd> &symm, const box<Dim, Simd> &bbox,
+                                const point<Dim, Simd> &end)
     : id_(id), num_sites_(num_sites), symm_(symm), bbox_(bbox), end_(end) {}
 
-template <int Dim> walk_node<Dim> *walk_node<Dim>::pivot_rep(const std::vector<point<Dim>> &steps, walk_node *buf) {
+template <int Dim, bool Simd>
+walk_node<Dim, Simd> *walk_node<Dim, Simd>::pivot_rep(const std::vector<point<Dim, Simd>> &steps, walk_node *buf) {
   int num_sites = steps.size();
   if (num_sites < 2) {
     throw std::invalid_argument("num_sites must be at least 2");
   }
-  walk_node<Dim> *root =
-      buf ? new (buf) walk_node(1, num_sites, transform(steps[0], steps[1]), box<Dim>(steps), steps[num_sites - 1])
-          : new walk_node(1, num_sites, transform(steps[0], steps[1]), box<Dim>(steps), steps[num_sites - 1]);
+  walk_node<Dim, Simd> *root =
+      buf ? new (buf)
+                walk_node(1, num_sites, transform(steps[0], steps[1]), box<Dim, Simd>(steps), steps[num_sites - 1])
+          : new walk_node(1, num_sites, transform(steps[0], steps[1]), box<Dim, Simd>(steps), steps[num_sites - 1]);
   auto node = root;
   for (int i = 0; i < num_sites - 2; ++i) {
     auto id = i + 2;
-    node->right_ = buf ? new (buf + id - 1)
-                             walk_node(id, num_sites - i - 1, transform(steps[i + 1], steps[i + 2]),
-                                       box(std::span<const point<Dim>>(steps).subspan(i + 1)), steps[num_sites - 1])
+    node->right_ = buf ? new (buf + id - 1) walk_node(id, num_sites - i - 1, transform(steps[i + 1], steps[i + 2]),
+                                                      box(std::span<const point<Dim, Simd>>(steps).subspan(i + 1)),
+                                                      steps[num_sites - 1])
                        : new walk_node(i + 2, num_sites - i - 1, transform(steps[i + 1], steps[i + 2]),
-                                       box(std::span<const point<Dim>>(steps).subspan(i + 1)),
+                                       box(std::span<const point<Dim, Simd>>(steps).subspan(i + 1)),
                                        steps[num_sites - 1]); // TODO: double-check this
     node->right_->parent_ = node;
     node = node->right_;
@@ -30,14 +32,16 @@ template <int Dim> walk_node<Dim> *walk_node<Dim>::pivot_rep(const std::vector<p
   return root;
 }
 
-template <int Dim>
-walk_node<Dim> *walk_node<Dim>::balanced_rep(const std::vector<point<Dim>> &steps, walk_node<Dim> *buf) {
-  return balanced_rep(steps, 1, transform<Dim>(), buf);
+template <int Dim, bool Simd>
+walk_node<Dim, Simd> *walk_node<Dim, Simd>::balanced_rep(const std::vector<point<Dim, Simd>> &steps,
+                                                         walk_node<Dim, Simd> *buf) {
+  return balanced_rep(steps, 1, transform<Dim, Simd>(), buf);
 }
 
-template <int Dim>
-walk_node<Dim> *walk_node<Dim>::balanced_rep(std::span<const point<Dim>> steps, int start,
-                                             const transform<Dim> &glob_symm, walk_node<Dim> *buf) {
+template <int Dim, bool Simd>
+walk_node<Dim, Simd> *walk_node<Dim, Simd>::balanced_rep(std::span<const point<Dim, Simd>> steps, int start,
+                                                         const transform<Dim, Simd> &glob_symm,
+                                                         walk_node<Dim, Simd> *buf) {
   int num_sites = steps.size();
   if (num_sites < 1) {
     throw std::invalid_argument("num_sites must be at least 1");
@@ -54,8 +58,8 @@ walk_node<Dim> *walk_node<Dim>::balanced_rep(std::span<const point<Dim>> steps,
   auto abs_symm = transform(steps[n - 1], steps[n]);
   auto glob_inv = glob_symm.inverse();
   auto rel_symm = glob_inv * abs_symm;
-  auto rel_end = glob_inv * (steps.back() - steps.front()) + point<Dim>::unit(0);
-  auto rel_box = point<Dim>::unit(0) + glob_inv * (box(steps) - point<Dim>::unit(0));
+  auto rel_end = glob_inv * (steps.back() - steps.front()) + point<Dim, Simd>::unit(0);
+  auto rel_box = point<Dim, Simd>::unit(0) + glob_inv * (box(steps) - point<Dim, Simd>::unit(0));
   int id = start + n - 1;
   walk_node *root = buf ? new (buf + id - 1) walk_node(id, num_sites, rel_symm, rel_box, rel_end)
                         : new walk_node(id, num_sites, rel_symm, rel_box, rel_end);
@@ -71,20 +75,20 @@ walk_node<Dim> *walk_node<Dim>::balanced_rep(std::span<const point<Dim>> steps,
   return root;
 }
 
-template <int Dim> walk_node<Dim> walk_node<Dim>::create_leaf() {
+template <int Dim, bool Simd> walk_node<Dim, Simd> walk_node<Dim, Simd>::create_leaf() {
   std::array<interval, Dim> intervals;
   intervals[0] = interval(1, 1);
   for (int i = 1; i < Dim; ++i) {
     intervals[i] = interval(0, 0);
   }
-  return walk_node(0, 1, transform<Dim>(), box<Dim>(intervals), point<Dim>::unit(0));
+  return walk_node(0, 1, transform<Dim, Simd>(), box<Dim, Simd>(intervals), point<Dim, Simd>::unit(0));
 }
 
-template <int Dim> walk_node<Dim> &walk_node<Dim>::leaf() {
-  static walk_node<Dim> leaf = create_leaf();
+template <int Dim, bool Simd> walk_node<Dim, Simd> &walk_node<Dim, Simd>::leaf() {
+  static walk_node<Dim, Simd> leaf = create_leaf();
   return leaf;
 }
 
-template <int Dim> walk_node<Dim>::~walk_node() = default;
+template <int Dim, bool Simd> walk_node<Dim, Simd>::~walk_node() = default;
 
 } // namespace pivot
diff --git a/src/walks/node/graphviz.hpp b/src/walks/node/graphviz.hpp
index 1102304..189b225 100644
--- a/src/walks/node/graphviz.hpp
+++ b/src/walks/node/graphviz.hpp
@@ -2,7 +2,7 @@
 
 namespace pivot {
 
-template <int Dim> void walk_node<Dim>::todot(const std::string &path) const {
+template <int Dim, bool Simd> void walk_node<Dim, Simd>::todot(const std::string &path) const {
   gvc_t &gvc = gvc_t::load();
   cgraph_t &cgraph = cgraph_t::load();
 
@@ -20,7 +20,7 @@ template <int Dim> void walk_node<Dim>::todot(const std::string &path) const {
   gvc.gvFreeContext(context);
 }
 
-template <int Dim> Agnode_t *walk_node<Dim>::todot(Agraph_t *g, const cgraph_t &cgraph) const {
+template <int Dim, bool Simd> Agnode_t *walk_node<Dim, Simd>::todot(Agraph_t *g, const cgraph_t &cgraph) const {
   auto name = std::to_string(id_);
   Agnode_t *node = cgraph.agnode(g, (char *)name.c_str(), 1);
   if (cgraph.agset((void *)node, (char *)"shape", (char *)"box") != 0) {
diff --git a/src/walks/node/pivot.hpp b/src/walks/node/pivot.hpp
index 6465da2..35103d3 100644
--- a/src/walks/node/pivot.hpp
+++ b/src/walks/node/pivot.hpp
@@ -7,7 +7,7 @@ namespace pivot {
 // Note: A detail missing from Clisby's paper regarding tree rotations is that parent pointers must be updated,
 // except when the rotation is called from shuffle_intersect.
 
-template <int Dim> walk_node<Dim> *walk_node<Dim>::rotate_left(bool set_parent) {
+template <int Dim, bool Simd> walk_node<Dim, Simd> *walk_node<Dim, Simd>::rotate_left(bool set_parent) {
   if (right_->is_leaf()) {
     throw std::invalid_argument("can't rotate left on a leaf node");
   }
@@ -41,7 +41,7 @@ template <int Dim> walk_node<Dim> *walk_node<Dim>::rotate_left(bool set_parent)
   return this;
 }
 
-template <int Dim> walk_node<Dim> *walk_node<Dim>::rotate_right(bool set_parent) {
+template <int Dim, bool Simd> walk_node<Dim, Simd> *walk_node<Dim, Simd>::rotate_right(bool set_parent) {
   if (left_->is_leaf()) {
     throw std::invalid_argument("can't rotate right on a leaf node");
   }
@@ -75,7 +75,7 @@ template <int Dim> walk_node<Dim> *walk_node<Dim>::rotate_right(bool set_parent)
   return this;
 }
 
-template <int Dim> void walk_node<Dim>::merge() {
+template <int Dim, bool Simd> void walk_node<Dim, Simd>::merge() {
   num_sites_ = left_->num_sites_ + right_->num_sites_;
 
   bbox_ = left_->bbox_ | (left_->end_ + symm_ * right_->bbox_);
@@ -84,7 +84,7 @@ template <int Dim> void walk_node<Dim>::merge() {
 
 /* USER LEVEL OPERATIONS */
 
-template <int Dim> walk_node<Dim> *walk_node<Dim>::shuffle_up(int id) {
+template <int Dim, bool Simd> walk_node<Dim, Simd> *walk_node<Dim, Simd>::shuffle_up(int id) {
   if (id < left_->num_sites_) {
     left_->shuffle_up(id);
     rotate_right();
@@ -96,7 +96,7 @@ template <int Dim> walk_node<Dim> *walk_node<Dim>::shuffle_up(int id) {
   return this;
 }
 
-template <int Dim> walk_node<Dim> *walk_node<Dim>::shuffle_down() {
+template <int Dim, bool Simd> walk_node<Dim, Simd> *walk_node<Dim, Simd>::shuffle_down() {
   int id = std::floor((num_sites_ + 1) / 2.0);
   if (id < left_->num_sites_) {
     rotate_right();
@@ -109,13 +109,14 @@ template <int Dim> walk_node<Dim> *walk_node<Dim>::shuffle_down() {
   return this;
 }
 
-template <int Dim> bool walk_node<Dim>::intersect() const {
-  return ::pivot::intersect<Dim>(left_, right_, point<Dim>(), left_->end_, transform<Dim>(), symm_);
+template <int Dim, bool Simd> bool walk_node<Dim, Simd>::intersect() const {
+  return ::pivot::intersect<Dim, Simd>(left_, right_, point<Dim, Simd>(), left_->end_, transform<Dim, Simd>(), symm_);
 }
 
-template <int Dim>
-bool intersect(const walk_node<Dim> *l_walk, const walk_node<Dim> *r_walk, const point<Dim> &l_anchor,
-               const point<Dim> &r_anchor, const transform<Dim> &l_symm, const transform<Dim> &r_symm) {
+template <int Dim, bool Simd>
+bool intersect(const walk_node<Dim, Simd> *l_walk, const walk_node<Dim, Simd> *r_walk, const point<Dim, Simd> &l_anchor,
+               const point<Dim, Simd> &r_anchor, const transform<Dim, Simd> &l_symm,
+               const transform<Dim, Simd> &r_symm) {
   auto l_box = l_anchor + l_symm * l_walk->bbox_;
   auto r_box = r_anchor + r_symm * r_walk->bbox_;
   if ((l_box & r_box).empty()) {
@@ -137,27 +138,29 @@ bool intersect(const walk_node<Dim> *l_walk, const walk_node<Dim> *r_walk, const
   }
 }
 
-template <int Dim> bool walk_node<Dim>::shuffle_intersect(const transform<Dim> &t, std::optional<bool> is_left_child) {
+template <int Dim, bool Simd>
+bool walk_node<Dim, Simd>::shuffle_intersect(const transform<Dim, Simd> &t, std::optional<bool> is_left_child) {
   return shuffle_intersect(t, std::nullopt, is_left_child);
 }
 
-template <int Dim>
-bool walk_node<Dim>::shuffle_intersect(const transform<Dim> &t, std::optional<bool> was_left_child,
-                                       std::optional<bool> is_left_child) {
+template <int Dim, bool Simd>
+bool walk_node<Dim, Simd>::shuffle_intersect(const transform<Dim, Simd> &t, std::optional<bool> was_left_child,
+                                             std::optional<bool> is_left_child) {
   /* BASE CASE */
   if (was_left_child.has_value()) {
     if (was_left_child.value()) {
-      if (::pivot::intersect(left_, right_->right_, point<Dim>(), left_->end_ + symm_ * t * right_->left_->end_,
-                             transform<Dim>(), symm_ * t * right_->symm_)) {
+      if (::pivot::intersect(left_, right_->right_, point<Dim, Simd>(), left_->end_ + symm_ * t * right_->left_->end_,
+                             transform<Dim, Simd>(), symm_ * t * right_->symm_)) {
         return true;
       }
     } else {
-      if (::pivot::intersect(left_->left_, right_, point<Dim>(), left_->end_, transform<Dim>(), symm_ * t)) {
+      if (::pivot::intersect(left_->left_, right_, point<Dim, Simd>(), left_->end_, transform<Dim, Simd>(),
+                             symm_ * t)) {
         return true;
       }
     }
   } else {
-    if (::pivot::intersect(left_, right_, point<Dim>(), left_->end_, transform<Dim>(), symm_ * t)) {
+    if (::pivot::intersect(left_, right_, point<Dim, Simd>(), left_->end_, transform<Dim, Simd>(), symm_ * t)) {
       return true;
     }
   }
diff --git a/src/walks/node/walk_node.cpp b/src/walks/node/walk_node.cpp
index c301d3f..fcf3ac6 100644
--- a/src/walks/node/walk_node.cpp
+++ b/src/walks/node/walk_node.cpp
@@ -10,7 +10,7 @@
 
 namespace pivot {
 
-template <int Dim> bool walk_node<Dim>::operator==(const walk_node &other) const {
+template <int Dim, bool Simd> bool walk_node<Dim, Simd>::operator==(const walk_node &other) const {
   if (is_leaf() && other.is_leaf()) {
     return true;
   }
@@ -18,10 +18,12 @@ template <int Dim> bool walk_node<Dim>::operator==(const walk_node &other) const
          end_ == other.end_ && *left_ == *other.left_ && *right_ == *other.right_;
 }
 
-template <int Dim> bool walk_node<Dim>::is_leaf() const { return left_ == nullptr && right_ == nullptr; }
+template <int Dim, bool Simd> bool walk_node<Dim, Simd>::is_leaf() const {
+  return left_ == nullptr && right_ == nullptr;
+}
 
-template <int Dim> std::vector<point<Dim>> walk_node<Dim>::steps() const {
-  std::vector<point<Dim>> result;
+template <int Dim, bool Simd> std::vector<point<Dim, Simd>> walk_node<Dim, Simd>::steps() const {
+  std::vector<point<Dim, Simd>> result;
   if (is_leaf()) {
     result.push_back(end_);
     return result;
diff --git a/src/walks/walk.cpp b/src/walks/walk.cpp
index 52600a7..0f4a184 100644
--- a/src/walks/walk.cpp
+++ b/src/walks/walk.cpp
@@ -6,8 +6,8 @@
 
 namespace pivot {
 
-template <int Dim>
-walk<Dim>::walk(const std::vector<point<Dim>> &steps, std::optional<unsigned int> seed)
+template <int Dim, bool Simd>
+walk<Dim, Simd>::walk(const std::vector<point<Dim, Simd>> &steps, std::optional<unsigned int> seed)
     : steps_(steps), occupied_(steps.size(), point_hash(steps.size())) {
   for (int i = 0; i < num_steps(); ++i) {
     occupied_[steps_[i]] = i;
@@ -17,19 +17,21 @@ walk<Dim>::walk(const std::vector<point<Dim>> &steps, std::optional<unsigned int
   dist_ = std::uniform_int_distribution<int>(0, num_steps() - 1);
 }
 
-template <int Dim>
-walk<Dim>::walk(int num_steps, std::optional<unsigned int> seed) : walk(line<Dim>(num_steps), seed) {}
+template <int Dim, bool Simd>
+walk<Dim, Simd>::walk(int num_steps, std::optional<unsigned int> seed) : walk(line<Dim, false>(num_steps), seed) {}
 
-template <int Dim>
-walk<Dim>::walk(const std::string &path, std::optional<unsigned int> seed) : walk(from_csv<Dim>(path), seed) {}
+template <int Dim, bool Simd>
+walk<Dim, Simd>::walk(const std::string &path, std::optional<unsigned int> seed)
+    : walk(from_csv<Dim, false>(path), seed) {}
 
-template <int Dim>
-std::optional<std::vector<point<Dim>>> walk<Dim>::try_pivot(int step, const transform<Dim> &trans) const {
+template <int Dim, bool Simd>
+std::optional<std::vector<point<Dim, Simd>>> walk<Dim, Simd>::try_pivot(int step,
+                                                                        const transform<Dim, Simd> &trans) const {
   if (trans.is_identity()) {
     return {};
   }
 
-  std::vector<point<Dim>> new_points(num_steps() - step - 1);
+  std::vector<point<Dim, Simd>> new_points(num_steps() - step - 1);
   for (int i = step + 1; i < num_steps(); ++i) {
     auto q = pivot_point(step, i, trans);
     auto it = occupied_.find(q);
@@ -41,13 +43,14 @@ std::optional<std::vector<point<Dim>>> walk<Dim>::try_pivot(int step, const tran
   return new_points;
 }
 
-template <int Dim> std::pair<int, std::optional<std::vector<point<Dim>>>> walk<Dim>::try_rand_pivot() const {
+template <int Dim, bool Simd>
+std::pair<int, std::optional<std::vector<point<Dim, Simd>>>> walk<Dim, Simd>::try_rand_pivot() const {
   auto step = dist_(rng_);
-  auto r = transform<Dim>::rand(rng_);
+  auto r = transform<Dim, Simd>::rand(rng_);
   return {step, try_pivot(step, r)};
 }
 
-template <int Dim> bool walk<Dim>::rand_pivot(bool fast) {
+template <int Dim, bool Simd> bool walk<Dim, Simd>::rand_pivot(bool fast) {
   if (fast) { // TODO: implement Kennedy algorithm
     throw std::invalid_argument("fast pivot not implemented for naive walk");
   }
@@ -60,14 +63,14 @@ template <int Dim> bool walk<Dim>::rand_pivot(bool fast) {
   return true;
 }
 
-template <int Dim> bool walk<Dim>::rand_pivot(int num_workers) {
+template <int Dim, bool Simd> bool walk<Dim, Simd>::rand_pivot(int num_workers) {
   if (num_workers == 0) {
     return rand_pivot();
   }
 
   std::vector<int> steps(num_workers);
-  std::vector<std::optional<std::vector<point<Dim>>>> proposals(num_workers);
-  std::vector<std::future<std::pair<int, std::optional<std::vector<point<Dim>>>>>> futures(num_workers);
+  std::vector<std::optional<std::vector<point<Dim, Simd>>>> proposals(num_workers);
+  std::vector<std::future<std::pair<int, std::optional<std::vector<point<Dim, Simd>>>>>> futures(num_workers);
   for (int i = 0; i < num_workers; ++i) {
     futures[i] = std::async(&walk::try_rand_pivot, this);
   }
@@ -89,7 +92,7 @@ template <int Dim> bool walk<Dim>::rand_pivot(int num_workers) {
   return success;
 }
 
-template <int Dim> bool walk<Dim>::self_avoiding() const {
+template <int Dim, bool Simd> bool walk<Dim, Simd>::self_avoiding() const {
   for (int i = 0; i < num_steps(); ++i) {
     for (int j = i + 1; j < num_steps(); ++j) {
       if (steps_[i] == steps_[j]) {
@@ -100,9 +103,11 @@ template <int Dim> bool walk<Dim>::self_avoiding() const {
   return true;
 }
 
-template <int Dim> void walk<Dim>::export_csv(const std::string &path) const { return to_csv(path, steps_); }
+template <int Dim, bool Simd> void walk<Dim, Simd>::export_csv(const std::string &path) const {
+  return to_csv(path, steps_);
+}
 
-template <int Dim> void walk<Dim>::do_pivot(int step, std::vector<point<Dim>> &new_points) {
+template <int Dim, bool Simd> void walk<Dim, Simd>::do_pivot(int step, std::vector<point<Dim, Simd>> &new_points) {
   for (auto it = steps_.begin() + step + 1; it != steps_.end(); ++it) {
     occupied_.erase(*it);
   }
@@ -112,12 +117,13 @@ template <int Dim> void walk<Dim>::do_pivot(int step, std::vector<point<Dim>> &n
   }
 }
 
-template <int Dim> point<Dim> walk<Dim>::pivot_point(int step, int i, const transform<Dim> &trans) const {
+template <int Dim, bool Simd>
+point<Dim, Simd> walk<Dim, Simd>::pivot_point(int step, int i, const transform<Dim, Simd> &trans) const {
   auto p = steps_[step];
   return p + trans * (steps_[i] - p);
 }
 
-#define WALK_INST(z, n, data) template class walk<n>;
+#define WALK_INST(z, n, data) template class walk<n, false>;
 
 // cppcheck-suppress syntaxError
 BOOST_PP_REPEAT_FROM_TO(1, DIMS_UB, WALK_INST, ~)
diff --git a/src/walks/walk_tree.cpp b/src/walks/walk_tree.cpp
index 44a12c8..39380b6 100644
--- a/src/walks/walk_tree.cpp
+++ b/src/walks/walk_tree.cpp
@@ -13,34 +13,35 @@ namespace pivot {
 
 /* CONSTRUCTORS, DESTRUCTOR */
 
-template <int Dim>
-walk_tree<Dim>::walk_tree(int num_sites, std::optional<unsigned int> seed, bool balanced)
-    : walk_tree(line<Dim>(num_sites), seed, balanced) {}
+template <int Dim, bool Simd>
+walk_tree<Dim, Simd>::walk_tree(int num_sites, std::optional<unsigned int> seed, bool balanced)
+    : walk_tree(line<Dim, Simd>(num_sites), seed, balanced) {}
 
-template <int Dim>
-walk_tree<Dim>::walk_tree(const std::string &path, std::optional<unsigned int> seed, bool balanced)
-    : walk_tree(from_csv<Dim>(path), seed, balanced) {}
+template <int Dim, bool Simd>
+walk_tree<Dim, Simd>::walk_tree(const std::string &path, std::optional<unsigned int> seed, bool balanced)
+    : walk_tree(from_csv<Dim, Simd>(path), seed, balanced) {}
 
-template <int Dim>
-walk_tree<Dim>::walk_tree(const std::vector<point<Dim>> &steps, std::optional<unsigned int> seed, bool balanced) {
+template <int Dim, bool Simd>
+walk_tree<Dim, Simd>::walk_tree(const std::vector<point<Dim, Simd>> &steps, std::optional<unsigned int> seed,
+                                bool balanced) {
   if (steps.size() < 2) {
     throw std::invalid_argument("walk must have at least 2 sites (1 step)");
   }
   buf_ = nullptr;
   if (balanced) {
-    size_t buf_size = sizeof(walk_node<Dim>) * (steps.size() - 1);
-    constexpr auto alignment = std::align_val_t(alignof(walk_node<Dim>));
-    buf_ = static_cast<walk_node<Dim> *>(::operator new[](buf_size, alignment));
+    size_t buf_size = sizeof(walk_node<Dim, Simd>) * (steps.size() - 1);
+    constexpr auto alignment = std::align_val_t(alignof(walk_node<Dim, Simd>));
+    buf_ = static_cast<walk_node<Dim, Simd> *>(::operator new[](buf_size, alignment));
   }
-  root_ = balanced ? std::unique_ptr<walk_node<Dim>>(walk_node<Dim>::balanced_rep(steps, buf_))
-                   : std::unique_ptr<walk_node<Dim>>(walk_node<Dim>::pivot_rep(steps, buf_));
+  root_ = balanced ? std::unique_ptr<walk_node<Dim, Simd>>(walk_node<Dim, Simd>::balanced_rep(steps, buf_))
+                   : std::unique_ptr<walk_node<Dim, Simd>>(walk_node<Dim, Simd>::pivot_rep(steps, buf_));
 
   rng_ = std::mt19937(seed.value_or(std::random_device()()));
   dist_ = std::uniform_int_distribution<int>(1, steps.size() - 1);
 }
 
-template <int Dim> walk_tree<Dim>::~walk_tree() {
-  std::stack<walk_node<Dim> *> nodes;
+template <int Dim, bool Simd> walk_tree<Dim, Simd>::~walk_tree() {
+  std::stack<walk_node<Dim, Simd> *> nodes;
   nodes.push(root_.release());
   while (!nodes.empty()) {
     auto curr = nodes.top();
@@ -55,32 +56,32 @@ template <int Dim> walk_tree<Dim>::~walk_tree() {
   }
 
   if (buf_) {
-    ::operator delete[](buf_, std::align_val_t(alignof(walk_node<Dim>)));
+    ::operator delete[](buf_, std::align_val_t(alignof(walk_node<Dim, Simd>)));
   }
 }
 
 /* GETTERS, SETTERS, SIMPLE UTILITIES */
 
-template <int Dim> walk_node<Dim> *walk_tree<Dim>::root() const { return root_.get(); }
+template <int Dim, bool Simd> walk_node<Dim, Simd> *walk_tree<Dim, Simd>::root() const { return root_.get(); }
 
-template <int Dim> point<Dim> walk_tree<Dim>::endpoint() const { return root_->endpoint(); }
+template <int Dim, bool Simd> point<Dim, Simd> walk_tree<Dim, Simd>::endpoint() const { return root_->endpoint(); }
 
-template <int Dim> bool walk_tree<Dim>::is_leaf() const { return root_->is_leaf(); }
+template <int Dim, bool Simd> bool walk_tree<Dim, Simd>::is_leaf() const { return root_->is_leaf(); }
 
 /* PRIMITIVE OPERATIONS */
 
-template <int Dim> walk_node<Dim> &walk_tree<Dim>::find_node(int n) {
+template <int Dim, bool Simd> walk_node<Dim, Simd> &walk_tree<Dim, Simd>::find_node(int n) {
   if (!buf_) {
     throw std::runtime_error("find_node can only be used on trees initialized with balanced=true");
   }
-  walk_node<Dim> &result = buf_[n - 1];
+  walk_node<Dim, Simd> &result = buf_[n - 1];
   assert(result.id_ == n);
   return result;
 }
 
 /* HIGH-LEVEL FUNCTIONS */
 
-template <int Dim> bool walk_tree<Dim>::try_pivot(int n, const transform<Dim> &r) {
+template <int Dim, bool Simd> bool walk_tree<Dim, Simd>::try_pivot(int n, const transform<Dim, Simd> &r) {
   if (r.is_identity()) {
     return false;
   }
@@ -98,13 +99,13 @@ template <int Dim> bool walk_tree<Dim>::try_pivot(int n, const transform<Dim> &r
   return success;
 }
 
-template <int Dim> bool walk_tree<Dim>::try_pivot_fast(int n, const transform<Dim> &t) {
+template <int Dim, bool Simd> bool walk_tree<Dim, Simd>::try_pivot_fast(int n, const transform<Dim, Simd> &t) {
   if (t.is_identity()) {
     return false;
   }
 
-  walk_node<Dim> *w = &find_node(n); // TODO: a pointer seems to be needed, but why?
-  walk_node<Dim> w_copy(*w);
+  walk_node<Dim, Simd> *w = &find_node(n); // TODO: a pointer seems to be needed, but why?
+  walk_node<Dim, Simd> w_copy(*w);
   auto success = !w_copy.shuffle_intersect(t, w->is_left_child());
   if (success) {
     root_->shuffle_up(n);
@@ -115,17 +116,19 @@ template <int Dim> bool walk_tree<Dim>::try_pivot_fast(int n, const transform<Di
   return success;
 }
 
-template <int Dim> bool walk_tree<Dim>::rand_pivot(bool fast) {
+template <int Dim, bool Simd> bool walk_tree<Dim, Simd>::rand_pivot(bool fast) {
   auto site = dist_(rng_);
-  auto r = transform<Dim>::rand(rng_);
+  auto r = transform<Dim, Simd>::rand(rng_);
   return fast ? try_pivot_fast(site, r) : try_pivot(site, r);
 }
 
 /* OTHER FUNCTIONS */
 
-template <int Dim> std::vector<point<Dim>> walk_tree<Dim>::steps() const { return root_->steps(); }
+template <int Dim, bool Simd> std::vector<point<Dim, Simd>> walk_tree<Dim, Simd>::steps() const {
+  return root_->steps();
+}
 
-template <int Dim> bool walk_tree<Dim>::self_avoiding() const {
+template <int Dim, bool Simd> bool walk_tree<Dim, Simd>::self_avoiding() const {
   auto steps = this->steps();
   for (size_t i = 0; i < steps.size(); ++i) {
     for (size_t j = i + 1; j < steps.size(); ++j) {
@@ -137,9 +140,11 @@ template <int Dim> bool walk_tree<Dim>::self_avoiding() const {
   return true;
 }
 
-template <int Dim> void walk_tree<Dim>::export_csv(const std::string &path) const { return to_csv(path, steps()); }
+template <int Dim, bool Simd> void walk_tree<Dim, Simd>::export_csv(const std::string &path) const {
+  return to_csv(path, steps());
+}
 
-template <int Dim> void walk_tree<Dim>::todot(const std::string &path) const { root_->todot(path); }
+template <int Dim, bool Simd> void walk_tree<Dim, Simd>::todot(const std::string &path) const { root_->todot(path); }
 
 /* TEMPLATE INSTANTIATION */