Update clang-format configuration

24 files changed, 2314 insertions, 2388 deletions

diff --git a/.clang-format b/.clang-format
index a728264..8b29ee3 100644
--- a/.clang-format
+++ b/.clang-format
@@ -1,90 +1,25 @@
 ---
-Language:        Cpp
-# BasedOnStyle:  Mozilla
-AccessModifierOffset: -4
-AlignAfterOpenBracket: Align
 AlignConsecutiveAssignments: true
 AlignConsecutiveDeclarations: true
 AlignEscapedNewlinesLeft: true
-AlignOperands:   true
-AlignTrailingComments: true
-AllowAllParametersOfDeclarationOnNextLine: false
-AllowShortBlocksOnASingleLine: true
-AllowShortCaseLabelsOnASingleLine: true
-AllowShortFunctionsOnASingleLine: Inline
-AllowShortIfStatementsOnASingleLine: false
-AllowShortLoopsOnASingleLine: false
-AlwaysBreakAfterDefinitionReturnType: TopLevel
-AlwaysBreakAfterReturnType: TopLevelDefinitions
-AlwaysBreakBeforeMultilineStrings: false
-AlwaysBreakTemplateDeclarations: true
-BinPackArguments: false
-BinPackParameters: false
-BraceWrapping:   
-  AfterClass:      true
-  AfterControlStatement: false
-  AfterEnum:       false
-  AfterFunction:   true
-  AfterNamespace:  false
-  AfterObjCDeclaration: false
-  AfterStruct:     true
-  AfterUnion:      true
-  BeforeCatch:     false
-  BeforeElse:      false
-  IndentBraces:    false
-BreakBeforeBinaryOperators: None
+BasedOnStyle: Mozilla
+BraceWrapping:
+  AfterNamespace: false
+  AfterClass: true
+  AfterEnum: false
+  AfterExternBlock: false
+  AfterFunction: true
+  AfterStruct: false
+  SplitEmptyFunction: false
+  SplitEmptyRecord: false
 BreakBeforeBraces: Custom
-BreakBeforeTernaryOperators: true
-BreakConstructorInitializersBeforeComma: true
-ColumnLimit:     80
-CommentPragmas:  '^ IWYU pragma:'
-ConstructorInitializerAllOnOneLineOrOnePerLine: false
-ConstructorInitializerIndentWidth: 2
-ContinuationIndentWidth: 2
 Cpp11BracedListStyle: true
-DerivePointerAlignment: false
-DisableFormat:   false
-ExperimentalAutoDetectBinPacking: false
-ForEachMacros:   [ foreach, Q_FOREACH, BOOST_FOREACH ]
-IncludeCategories: 
-  - Regex:           '^"(llvm|llvm-c|clang|clang-c)/'
-    Priority:        2
-  - Regex:           '^(<|"(gtest|isl|json)/)'
-    Priority:        3
-  - Regex:           '.*'
-    Priority:        1
+FixNamespaceComments: true
 IndentCaseLabels: false
-IndentWidth:     4
-IndentWrappedFunctionNames: false
-KeepEmptyLinesAtTheStartOfBlocks: true
-MacroBlockBegin: ''
-MacroBlockEnd:   ''
-MaxEmptyLinesToKeep: 1
-NamespaceIndentation: None
-ObjCBlockIndentWidth: 2
-ObjCSpaceAfterProperty: true
-ObjCSpaceBeforeProtocolList: false
-PenaltyBreakBeforeFirstCallParameter: 19
-PenaltyBreakComment: 300
-PenaltyBreakFirstLessLess: 120
-PenaltyBreakString: 1000
-PenaltyExcessCharacter: 1000000
-PenaltyReturnTypeOnItsOwnLine: 200
-PointerAlignment: Left
-ReflowComments:  true
-SortIncludes:    true
-SpaceAfterCStyleCast: false
-SpaceBeforeAssignmentOperators: true
-SpaceBeforeParens: ControlStatements
-SpaceInEmptyParentheses: false
-SpacesBeforeTrailingComments: 1
-SpacesInAngles:  false
-SpacesInContainerLiterals: true
-SpacesInCStyleCastParentheses: false
-SpacesInParentheses: false
-SpacesInSquareBrackets: false
-Standard:        Cpp11
-TabWidth:        4
-UseTab:          ForIndentation
+IndentPPDirectives: AfterHash
+KeepEmptyLinesAtTheStartOfBlocks: false
+SpacesInContainerLiterals: false
+StatementMacros:
+  - _Pragma
+  - __pragma
 ...
-
diff --git a/benchmark/bench_bitvec.cpp b/benchmark/bench_bitvec.cpp
index 6580af9..5c885b1 100644
--- a/benchmark/bench_bitvec.cpp
+++ b/benchmark/bench_bitvec.cpp
@@ -29,324 +29,303 @@
 #include <string>
 #include <vector>
 
-template <class T, size_t N>
-struct BenchAnd
-{
-	Duration operator()(Context& ctx)
-	{
-		const T v = make_random_bitvec<T, N>(ctx);
-		T       r = make_random_bitvec<T, N>(ctx);
-
-		return run_bench([&](auto) { r &= v; });
-	}
+template<class T, size_t N>
+struct BenchAnd {
+  Duration operator()(Context& ctx)
+  {
+    const T v = make_random_bitvec<T, N>(ctx);
+    T       r = make_random_bitvec<T, N>(ctx);
+
+    return run_bench([&](auto) { r &= v; });
+  }
 };
 
-template <class T, size_t N>
-struct BenchOr
-{
-	Duration operator()(Context& ctx)
-	{
-		const T v = make_random_bitvec<T, N>(ctx);
-		T       r = make_random_bitvec<T, N>(ctx);
+template<class T, size_t N>
+struct BenchOr {
+  Duration operator()(Context& ctx)
+  {
+    const T v = make_random_bitvec<T, N>(ctx);
+    T       r = make_random_bitvec<T, N>(ctx);
 
-		return run_bench([&](auto) { r |= v; });
-	}
+    return run_bench([&](auto) { r |= v; });
+  }
 };
 
-template <class T, size_t N>
-struct BenchXor
-{
-	Duration operator()(Context& ctx)
-	{
-		const T v = make_random_bitvec<T, N>(ctx);
-		T       r = make_random_bitvec<T, N>(ctx);
+template<class T, size_t N>
+struct BenchXor {
+  Duration operator()(Context& ctx)
+  {
+    const T v = make_random_bitvec<T, N>(ctx);
+    T       r = make_random_bitvec<T, N>(ctx);
 
-		return run_bench([&](auto) { r ^= v; });
-	}
+    return run_bench([&](auto) { r ^= v; });
+  }
 };
 
-template <class T, size_t N>
-struct BenchSetOne
-{
-	Duration operator()(Context& ctx)
-	{
-		T v = make_random_bitvec<T, N>(ctx);
-		return run_bench([&](const auto i) { v.set(i % N); });
-	}
+template<class T, size_t N>
+struct BenchSetOne {
+  Duration operator()(Context& ctx)
+  {
+    T v = make_random_bitvec<T, N>(ctx);
+    return run_bench([&](const auto i) { v.set(i % N); });
+  }
 };
 
-template <class T, size_t N>
-struct BenchSetAll
-{
-	Duration operator()(Context& ctx)
-	{
-		T v = make_random_bitvec<T, N>(ctx);
-		return run_bench([&](auto) { v.set(); });
-	}
+template<class T, size_t N>
+struct BenchSetAll {
+  Duration operator()(Context& ctx)
+  {
+    T v = make_random_bitvec<T, N>(ctx);
+    return run_bench([&](auto) { v.set(); });
+  }
 };
 
-template <class T, size_t N>
-struct BenchResetOne
-{
-	Duration operator()(Context& ctx)
-	{
-		T v = make_random_bitvec<T, N>(ctx);
-		return run_bench([&](const auto i) { v.reset(i % N); });
-	}
+template<class T, size_t N>
+struct BenchResetOne {
+  Duration operator()(Context& ctx)
+  {
+    T v = make_random_bitvec<T, N>(ctx);
+    return run_bench([&](const auto i) { v.reset(i % N); });
+  }
 };
 
-template <class T, size_t N>
-struct BenchResetAll
-{
-	Duration operator()(Context& ctx)
-	{
-		T v = make_random_bitvec<T, N>(ctx);
-		return run_bench([&](auto) { v.reset(); });
-	}
+template<class T, size_t N>
+struct BenchResetAll {
+  Duration operator()(Context& ctx)
+  {
+    T v = make_random_bitvec<T, N>(ctx);
+    return run_bench([&](auto) { v.reset(); });
+  }
 };
 
-template <class T, size_t N>
-struct BenchFlipOne
-{
-	Duration operator()(Context& ctx)
-	{
-		T v = make_random_bitvec<T, N>(ctx);
-		return run_bench([&](const auto i) { v.flip(i % N); });
-	}
+template<class T, size_t N>
+struct BenchFlipOne {
+  Duration operator()(Context& ctx)
+  {
+    T v = make_random_bitvec<T, N>(ctx);
+    return run_bench([&](const auto i) { v.flip(i % N); });
+  }
 };
 
-template <class T, size_t N>
-struct BenchFlipAll
-{
-	Duration operator()(Context& ctx)
-	{
-		T v = make_random_bitvec<T, N>(ctx);
-		return run_bench([&](auto) { v.flip(); });
-	}
+template<class T, size_t N>
+struct BenchFlipAll {
+  Duration operator()(Context& ctx)
+  {
+    T v = make_random_bitvec<T, N>(ctx);
+    return run_bench([&](auto) { v.flip(); });
+  }
 };
 
-template <class T, size_t N>
-struct BenchNone
-{
-	Duration operator()(Context& ctx)
-	{
-		T    v = make_random_bitvec<T, N>(ctx);
-		bool r = false;
-		(void)r;
-
-		return run_bench([&](auto) { r = v.none(); });
-	}
+template<class T, size_t N>
+struct BenchNone {
+  Duration operator()(Context& ctx)
+  {
+    T    v = make_random_bitvec<T, N>(ctx);
+    bool r = false;
+    (void)r;
+
+    return run_bench([&](auto) { r = v.none(); });
+  }
 };
 
-template <class T, size_t N>
-struct BenchCount
-{
-	Duration operator()(Context& ctx)
-	{
-		T               v = make_random_bitvec<T, N>(ctx);
-		volatile size_t r = 0;
-		(void)r;
-
-		return run_bench([&](auto) { r = v.count(); });
-	}
+template<class T, size_t N>
+struct BenchCount {
+  Duration operator()(Context& ctx)
+  {
+    T               v = make_random_bitvec<T, N>(ctx);
+    volatile size_t r = 0;
+    (void)r;
+
+    return run_bench([&](auto) { r = v.count(); });
+  }
 };
 
-template <class T, size_t N>
-struct BenchLeftShift
-{
-	Duration operator()(Context& ctx)
-	{
-		const T v = make_random_bitvec<T, N>(ctx);
-		T       r = v;
-		return run_bench([&](const auto i) { r <<= (i % N); });
-	}
+template<class T, size_t N>
+struct BenchLeftShift {
+  Duration operator()(Context& ctx)
+  {
+    const T v = make_random_bitvec<T, N>(ctx);
+    T       r = v;
+    return run_bench([&](const auto i) { r <<= (i % N); });
+  }
 };
 
-template <class T, size_t N>
-struct BenchRightShift
-{
-	Duration operator()(Context& ctx)
-	{
-		const T v = make_random_bitvec<T, N>(ctx);
-		T       r = v;
-		return run_bench([&](const auto i) { r >>= (i % N); });
-	}
+template<class T, size_t N>
+struct BenchRightShift {
+  Duration operator()(Context& ctx)
+  {
+    const T v = make_random_bitvec<T, N>(ctx);
+    T       r = v;
+    return run_bench([&](const auto i) { r >>= (i % N); });
+  }
 };
 
-template <class T, size_t N>
-struct BenchLeftRotate
-{
-	Duration operator()(Context& ctx)
-	{
-		T v = make_random_bitvec<T, N>(ctx);
-		return run_bench([&](const auto i) { v.rotl(i % N); });
-	}
+template<class T, size_t N>
+struct BenchLeftRotate {
+  Duration operator()(Context& ctx)
+  {
+    T v = make_random_bitvec<T, N>(ctx);
+    return run_bench([&](const auto i) { v.rotl(i % N); });
+  }
 };
 
-template <class T, size_t N>
-struct BenchRightRotate
-{
-	Duration operator()(Context& ctx)
-	{
-		T             v = make_random_bitvec<T, N>(ctx);
-		volatile bool r = false;
-		(void)r;
-
-		return run_bench([&](const auto i) {
-			v.rotr(i % N);
-			r = v.test(0);
-		});
-	}
+template<class T, size_t N>
+struct BenchRightRotate {
+  Duration operator()(Context& ctx)
+  {
+    T             v = make_random_bitvec<T, N>(ctx);
+    volatile bool r = false;
+    (void)r;
+
+    return run_bench([&](const auto i) {
+      v.rotr(i % N);
+      r = v.test(0);
+    });
+  }
 };
 
-template <class T, size_t N>
-struct BenchFindFirst
-{
-	Duration operator()(Context& ctx)
-	{
-		T               v = make_random_bitvec<T, N>(ctx);
-		volatile size_t r = 0;
-		(void)r;
-
-		return run_bench([&](auto) { r = v.find_first(); });
-	}
+template<class T, size_t N>
+struct BenchFindFirst {
+  Duration operator()(Context& ctx)
+  {
+    T               v = make_random_bitvec<T, N>(ctx);
+    volatile size_t r = 0;
+    (void)r;
+
+    return run_bench([&](auto) { r = v.find_first(); });
+  }
 };
 
-template <class T, size_t N>
-struct BenchGrayCode
-{
-	Duration operator()(Context& ctx)
-	{
-		const T v = make_random_bitvec<T, N>(ctx);
-		T       r = v;
+template<class T, size_t N>
+struct BenchGrayCode {
+  Duration operator()(Context& ctx)
+  {
+    const T v = make_random_bitvec<T, N>(ctx);
+    T       r = v;
 
-		return run_bench([&](auto) { chilbert::detail::gray_code(r); });
-	}
+    return run_bench([&](auto) { chilbert::detail::gray_code(r); });
+  }
 };
 
-template <class T, size_t N>
-struct BenchGrayCodeInv
-{
-	Duration operator()(Context& ctx)
-	{
-		const T v = make_random_bitvec<T, N>(ctx);
-		T       r = v;
+template<class T, size_t N>
+struct BenchGrayCodeInv {
+  Duration operator()(Context& ctx)
+  {
+    const T v = make_random_bitvec<T, N>(ctx);
+    T       r = v;
 
-		return run_bench([&](auto) { chilbert::detail::gray_code_inv(r); });
-	}
+    return run_bench([&](auto) { chilbert::detail::gray_code_inv(r); });
+  }
 };
 
-template <class T, size_t N>
-struct BenchComparison
-{
-	Duration operator()(Context& ctx)
-	{
-		std::vector<T> vecs;
-		for (size_t i = 0; i < 32; ++i) {
-			vecs.emplace_back(make_random_bitvec<T, N>(ctx));
-		}
-
-		volatile bool r = false;
-
-		return run_bench([&](const auto i) {
-			r |= vecs[i % vecs.size()] < vecs[(i + 1) % vecs.size()];
-		});
-	}
+template<class T, size_t N>
+struct BenchComparison {
+  Duration operator()(Context& ctx)
+  {
+    std::vector<T> vecs;
+    for (size_t i = 0; i < 32; ++i) {
+      vecs.emplace_back(make_random_bitvec<T, N>(ctx));
+    }
+
+    volatile bool r = false;
+
+    return run_bench([&](const auto i) {
+      r |= vecs[i % vecs.size()] < vecs[(i + 1) % vecs.size()];
+    });
+  }
 };
 
-template <class T, size_t N>
-struct BenchIteration
-{
-	Duration operator()(Context& ctx)
-	{
-		T             v = make_random_bitvec<T, N>(ctx);
-		auto          i = v.begin();
-		volatile bool r = false;
-		return run_bench([&](const auto) {
-			r = (++i == v.end());
-			if (r) {
-				i = v.begin();
-			}
-		});
-	}
+template<class T, size_t N>
+struct BenchIteration {
+  Duration operator()(Context& ctx)
+  {
+    T             v = make_random_bitvec<T, N>(ctx);
+    auto          i = v.begin();
+    volatile bool r = false;
+    return run_bench([&](const auto) {
+      r = (++i == v.end());
+      if (r) {
+        i = v.begin();
+      }
+    });
+  }
 };
 
 /// Run benchmark for size N
-template <template <class U, size_t M> class Bench, size_t N>
+template<template<class U, size_t M> class Bench, size_t N>
 void
 bench_row(Context& ctx, std::ostream& os)
 {
-	std::array<Duration, 4> results = {
-	  ((N <= chilbert::SmallBitVec::bits_per_rack)
-	     ? Bench<chilbert::SmallBitVec, N>{}(ctx)
-	     : Duration{}),
-	  Bench<chilbert::StaticBitVec<N>, N>{}(ctx),
-	  Bench<chilbert::BoundedBitVec<N>, N>{}(ctx),
-	  Bench<chilbert::DynamicBitVec, N>{}(ctx),
-	};
-
-	write_row(os, N, results);
+  std::array<Duration, 4> results = {
+    ((N <= chilbert::SmallBitVec::bits_per_rack)
+       ? Bench<chilbert::SmallBitVec, N>{}(ctx)
+       : Duration{}),
+    Bench<chilbert::StaticBitVec<N>, N>{}(ctx),
+    Bench<chilbert::BoundedBitVec<N>, N>{}(ctx),
+    Bench<chilbert::DynamicBitVec, N>{}(ctx),
+  };
+
+  write_row(os, N, results);
 }
 
 /// Terminate recursion
-template <template <class U, size_t M> class Bench>
+template<template<class U, size_t M> class Bench>
 void
 bench_rec(Context&, std::ostream&)
-{
-}
+{}
 
 /// Run benchmark for sizes N, Ns... (recursive helper)
-template <template <class U, size_t M> class Bench, size_t N, size_t... Ns>
+template<template<class U, size_t M> class Bench, size_t N, size_t... Ns>
 void
 bench_rec(Context& ctx, std::ostream& os)
 {
-	bench_row<Bench, N>(ctx, os);
-	bench_rec<Bench, Ns...>(ctx, os);
+  bench_row<Bench, N>(ctx, os);
+  bench_rec<Bench, Ns...>(ctx, os);
 }
 
 /// Run benchmark
-template <template <class U, size_t M> class Bench>
+template<template<class U, size_t M> class Bench>
 void
 bench(Context& ctx, const std::string& name)
 {
-	std::ofstream out("bitvec_" + name + ".txt");
-	out << "n\tsmall\tstatic\tbounded\tdynamic\n";
-	bench_rec<Bench, 16, 32, 64, 128, 256, 512>(ctx, out);
+  std::ofstream out("bitvec_" + name + ".txt");
+  out << "n\tsmall\tstatic\tbounded\tdynamic\n";
+  bench_rec<Bench, 16, 32, 64, 128, 256, 512>(ctx, out);
 }
 
 int
 main()
 {
-	Context ctx;
+  Context ctx;
 
-	bench<BenchAnd>(ctx, "and");
-	bench<BenchOr>(ctx, "or");
-	bench<BenchXor>(ctx, "xor");
-	bench<BenchNone>(ctx, "none");
+  bench<BenchAnd>(ctx, "and");
+  bench<BenchOr>(ctx, "or");
+  bench<BenchXor>(ctx, "xor");
+  bench<BenchNone>(ctx, "none");
 
-	bench<BenchSetOne>(ctx, "set_one");
-	bench<BenchSetAll>(ctx, "set_all");
+  bench<BenchSetOne>(ctx, "set_one");
+  bench<BenchSetAll>(ctx, "set_all");
 
-	bench<BenchResetOne>(ctx, "reset_one");
-	bench<BenchResetAll>(ctx, "reset_all");
+  bench<BenchResetOne>(ctx, "reset_one");
+  bench<BenchResetAll>(ctx, "reset_all");
 
-	bench<BenchFlipOne>(ctx, "flip_one");
-	bench<BenchFlipAll>(ctx, "flip_all");
+  bench<BenchFlipOne>(ctx, "flip_one");
+  bench<BenchFlipAll>(ctx, "flip_all");
 
-	bench<BenchCount>(ctx, "count");
-	bench<BenchFindFirst>(ctx, "find_first");
+  bench<BenchCount>(ctx, "count");
+  bench<BenchFindFirst>(ctx, "find_first");
 
-	bench<BenchLeftShift>(ctx, "left_shift");
-	bench<BenchRightShift>(ctx, "right_shift");
+  bench<BenchLeftShift>(ctx, "left_shift");
+  bench<BenchRightShift>(ctx, "right_shift");
 
-	bench<BenchLeftRotate>(ctx, "left_rotate");
-	bench<BenchRightRotate>(ctx, "right_rotate");
+  bench<BenchLeftRotate>(ctx, "left_rotate");
+  bench<BenchRightRotate>(ctx, "right_rotate");
 
-	bench<BenchGrayCode>(ctx, "gray_code");
-	bench<BenchGrayCodeInv>(ctx, "gray_code_inv");
+  bench<BenchGrayCode>(ctx, "gray_code");
+  bench<BenchGrayCodeInv>(ctx, "gray_code_inv");
 
-	bench<BenchComparison>(ctx, "comparison");
-	bench<BenchIteration>(ctx, "iteration");
+  bench<BenchComparison>(ctx, "comparison");
+  bench<BenchIteration>(ctx, "iteration");
 
-	return 0;
+  return 0;
 }
diff --git a/benchmark/bench_hilbert.cpp b/benchmark/bench_hilbert.cpp
index 52c12dc..5793d29 100644
--- a/benchmark/bench_hilbert.cpp
+++ b/benchmark/bench_hilbert.cpp
@@ -29,85 +29,80 @@
 #include <fstream>
 #include <string>
 
-template <class H, size_t M, size_t D>
-struct BenchCoordsToIndex
-{
-	Duration operator()(Context& ctx)
-	{
-		const auto p  = make_random_point<M, D>(ctx);
-		H          ha = make_zero_bitvec<H, D * M>();
-
-		return run_bench([&](auto) { chilbert::coords_to_index(p, M, D, ha); });
-	}
+template<class H, size_t M, size_t D>
+struct BenchCoordsToIndex {
+  Duration operator()(Context& ctx)
+  {
+    const auto p  = make_random_point<M, D>(ctx);
+    H          ha = make_zero_bitvec<H, D * M>();
+
+    return run_bench([&](auto) { chilbert::coords_to_index(p, M, D, ha); });
+  }
 };
 
-template <class H, size_t M, size_t D>
-struct BenchIndexToCoords
-{
-	Duration operator()(Context& ctx)
-	{
-		auto    p  = make_random_point<M, D>(ctx);
-		const H ha = make_random_bitvec<H, D * M>(ctx);
+template<class H, size_t M, size_t D>
+struct BenchIndexToCoords {
+  Duration operator()(Context& ctx)
+  {
+    auto    p  = make_random_point<M, D>(ctx);
+    const H ha = make_random_bitvec<H, D * M>(ctx);
 
-		return run_bench([&](auto) { chilbert::index_to_coords(p, M, D, ha); });
-	}
+    return run_bench([&](auto) { chilbert::index_to_coords(p, M, D, ha); });
+  }
 };
 
 /// Run benchmark for size N
-template <template <class H, size_t M, size_t D> class Bench,
-          size_t M,
-          size_t D>
+template<template<class H, size_t M, size_t D> class Bench, size_t M, size_t D>
 void
 bench_row(Context& ctx, std::ostream& os)
 {
-	std::array<Duration, 4> results = {
-	  ((M * D <= chilbert::SmallBitVec::bits_per_rack)
-	     ? Bench<chilbert::SmallBitVec, M, D>{}(ctx)
-	     : Duration{}),
-	  Bench<chilbert::StaticBitVec<M * D>, M, D>{}(ctx),
-	  Bench<chilbert::BoundedBitVec<M * D>, M, D>{}(ctx),
-	  Bench<chilbert::DynamicBitVec, M, D>{}(ctx),
-	};
-
-	write_row(os, D, results);
+  std::array<Duration, 4> results = {
+    ((M * D <= chilbert::SmallBitVec::bits_per_rack)
+       ? Bench<chilbert::SmallBitVec, M, D>{}(ctx)
+       : Duration{}),
+    Bench<chilbert::StaticBitVec<M * D>, M, D>{}(ctx),
+    Bench<chilbert::BoundedBitVec<M * D>, M, D>{}(ctx),
+    Bench<chilbert::DynamicBitVec, M, D>{}(ctx),
+  };
+
+  write_row(os, D, results);
 }
 
 /// Terminate recursion
-template <template <class H, size_t M, size_t D> class Bench, size_t M>
+template<template<class H, size_t M, size_t D> class Bench, size_t M>
 void
 bench_rec(Context&, std::ostream&)
-{
-}
+{}
 
 /// Run benchmark for sizes N, Ns... (recursive helper)
-template <template <class H, size_t M, size_t D> class Bench,
-          size_t M,
-          size_t D,
-          size_t... Ds>
+template<template<class H, size_t M, size_t D> class Bench,
+         size_t M,
+         size_t D,
+         size_t... Ds>
 void
 bench_rec(Context& ctx, std::ostream& os)
 {
-	bench_row<Bench, M, D>(ctx, os);
-	bench_rec<Bench, M, Ds...>(ctx, os);
+  bench_row<Bench, M, D>(ctx, os);
+  bench_rec<Bench, M, Ds...>(ctx, os);
 }
 
 /// Run benchmark
-template <template <class H, size_t M, size_t D> class Bench>
+template<template<class H, size_t M, size_t D> class Bench>
 void
 bench(Context& ctx, const std::string& name)
 {
-	std::ofstream out("hilbert_" + name + ".txt");
-	out << "d\tsmall\tstatic\tbounded\tdynamic\n";
-	bench_rec<Bench, 8, 2, 4, 8, 16, 32, 64>(ctx, out);
+  std::ofstream out("hilbert_" + name + ".txt");
+  out << "d\tsmall\tstatic\tbounded\tdynamic\n";
+  bench_rec<Bench, 8, 2, 4, 8, 16, 32, 64>(ctx, out);
 }
 
 int
 main()
 {
-	Context ctx;
+  Context ctx;
 
-	bench<BenchCoordsToIndex>(ctx, "coords_to_index");
-	bench<BenchIndexToCoords>(ctx, "index_to_coords");
+  bench<BenchCoordsToIndex>(ctx, "coords_to_index");
+  bench<BenchIndexToCoords>(ctx, "index_to_coords");
 
-	return 0;
+  return 0;
 }
diff --git a/benchmark/bench_utils.hpp b/benchmark/bench_utils.hpp
index a3f7f81..620cb27 100644
--- a/benchmark/bench_utils.hpp
+++ b/benchmark/bench_utils.hpp
@@ -28,37 +28,37 @@
 using Duration = std::chrono::duration<double, std::micro>;
 
 /// Write a TSV row to `os` with `n` as the first column followed by `results`
-template <class T, size_t M>
+template<class T, size_t M>
 void
 write_row(std::ostream& os, const size_t n, const std::array<T, M>& results)
 {
-	os << n;
-	for (const auto t : results) {
-		if (t == Duration::zero()) {
-			os << "\tNaN";
-		} else {
-			os << '\t' << t.count();
-		}
-	}
-	os << std::endl;
+  os << n;
+  for (const auto t : results) {
+    if (t == Duration::zero()) {
+      os << "\tNaN";
+    } else {
+      os << '\t' << t.count();
+    }
+  }
+  os << std::endl;
 }
 
 /// Repeatedly run an operation and return the average time
-template <class Operation>
+template<class Operation>
 Duration
 run_bench(const Operation& op)
 {
-	static constexpr auto bench_duration = std::chrono::milliseconds{10};
+  static constexpr auto bench_duration = std::chrono::milliseconds{10};
 
-	const auto t_start = std::chrono::steady_clock{}.now();
-	auto       t_now   = t_start;
-	size_t     count   = 0;
-	for (; t_now < t_start + bench_duration; ++count) {
-		op(count);
-		t_now = std::chrono::steady_clock{}.now();
-	}
+  const auto t_start = std::chrono::steady_clock{}.now();
+  auto       t_now   = t_start;
+  size_t     count   = 0;
+  for (; t_now < t_start + bench_duration; ++count) {
+    op(count);
+    t_now = std::chrono::steady_clock{}.now();
+  }
 
-	return (t_now - t_start) / count;
+  return (t_now - t_start) / count;
 }
 
 #endif
diff --git a/include/chilbert/BoundedBitVec.hpp b/include/chilbert/BoundedBitVec.hpp
index 4f9ef60..6fe5d2a 100644
--- a/include/chilbert/BoundedBitVec.hpp
+++ b/include/chilbert/BoundedBitVec.hpp
@@ -41,81 +41,103 @@ namespace chilbert {
  *
  * @tparam MaxN Maximum number of bits.
  */
-template <size_t MaxN>
+template<size_t MaxN>
 class BoundedBitVec : public detail::MultiBitVec<BoundedBitVec<MaxN>>
 {
 public:
-	using Rack = typename detail::MultiBitVec<BoundedBitVec<MaxN>>::Rack;
+  using Rack = typename detail::MultiBitVec<BoundedBitVec<MaxN>>::Rack;
 
-	using detail::MultiBitVec<BoundedBitVec<MaxN>>::bits_per_rack;
+  using detail::MultiBitVec<BoundedBitVec<MaxN>>::bits_per_rack;
 
-	BoundedBitVec() = default;
+  BoundedBitVec() = default;
 
-	explicit BoundedBitVec(const size_t bits)
-	  : m_size{bits}
-	{
-		assert(bits <= MaxN);
-	}
+  explicit BoundedBitVec(const size_t bits)
+    : m_size{bits}
+  {
+    assert(bits <= MaxN);
+  }
 
-	BoundedBitVec(const size_t bits, const Rack value)
-	  : BoundedBitVec{bits}
-	{
-		m_racks[0] = value;
-	}
+  BoundedBitVec(const size_t bits, const Rack value)
+    : BoundedBitVec{bits}
+  {
+    m_racks[0] = value;
+  }
 
-	/// Return the size in bits
-	size_t size() const { return m_size; }
+  /// Return the size in bits
+  size_t size() const { return m_size; }
 
-	/// Return a reference to the `index`th rack
+  /// Return a reference to the `index`th rack
 #ifndef NDEBUG
-	const auto& rack(const size_t index) const { return m_racks.at(index); }
-	auto&       rack(const size_t index) { return m_racks.at(index); }
+  const auto& rack(const size_t index) const
+  {
+    return m_racks.at(index);
+  }
+  auto& rack(const size_t index)
+  {
+    return m_racks.at(index);
+  }
 #else
-	const auto& rack(const size_t index) const { return m_racks[index]; }
-	auto&       rack(const size_t index) { return m_racks[index]; }
+  const auto& rack(const size_t index) const
+  {
+    return m_racks[index];
+  }
+  auto& rack(const size_t index)
+  {
+    return m_racks[index];
+  }
 #endif
 
-	/// Return a raw pointer to the racks
-	Rack*       data() { return m_racks.data(); }
-	const Rack* data() const { return m_racks.data(); }
-
-	/// Return the total size of all racks in bytes
-	size_t data_size() const { return num_racks() * sizeof(Rack); }
-
-	/// Return the number of racks
-	size_t num_racks() const { return calculate_num_racks(m_size); }
+  /// Return a raw pointer to the racks
+  Rack* data()
+  {
+    return m_racks.data();
+  }
+  const Rack* data() const
+  {
+    return m_racks.data();
+  }
+
+  /// Return the total size of all racks in bytes
+  size_t data_size() const
+  {
+    return num_racks() * sizeof(Rack);
+  }
+
+  /// Return the number of racks
+  size_t num_racks() const
+  {
+    return calculate_num_racks(m_size);
+  }
 
 private:
-	static constexpr size_t calculate_num_racks(const size_t bits)
-	{
-		return (std::max(bits, size_t(1)) + bits_per_rack - 1) / bits_per_rack;
-	}
+  static constexpr size_t calculate_num_racks(const size_t bits)
+  {
+    return (std::max(bits, size_t(1)) + bits_per_rack - 1) / bits_per_rack;
+  }
 
-	std::array<Rack, calculate_num_racks(MaxN)> m_racks{};
-	size_t                                      m_size{};
+  std::array<Rack, calculate_num_racks(MaxN)> m_racks{};
+  size_t                                      m_size{};
 };
 
 namespace detail {
 
-template <size_t MaxN>
-struct is_bitvec<BoundedBitVec<MaxN>>
-{
-	constexpr static bool value = true;
+template<size_t MaxN>
+struct is_bitvec<BoundedBitVec<MaxN>> {
+  constexpr static bool value = true;
 };
 
-template <size_t MaxN>
+template<size_t MaxN>
 void
 gray_code(BoundedBitVec<MaxN>& value)
 {
-	gray_code(static_cast<detail::MultiBitVec<BoundedBitVec<MaxN>>&>(value));
+  gray_code(static_cast<detail::MultiBitVec<BoundedBitVec<MaxN>>&>(value));
 }
 
-template <size_t MaxN>
+template<size_t MaxN>
 void
 gray_code_inv(BoundedBitVec<MaxN>& value)
 {
-	gray_code_inv(
-	  static_cast<detail::MultiBitVec<BoundedBitVec<MaxN>>&>(value));
+  gray_code_inv(static_cast<detail::MultiBitVec<BoundedBitVec<MaxN>>&>(value));
 }
 
 } // namespace detail
diff --git a/include/chilbert/DynamicBitVec.hpp b/include/chilbert/DynamicBitVec.hpp
index fe7a535..895860c 100644
--- a/include/chilbert/DynamicBitVec.hpp
+++ b/include/chilbert/DynamicBitVec.hpp
@@ -41,117 +41,113 @@ namespace chilbert {
 class DynamicBitVec : public detail::MultiBitVec<DynamicBitVec>
 {
 public:
-	struct RacksDeleter
-	{
-		void operator()(Rack* const racks) { free(racks); }
-	};
-
-	struct NullDeleter
-	{
-		void operator()(const Rack* const) {}
-	};
-
-	using RacksPtr      = std::unique_ptr<Rack[], RacksDeleter>;
-	using ConstRacksPtr = std::unique_ptr<const Rack[], NullDeleter>;
-
-	explicit DynamicBitVec(const size_t bits)
-	  : m_racks{make_racks(calculate_num_racks(bits))}
-	  , m_size{bits}
-	{
-	}
-
-	DynamicBitVec(const size_t bits, const Rack value)
-	  : DynamicBitVec{bits}
-	{
-		m_racks[0] = value;
-	}
-
-	DynamicBitVec(const DynamicBitVec& vec)
-	  : m_racks{make_racks(vec.num_racks())}
-	  , m_size{vec.m_size}
-	{
-		if (vec.data()) {
-			memcpy(data(), vec.data(), data_size());
-		}
-	}
-
-	DynamicBitVec(DynamicBitVec&&) = default;
-
-	DynamicBitVec& operator=(const DynamicBitVec& vec)
-	{
-		if (&vec != this) {
-			if (num_racks() < vec.num_racks()) {
-				m_racks = make_racks(vec.num_racks());
-				m_size  = vec.m_size;
-				memcpy(data(), vec.data(), data_size());
-			} else if (vec.num_racks() > 0) {
-				m_size = vec.m_size;
-				memcpy(data(), vec.data(), data_size());
-			} else {
-				m_size = 0;
-				m_racks.reset();
-			}
-		}
-
-		return *this;
-	}
-
-	DynamicBitVec& operator=(DynamicBitVec&&) = default;
-
-	~DynamicBitVec() = default;
-
-	/// Return the size in bits
-	size_t size() const { return m_size; }
-
-	/// Return a reference to the `index`th rack
-	const Rack& rack(const size_t index) const { return m_racks[index]; }
-	Rack&       rack(const size_t index) { return m_racks[index]; }
-
-	/// Return a raw pointer to the racks
-	Rack*       data() { return m_racks.get(); }
-	const Rack* data() const { return m_racks.get(); }
-
-	/// Return the total size of all racks in bytes
-	size_t data_size() const { return num_racks() * sizeof(Rack); }
-
-	/// Return the number of racks
-	size_t num_racks() const { return calculate_num_racks(m_size); }
+  struct RacksDeleter {
+    void operator()(Rack* const racks) { free(racks); }
+  };
+
+  struct NullDeleter {
+    void operator()(const Rack* const) {}
+  };
+
+  using RacksPtr      = std::unique_ptr<Rack[], RacksDeleter>;
+  using ConstRacksPtr = std::unique_ptr<const Rack[], NullDeleter>;
+
+  explicit DynamicBitVec(const size_t bits)
+    : m_racks{make_racks(calculate_num_racks(bits))}
+    , m_size{bits}
+  {}
+
+  DynamicBitVec(const size_t bits, const Rack value)
+    : DynamicBitVec{bits}
+  {
+    m_racks[0] = value;
+  }
+
+  DynamicBitVec(const DynamicBitVec& vec)
+    : m_racks{make_racks(vec.num_racks())}
+    , m_size{vec.m_size}
+  {
+    if (vec.data()) {
+      memcpy(data(), vec.data(), data_size());
+    }
+  }
+
+  DynamicBitVec(DynamicBitVec&&) = default;
+
+  DynamicBitVec& operator=(const DynamicBitVec& vec)
+  {
+    if (&vec != this) {
+      if (num_racks() < vec.num_racks()) {
+        m_racks = make_racks(vec.num_racks());
+        m_size  = vec.m_size;
+        memcpy(data(), vec.data(), data_size());
+      } else if (vec.num_racks() > 0) {
+        m_size = vec.m_size;
+        memcpy(data(), vec.data(), data_size());
+      } else {
+        m_size = 0;
+        m_racks.reset();
+      }
+    }
+
+    return *this;
+  }
+
+  DynamicBitVec& operator=(DynamicBitVec&&) = default;
+
+  ~DynamicBitVec() = default;
+
+  /// Return the size in bits
+  size_t size() const { return m_size; }
+
+  /// Return a reference to the `index`th rack
+  const Rack& rack(const size_t index) const { return m_racks[index]; }
+  Rack&       rack(const size_t index) { return m_racks[index]; }
+
+  /// Return a raw pointer to the racks
+  Rack*       data() { return m_racks.get(); }
+  const Rack* data() const { return m_racks.get(); }
+
+  /// Return the total size of all racks in bytes
+  size_t data_size() const { return num_racks() * sizeof(Rack); }
+
+  /// Return the number of racks
+  size_t num_racks() const { return calculate_num_racks(m_size); }
 
 private:
-	static size_t calculate_num_racks(const size_t bits)
-	{
-		return (std::max(bits, size_t(1)) + bits_per_rack - 1) / bits_per_rack;
-	}
-
-	static RacksPtr make_racks(const size_t n)
-	{
-		return RacksPtr{static_cast<Rack*>(calloc(n, sizeof(Rack)))};
-	}
-
-	RacksPtr m_racks;
-	size_t   m_size;
+  static size_t calculate_num_racks(const size_t bits)
+  {
+    return (std::max(bits, size_t(1)) + bits_per_rack - 1) / bits_per_rack;
+  }
+
+  static RacksPtr make_racks(const size_t n)
+  {
+    return RacksPtr{static_cast<Rack*>(calloc(n, sizeof(Rack)))};
+  }
+
+  RacksPtr m_racks;
+  size_t   m_size;
 };
 
 namespace detail {
 
-template <>
-struct is_bitvec<DynamicBitVec>
-{
-	constexpr static bool value = true;
+template<>
+struct is_bitvec<DynamicBitVec> {
+  constexpr static bool value = true;
 };
 
-template <>
+template<>
 inline void
 gray_code(DynamicBitVec& value)
 {
-	gray_code(static_cast<MultiBitVec<DynamicBitVec>&>(value));
+  gray_code(static_cast<MultiBitVec<DynamicBitVec>&>(value));
 }
 
-template <>
+template<>
 inline void
 gray_code_inv(DynamicBitVec& value)
 {
-	gray_code_inv(static_cast<MultiBitVec<DynamicBitVec>&>(value));
+  gray_code_inv(static_cast<MultiBitVec<DynamicBitVec>&>(value));
 }
 
 } // namespace detail
diff --git a/include/chilbert/SmallBitVec.hpp b/include/chilbert/SmallBitVec.hpp
index 6f04965..485ca9e 100644
--- a/include/chilbert/SmallBitVec.hpp
+++ b/include/chilbert/SmallBitVec.hpp
@@ -33,336 +33,316 @@ namespace chilbert {
 class SmallBitVec
 {
 public:
-	using Rack = uintptr_t;
-
-	static constexpr size_t bits_per_rack = sizeof(Rack) * CHAR_BIT;
-
-	/** Mask for a bit that can be incremented like an index.
-	 *
-	 * This enables fast iteration while setting or resetting bits since it is
-	 * internally stored as a mask which avoids repeated index math.
-	 */
-	class Mask
-	{
-	public:
-		void operator++() { m_mask <<= 1U; }
-		void operator--() { m_mask >>= 1U; }
-
-		bool operator==(const Mask& mask) const
-		{
-			return m_mask == mask.m_mask;
-		}
-
-		bool operator!=(const Mask& mask) const
-		{
-			return m_mask == mask.m_mask;
-		}
-
-	private:
-		friend class SmallBitVec;
-
-		explicit Mask(const size_t index)
-		  : m_mask{index < bits_per_rack ? Rack{1} << index : 0}
-		{
-		}
-
-		Rack m_mask;
-	};
-
-	explicit SmallBitVec(const size_t bits)
-	  : m_size{bits}
-	{
-		assert(bits <= bits_per_rack);
-	}
-
-	SmallBitVec(const size_t bits, const Rack value)
-	  : m_rack{value}
-	  , m_size{bits}
-	{
-		assert(bits <= bits_per_rack);
-	}
-
-	/// Return the size in bits
-	size_t size() const { return m_size; }
-
-	/// Set all bits to one
-	SmallBitVec& set()
-	{
-		m_rack = ~Rack{0} >> (bits_per_rack - m_size);
-		return *this;
-	}
-
-	/// Set all bits to zero
-	SmallBitVec& reset()
-	{
-		m_rack = 0;
-		return *this;
-	}
-
-	/// Return the value of the bit covered by `mask`
-	bool test(const Mask mask) const { return m_rack & mask.m_mask; }
-
-	/// Return the value of the `index`th bit
-	bool test(const size_t index) const { return test(mask(index)); }
-
-	/// Set the bit covered by `mask` to 1
-	SmallBitVec& set(const Mask mask)
-	{
-		m_rack |= mask.m_mask;
-		return *this;
-	}
-
-	/// Set the `index`th bit to 1
-	SmallBitVec& set(const size_t index) { return set(mask(index)); }
-
-	/// Reset the bit covered by `mask` to 0
-	SmallBitVec& reset(const Mask mask)
-	{
-		m_rack &= ~mask.m_mask;
-		return *this;
-	}
-
-	/// Reset the `index`th bit to 0
-	SmallBitVec& reset(const size_t index) { return reset(mask(index)); }
-
-	/// Set the bit covered by `mask` to `value`
-	SmallBitVec& set(const Mask mask, const bool value)
-	{
-		m_rack ^= (-Rack{value} ^ m_rack) & mask.m_mask;
-		return *this;
-	}
-
-	/// Set the `index`th bit to `value`
-	SmallBitVec& set(const size_t index, const bool value)
-	{
-		return set(mask(index), value);
-	}
-
-	/// Flip the value of the bit covered by `mask`
-	SmallBitVec& flip(const Mask mask)
-	{
-		m_rack ^= mask.m_mask;
-		return *this;
-	}
-
-	/// Flip the value of the `index`th bit
-	SmallBitVec& flip(const size_t index) { return flip(mask(index)); }
-
-	bool operator==(const SmallBitVec& vec) const
-	{
-		return m_rack == vec.m_rack;
-	}
-
-	bool operator!=(const SmallBitVec& vec) const
-	{
-		return m_rack != vec.m_rack;
-	}
-
-	bool operator<(const SmallBitVec& vec) const { return m_rack < vec.m_rack; }
-
-	SmallBitVec& operator=(const Rack i)
-	{
-		m_rack = i;
-		return *this;
-	}
-
-	SmallBitVec& operator&=(const SmallBitVec& vec)
-	{
-		m_rack &= vec.m_rack;
-		return *this;
-	}
-
-	SmallBitVec& operator|=(const SmallBitVec& vec)
-	{
-		m_rack |= vec.m_rack;
-		return *this;
-	}
-
-	SmallBitVec& operator^=(const SmallBitVec& vec)
-	{
-		m_rack ^= vec.m_rack;
-		return *this;
-	}
-
-	SmallBitVec& operator^=(const Rack i)
-	{
-		m_rack ^= i;
-		return *this;
-	}
-
-	SmallBitVec& operator<<=(const size_t bits)
-	{
-		assert(bits < size());
-		m_rack <<= bits;
-		return *this;
-	}
-
-	SmallBitVec& operator>>=(const size_t bits)
-	{
-		assert(bits < size());
-		m_rack >>= bits;
-		return *this;
-	}
-
-	/// Right-rotate by `bits` positions
-	SmallBitVec& rotr(const size_t bits)
-	{
-		if (bits > 0 && bits < size()) {
-			assert(bits <= bits_per_rack);
-			m_rack = (m_rack >> bits) | (m_rack << (size() - bits));
-			m_rack &= (~Rack{0} >> (bits_per_rack - size()));
-		}
-		return *this;
-	}
-
-	/// Left-rotate by `bits` positions
-	SmallBitVec& rotl(const size_t bits)
-	{
-		if (bits > 0 && bits < size()) {
-			assert(bits <= bits_per_rack);
-			m_rack = (m_rack << bits) | (m_rack >> (size() - bits));
-			m_rack &= (~Rack{0} >> (bits_per_rack - size()));
-		}
-		return *this;
-	}
-
-	/// Return true iff all bits are zero
-	bool none() const { return m_rack == 0; }
-
-	/// Return 1 + the index of the first set bit, or 0 if there are none
-	size_t find_first() const
-	{
-		return static_cast<size_t>(detail::find_first(m_rack));
-	}
-
-	/// Return the number of set bits
-	size_t count() const
-	{
-		return static_cast<size_t>(detail::pop_count(m_rack));
-	}
-
-	/// Flip all bits (one's complement)
-	SmallBitVec& flip()
-	{
-		m_rack = ~m_rack;
-		return *this;
-	}
-
-	/// Return the first rack
-	Rack& rack() { return m_rack; }
-	Rack  rack() const { return m_rack; }
-
-	/// Return a raw pointer to the racks
-	Rack*       data() { return &m_rack; }
-	const Rack* data() const { return &m_rack; }
-
-	/// Return the number of racks
-	size_t num_racks() const { return 1; }
-
-	template <class BitVec>
-	class iterator_base : public Mask
-	{
-	public:
-		iterator_base& operator++()
-		{
-			Mask::operator++();
-			return *this;
-		}
-
-		iterator_base& operator--()
-		{
-			Mask::operator--();
-			return *this;
-		}
-
-		bool operator==(const iterator_base& rhs) const
-		{
-			return m_vec == rhs.m_vec && Mask::operator==(rhs);
-		}
-
-		bool operator!=(const iterator_base& rhs) const
-		{
-			return !operator==(rhs);
-		}
-
-		bool operator*() const { return m_vec->test(*this); }
-
-	protected:
-		iterator_base(BitVec& vec, const size_t index)
-		  : Mask{index}
-		  , m_vec{&vec}
-		{
-		}
-
-		BitVec* m_vec;
-	};
-
-	class iterator : public iterator_base<SmallBitVec>
-	{
-	public:
-		void set() { m_vec->set(*this); }
-		void reset() { m_vec->reset(*this); }
-
-	private:
-		friend class SmallBitVec;
-
-		iterator(SmallBitVec& vec, const size_t index)
-		  : iterator_base{vec, index}
-		{
-		}
-	};
-
-	class const_iterator : public iterator_base<const SmallBitVec>
-	{
-	private:
-		friend class SmallBitVec;
-
-		const_iterator(const SmallBitVec& vec, const size_t index)
-		  : iterator_base{vec, index}
-		{
-		}
-	};
-
-	Mask mask(const size_t i = 0) const
-	{
-		assert(i <= size());
-		return Mask{i};
-	}
-
-	iterator begin(const size_t i = 0) { return {*this, i}; }
-	iterator end() { return {*this, size()}; }
-
-	const_iterator begin(const size_t i = 0) const { return {*this, i}; }
-	const_iterator end() const { return {*this, size()}; }
+  using Rack = uintptr_t;
+
+  static constexpr size_t bits_per_rack = sizeof(Rack) * CHAR_BIT;
+
+  /** Mask for a bit that can be incremented like an index.
+   *
+   * This enables fast iteration while setting or resetting bits since it is
+   * internally stored as a mask which avoids repeated index math.
+   */
+  class Mask
+  {
+  public:
+    void operator++() { m_mask <<= 1U; }
+    void operator--() { m_mask >>= 1U; }
+
+    bool operator==(const Mask& mask) const { return m_mask == mask.m_mask; }
+
+    bool operator!=(const Mask& mask) const { return m_mask == mask.m_mask; }
+
+  private:
+    friend class SmallBitVec;
+
+    explicit Mask(const size_t index)
+      : m_mask{index < bits_per_rack ? Rack{1} << index : 0}
+    {}
+
+    Rack m_mask;
+  };
+
+  explicit SmallBitVec(const size_t bits)
+    : m_size{bits}
+  {
+    assert(bits <= bits_per_rack);
+  }
+
+  SmallBitVec(const size_t bits, const Rack value)
+    : m_rack{value}
+    , m_size{bits}
+  {
+    assert(bits <= bits_per_rack);
+  }
+
+  /// Return the size in bits
+  size_t size() const { return m_size; }
+
+  /// Set all bits to one
+  SmallBitVec& set()
+  {
+    m_rack = ~Rack{0} >> (bits_per_rack - m_size);
+    return *this;
+  }
+
+  /// Set all bits to zero
+  SmallBitVec& reset()
+  {
+    m_rack = 0;
+    return *this;
+  }
+
+  /// Return the value of the bit covered by `mask`
+  bool test(const Mask mask) const { return m_rack & mask.m_mask; }
+
+  /// Return the value of the `index`th bit
+  bool test(const size_t index) const { return test(mask(index)); }
+
+  /// Set the bit covered by `mask` to 1
+  SmallBitVec& set(const Mask mask)
+  {
+    m_rack |= mask.m_mask;
+    return *this;
+  }
+
+  /// Set the `index`th bit to 1
+  SmallBitVec& set(const size_t index) { return set(mask(index)); }
+
+  /// Reset the bit covered by `mask` to 0
+  SmallBitVec& reset(const Mask mask)
+  {
+    m_rack &= ~mask.m_mask;
+    return *this;
+  }
+
+  /// Reset the `index`th bit to 0
+  SmallBitVec& reset(const size_t index) { return reset(mask(index)); }
+
+  /// Set the bit covered by `mask` to `value`
+  SmallBitVec& set(const Mask mask, const bool value)
+  {
+    m_rack ^= (-Rack{value} ^ m_rack) & mask.m_mask;
+    return *this;
+  }
+
+  /// Set the `index`th bit to `value`
+  SmallBitVec& set(const size_t index, const bool value)
+  {
+    return set(mask(index), value);
+  }
+
+  /// Flip the value of the bit covered by `mask`
+  SmallBitVec& flip(const Mask mask)
+  {
+    m_rack ^= mask.m_mask;
+    return *this;
+  }
+
+  /// Flip the value of the `index`th bit
+  SmallBitVec& flip(const size_t index) { return flip(mask(index)); }
+
+  bool operator==(const SmallBitVec& vec) const { return m_rack == vec.m_rack; }
+
+  bool operator!=(const SmallBitVec& vec) const { return m_rack != vec.m_rack; }
+
+  bool operator<(const SmallBitVec& vec) const { return m_rack < vec.m_rack; }
+
+  SmallBitVec& operator=(const Rack i)
+  {
+    m_rack = i;
+    return *this;
+  }
+
+  SmallBitVec& operator&=(const SmallBitVec& vec)
+  {
+    m_rack &= vec.m_rack;
+    return *this;
+  }
+
+  SmallBitVec& operator|=(const SmallBitVec& vec)
+  {
+    m_rack |= vec.m_rack;
+    return *this;
+  }
+
+  SmallBitVec& operator^=(const SmallBitVec& vec)
+  {
+    m_rack ^= vec.m_rack;
+    return *this;
+  }
+
+  SmallBitVec& operator^=(const Rack i)
+  {
+    m_rack ^= i;
+    return *this;
+  }
+
+  SmallBitVec& operator<<=(const size_t bits)
+  {
+    assert(bits < size());
+    m_rack <<= bits;
+    return *this;
+  }
+
+  SmallBitVec& operator>>=(const size_t bits)
+  {
+    assert(bits < size());
+    m_rack >>= bits;
+    return *this;
+  }
+
+  /// Right-rotate by `bits` positions
+  SmallBitVec& rotr(const size_t bits)
+  {
+    if (bits > 0 && bits < size()) {
+      assert(bits <= bits_per_rack);
+      m_rack = (m_rack >> bits) | (m_rack << (size() - bits));
+      m_rack &= (~Rack{0} >> (bits_per_rack - size()));
+    }
+    return *this;
+  }
+
+  /// Left-rotate by `bits` positions
+  SmallBitVec& rotl(const size_t bits)
+  {
+    if (bits > 0 && bits < size()) {
+      assert(bits <= bits_per_rack);
+      m_rack = (m_rack << bits) | (m_rack >> (size() - bits));
+      m_rack &= (~Rack{0} >> (bits_per_rack - size()));
+    }
+    return *this;
+  }
+
+  /// Return true iff all bits are zero
+  bool none() const { return m_rack == 0; }
+
+  /// Return 1 + the index of the first set bit, or 0 if there are none
+  size_t find_first() const
+  {
+    return static_cast<size_t>(detail::find_first(m_rack));
+  }
+
+  /// Return the number of set bits
+  size_t count() const
+  {
+    return static_cast<size_t>(detail::pop_count(m_rack));
+  }
+
+  /// Flip all bits (one's complement)
+  SmallBitVec& flip()
+  {
+    m_rack = ~m_rack;
+    return *this;
+  }
+
+  /// Return the first rack
+  Rack& rack() { return m_rack; }
+  Rack  rack() const { return m_rack; }
+
+  /// Return a raw pointer to the racks
+  Rack*       data() { return &m_rack; }
+  const Rack* data() const { return &m_rack; }
+
+  /// Return the number of racks
+  size_t num_racks() const { return 1; }
+
+  template<class BitVec>
+  class iterator_base : public Mask
+  {
+  public:
+    iterator_base& operator++()
+    {
+      Mask::operator++();
+      return *this;
+    }
+
+    iterator_base& operator--()
+    {
+      Mask::operator--();
+      return *this;
+    }
+
+    bool operator==(const iterator_base& rhs) const
+    {
+      return m_vec == rhs.m_vec && Mask::operator==(rhs);
+    }
+
+    bool operator!=(const iterator_base& rhs) const { return !operator==(rhs); }
+
+    bool operator*() const { return m_vec->test(*this); }
+
+  protected:
+    iterator_base(BitVec& vec, const size_t index)
+      : Mask{index}
+      , m_vec{&vec}
+    {}
+
+    BitVec* m_vec;
+  };
+
+  class iterator : public iterator_base<SmallBitVec>
+  {
+  public:
+    void set() { m_vec->set(*this); }
+    void reset() { m_vec->reset(*this); }
+
+  private:
+    friend class SmallBitVec;
+
+    iterator(SmallBitVec& vec, const size_t index)
+      : iterator_base{vec, index}
+    {}
+  };
+
+  class const_iterator : public iterator_base<const SmallBitVec>
+  {
+  private:
+    friend class SmallBitVec;
+
+    const_iterator(const SmallBitVec& vec, const size_t index)
+      : iterator_base{vec, index}
+    {}
+  };
+
+  Mask mask(const size_t i = 0) const
+  {
+    assert(i <= size());
+    return Mask{i};
+  }
+
+  iterator begin(const size_t i = 0) { return {*this, i}; }
+  iterator end() { return {*this, size()}; }
+
+  const_iterator begin(const size_t i = 0) const { return {*this, i}; }
+  const_iterator end() const { return {*this, size()}; }
 
 private:
-	static_assert(8 * sizeof(Rack) == bits_per_rack, "");
-	static_assert((sizeof(Rack) == 4) || (sizeof(Rack) == 8), "");
+  static_assert(8 * sizeof(Rack) == bits_per_rack, "");
+  static_assert((sizeof(Rack) == 4) || (sizeof(Rack) == 8), "");
 
-	Rack   m_rack{};
-	size_t m_size{};
+  Rack   m_rack{};
+  size_t m_size{};
 };
 
 namespace detail {
 
-template <>
-struct is_bitvec<SmallBitVec>
-{
-	constexpr static bool value = true;
+template<>
+struct is_bitvec<SmallBitVec> {
+  constexpr static bool value = true;
 };
 
-template <>
+template<>
 inline void
 gray_code(SmallBitVec& value)
 {
-	value.rack() ^= (value.rack() >> 1U);
+  value.rack() ^= (value.rack() >> 1U);
 }
 
-template <>
+template<>
 inline void
 gray_code_inv(SmallBitVec& value)
 {
-	gray_code_inv<SmallBitVec::Rack>(value.rack());
+  gray_code_inv<SmallBitVec::Rack>(value.rack());
 }
 
 } // namespace detail
diff --git a/include/chilbert/StaticBitVec.hpp b/include/chilbert/StaticBitVec.hpp
index 9aff3ad..1cd59b0 100644
--- a/include/chilbert/StaticBitVec.hpp
+++ b/include/chilbert/StaticBitVec.hpp
@@ -42,75 +42,95 @@ namespace chilbert {
  *
  * @tparam N Number of bits.
  */
-template <size_t N>
+template<size_t N>
 class StaticBitVec : public detail::MultiBitVec<StaticBitVec<N>>
 {
 public:
-	using Rack = typename detail::MultiBitVec<StaticBitVec<N>>::Rack;
+  using Rack = typename detail::MultiBitVec<StaticBitVec<N>>::Rack;
 
-	using detail::MultiBitVec<StaticBitVec<N>>::bits_per_rack;
+  using detail::MultiBitVec<StaticBitVec<N>>::bits_per_rack;
 
-	StaticBitVec() = default;
+  StaticBitVec() = default;
 
-	/// Constructor for compatibility with DynamicBitVec
-	explicit StaticBitVec(const size_t bits) { assert(bits == size()); }
+  /// Constructor for compatibility with DynamicBitVec
+  explicit StaticBitVec(const size_t bits) { assert(bits == size()); }
 
-	/// Constructor for compatibility with DynamicBitVec
-	StaticBitVec(const size_t bits, const Rack value)
-	  : StaticBitVec{bits}
-	{
-		m_racks[0] = value;
-	}
+  /// Constructor for compatibility with DynamicBitVec
+  StaticBitVec(const size_t bits, const Rack value)
+    : StaticBitVec{bits}
+  {
+    m_racks[0] = value;
+  }
 
-	/// Return the size in bits
-	size_t size() const { return N; }
+  /// Return the size in bits
+  size_t size() const { return N; }
 
-	/// Return a reference to the `index`th rack
+  /// Return a reference to the `index`th rack
 #ifndef NDEBUG
-	const auto& rack(const size_t index) const { return m_racks.at(index); }
-	auto&       rack(const size_t index) { return m_racks.at(index); }
+  const auto& rack(const size_t index) const
+  {
+    return m_racks.at(index);
+  }
+  auto& rack(const size_t index)
+  {
+    return m_racks.at(index);
+  }
 #else
-	const auto& rack(const size_t index) const { return m_racks[index]; }
-	auto&       rack(const size_t index) { return m_racks[index]; }
+  const auto& rack(const size_t index) const
+  {
+    return m_racks[index];
+  }
+  auto& rack(const size_t index)
+  {
+    return m_racks[index];
+  }
 #endif
 
-	/// Return a raw pointer to the racks
-	Rack*       data() { return m_racks.data(); }
-	const Rack* data() const { return m_racks.data(); }
-
-	/// Return the total size of all racks in bytes
-	static constexpr size_t data_size() { return num_racks() * sizeof(Rack); }
-
-	/// Return the number of racks
-	static constexpr size_t num_racks()
-	{
-		return (std::max(N, size_t(1)) + bits_per_rack - 1) / bits_per_rack;
-	}
+  /// Return a raw pointer to the racks
+  Rack* data()
+  {
+    return m_racks.data();
+  }
+  const Rack* data() const
+  {
+    return m_racks.data();
+  }
+
+  /// Return the total size of all racks in bytes
+  static constexpr size_t data_size()
+  {
+    return num_racks() * sizeof(Rack);
+  }
+
+  /// Return the number of racks
+  static constexpr size_t num_racks()
+  {
+    return (std::max(N, size_t(1)) + bits_per_rack - 1) / bits_per_rack;
+  }
 
 private:
-	std::array<Rack, num_racks()> m_racks{};
+  std::array<Rack, num_racks()> m_racks{};
 };
 
 namespace detail {
 
-template <size_t N>
-struct is_bitvec<StaticBitVec<N>>
-{
-	constexpr static bool value = true;
+template<size_t N>
+struct is_bitvec<StaticBitVec<N>> {
+  constexpr static bool value = true;
 };
 
-template <size_t N>
+template<size_t N>
 void
 gray_code(StaticBitVec<N>& value)
 {
-	gray_code(static_cast<MultiBitVec<StaticBitVec<N>>&>(value));
+  gray_code(static_cast<MultiBitVec<StaticBitVec<N>>&>(value));
 }
 
-template <size_t N>
+template<size_t N>
 void
 gray_code_inv(StaticBitVec<N>& value)
 {
-	gray_code_inv(static_cast<MultiBitVec<StaticBitVec<N>>&>(value));
+  gray_code_inv(static_cast<MultiBitVec<StaticBitVec<N>>&>(value));
 }
 
 } // namespace detail
diff --git a/include/chilbert/chilbert.hpp b/include/chilbert/chilbert.hpp
index c55fc14..36b94cc 100644
--- a/include/chilbert/chilbert.hpp
+++ b/include/chilbert/chilbert.hpp
@@ -33,8 +33,9 @@ namespace chilbert {
  * @param n Number of dimensions.
  * @param[out] h Hilbert Index.
  */
-template <class P, class H>
-inline void coords_to_index(const P& p, size_t m, size_t n, H& h);
+template<class P, class H>
+inline void
+coords_to_index(const P& p, size_t m, size_t n, H& h);
 
 /** Map the Hilbert Index `p` to a point.
  *
@@ -46,8 +47,9 @@ inline void coords_to_index(const P& p, size_t m, size_t n, H& h);
  * @param n Number of dimensions.
  * @param h Hilbert Index.
  */
-template <class P, class H>
-inline void index_to_coords(P& p, size_t m, size_t n, const H& h);
+template<class P, class H>
+inline void
+index_to_coords(P& p, size_t m, size_t n, const H& h);
 
 /** Map the point `p` to a Compact Hilbert Index.
  *
@@ -61,13 +63,14 @@ inline void index_to_coords(P& p, size_t m, size_t n, const H& h);
  * @param M Optional net precision (sum of `ms`), the size of `hc` in bits.
  * @param m Optional largest precision in `m`.
  */
-template <class P, class H>
-inline void coords_to_compact_index(const P&      p,
-                                    const size_t* ms,
-                                    size_t        n,
-                                    H&            hc,
-                                    size_t        M = 0,
-                                    size_t        m = 0);
+template<class P, class H>
+inline void
+coords_to_compact_index(const P&      p,
+                        const size_t* ms,
+                        size_t        n,
+                        H&            hc,
+                        size_t        M = 0,
+                        size_t        m = 0);
 
 /** Map the Compact Hilbert Index `hc` to a point.
  *
@@ -81,13 +84,14 @@ inline void coords_to_compact_index(const P&      p,
  * @param M Optional net precision (sum of `ms`), the size of `hc` in bits.
  * @param m Optional largest precision in `m`.
  */
-template <class P, class H>
-inline void compact_index_to_coords(P&            p,
-                                    const size_t* ms,
-                                    size_t        n,
-                                    const H&      hc,
-                                    size_t        M = 0,
-                                    size_t        m = 0);
+template<class P, class H>
+inline void
+compact_index_to_coords(P&            p,
+                        const size_t* ms,
+                        size_t        n,
+                        const H&      hc,
+                        size_t        M = 0,
+                        size_t        m = 0);
 
 } // namespace chilbert
 
diff --git a/include/chilbert/chilbert.ipp b/include/chilbert/chilbert.ipp
index 5c18570..f333630 100644
--- a/include/chilbert/chilbert.ipp
+++ b/include/chilbert/chilbert.ipp
@@ -41,34 +41,34 @@ namespace detail {
  */
 static constexpr size_t D0 = 1;
 
-template <class T>
+template<class T>
 size_t
 num_bits(const T&, std::enable_if_t<std::is_integral<T>::value>* = nullptr)
 {
-	return sizeof(T) * CHAR_BIT;
+  return sizeof(T) * CHAR_BIT;
 }
 
-template <class T>
+template<class T>
 size_t
 num_bits(const T& vec,
          std::enable_if_t<std::is_same<T, SmallBitVec>::value ||
                           std::is_same<T, DynamicBitVec>::value>* = nullptr)
 {
-	return vec.size();
+  return vec.size();
 }
 
-template <size_t N>
+template<size_t N>
 size_t
 num_bits(const StaticBitVec<N>&, void* = nullptr)
 {
-	return N;
+  return N;
 }
 
-template <size_t MaxN>
+template<size_t MaxN>
 size_t
 num_bits(const BoundedBitVec<MaxN>& vec, void* = nullptr)
 {
-	return vec.size();
+  return vec.size();
 }
 
 /** Copy a range of bits from one field to the start of another.
@@ -78,13 +78,13 @@ num_bits(const BoundedBitVec<MaxN>& vec, void* = nullptr)
  * @param i Start bit index in source
  * @param[out] w Destination bit field
  */
-template <class H, class I>
+template<class H, class I>
 inline void
 get_bits(const H& h, const size_t n, const size_t i, I& w)
 {
-	for (size_t j = 0; j < n; ++j) {
-		set_bit(w, j, test_bit(h, i + j));
-	}
+  for (size_t j = 0; j < n; ++j) {
+    set_bit(w, j, test_bit(h, i + j));
+  }
 }
 
 /** Set a range of bits in one field from the start of another.
@@ -94,13 +94,13 @@ get_bits(const H& h, const size_t n, const size_t i, I& w)
  * @param i Start bit index in destination
  * @param w Source bit field
  */
-template <class H, class I>
+template<class H, class I>
 inline void
 set_bits(H& h, const size_t n, const size_t i, const I& w)
 {
-	for (size_t j = 0; j < n; ++j) {
-		set_bit(h, i + j, test_bit(w, j));
-	}
+  for (size_t j = 0; j < n; ++j) {
+    set_bit(h, i + j, test_bit(w, j));
+  }
 }
 
 /** Set the leading bits in `l` to bit `i` from the corresponding value in `p`.
@@ -110,13 +110,13 @@ set_bits(H& h, const size_t n, const size_t i, const I& w)
  * @param i Bit position to copy from values in `p`.
  * @param[out] l Output bits.
  */
-template <class P, class I>
+template<class P, class I>
 inline void
 get_location(const P& p, const size_t n, const size_t i, I& l)
 {
-	for (size_t j = 0; j < n; ++j) {
-		set_bit(l, j, test_bit(p[j], i));
-	}
+  for (size_t j = 0; j < n; ++j) {
+    set_bit(l, j, test_bit(p[j], i));
+  }
 }
 
 /** Set bit `i` in values in `p` to the corresponding leadings bits in `l`.
@@ -126,81 +126,81 @@ get_location(const P& p, const size_t n, const size_t i, I& l)
  * @param i Bit position to set in values in `p`.
  * @param l Input bits.
  */
-template <class P, class I>
+template<class P, class I>
 inline void
 set_location(P& p, const size_t n, const size_t i, const I& l)
 {
-	for (size_t j = 0; j < n; ++j) {
-		set_bit(p[j], i, test_bit(l, j));
-	}
+  for (size_t j = 0; j < n; ++j) {
+    set_bit(p[j], i, test_bit(l, j));
+  }
 }
 
 // 'Transforms' a point.
-template <class I>
+template<class I>
 inline void
 transform(const I& e, const size_t d, const size_t n, I& a)
 {
-	(void)n;
-	assert(a.size() == n);
-	a ^= e;
-	a.rotr(d);
+  (void)n;
+  assert(a.size() == n);
+  a ^= e;
+  a.rotr(d);
 }
 
 // Inverse 'transforms' a point.
-template <class I>
+template<class I>
 inline void
 transform_inv(const I& e, const size_t d, const size_t n, I& a)
 {
-	assert(a.size() == n);
-	a.rotl(d);
-	a ^= e;
+  assert(a.size() == n);
+  a.rotl(d);
+  a ^= e;
 }
 
 // Update for method 1 (GrayCodeInv in the loop)
-template <class I>
+template<class I>
 inline void
 update1(const I& l, const I& t, const I& w, const size_t n, I& e, size_t& d)
 {
-	assert(d < n);
-	e = l;
-	e.flip(d); //#D d == n-1 ? 0 : d+1 );
-
-	// Update direction
-	d += 1 + t.find_first();
-	if (d >= n) {
-		d -= n;
-	}
-	if (d >= n) {
-		d -= n;
-	}
-	assert(d < n);
-
-	if (!w.test(0)) {
-		e.flip(d == 0 ? n - 1 : d - 1); //#D d );
-	}
+  assert(d < n);
+  e = l;
+  e.flip(d); //#D d == n-1 ? 0 : d+1 );
+
+  // Update direction
+  d += 1 + t.find_first();
+  if (d >= n) {
+    d -= n;
+  }
+  if (d >= n) {
+    d -= n;
+  }
+  assert(d < n);
+
+  if (!w.test(0)) {
+    e.flip(d == 0 ? n - 1 : d - 1); //#D d );
+  }
 }
 
 // Update for method 2 (GrayCodeInv out of loop)
-template <class I>
+template<class I>
 inline void
 update2(const I& l, const I& t, const size_t n, I& e, size_t& d)
 {
-	assert(d < n);
-	e = l;
-	e.flip(d); //#D d == n-1 ? 0 : d+1 );
-
-	// Update direction
-	d += 1 + t.find_first();
-	if (d >= n) {
-		d -= n;
-	}
-	if (d >= n) {
-		d -= n;
-	}
-	assert(d < n);
+  assert(d < n);
+  e = l;
+  e.flip(d); //#D d == n-1 ? 0 : d+1 );
+
+  // Update direction
+  d += 1 + t.find_first();
+  if (d >= n) {
+    d -= n;
+  }
+  if (d >= n) {
+    d -= n;
+  }
+  assert(d < n);
 }
 
-template <class P, class H, class I>
+template<class P, class H, class I>
 inline void
 coords_to_index(const P&      p,
                 const size_t  m,
@@ -209,91 +209,91 @@ coords_to_index(const P&      p,
                 I&&           scratch,
                 size_t* const ds = nullptr)
 {
-	I e{std::forward<I>(scratch)};
-	I l{e};
-	I t{e};
-	I w{e};
-
-	// Initialize
-	e.reset();
-	l.reset();
-	reset_bits(h);
-
-	// Work from MSB to LSB
-	size_t d  = D0;
-	size_t ho = m * n;
-	for (auto i = static_cast<intptr_t>(m - 1); i >= 0; --i) {
-		if (ds) {
-			ds[i] = d;
-		}
-
-		// Get corner of sub-hypercube where point lies.
-		get_location<P, I>(p, n, static_cast<size_t>(i), l);
-
-		// Mirror and reflect the location.
-		// t = T_{(e,d)}(l)
-		t = l;
-		transform<I>(e, d, n, t);
-
-		w = t;
-		if (static_cast<size_t>(i) < m - 1) {
-			w.flip(n - 1);
-		}
-
-		// Concatenate to the index.
-		ho -= n;
-		set_bits<H, I>(h, n, ho, w);
-
-		// Update the entry point and direction.
-		update2<I>(l, t, n, e, d);
-	}
-
-	gray_code_inv(h);
+  I e{std::forward<I>(scratch)};
+  I l{e};
+  I t{e};
+  I w{e};
+
+  // Initialize
+  e.reset();
+  l.reset();
+  reset_bits(h);
+
+  // Work from MSB to LSB
+  size_t d  = D0;
+  size_t ho = m * n;
+  for (auto i = static_cast<intptr_t>(m - 1); i >= 0; --i) {
+    if (ds) {
+      ds[i] = d;
+    }
+
+    // Get corner of sub-hypercube where point lies.
+    get_location<P, I>(p, n, static_cast<size_t>(i), l);
+
+    // Mirror and reflect the location.
+    // t = T_{(e,d)}(l)
+    t = l;
+    transform<I>(e, d, n, t);
+
+    w = t;
+    if (static_cast<size_t>(i) < m - 1) {
+      w.flip(n - 1);
+    }
+
+    // Concatenate to the index.
+    ho -= n;
+    set_bits<H, I>(h, n, ho, w);
+
+    // Update the entry point and direction.
+    update2<I>(l, t, n, e, d);
+  }
+
+  gray_code_inv(h);
 }
 
-template <class P, class H, class I>
+template<class P, class H, class I>
 inline void
 index_to_coords(P& p, const size_t m, const size_t n, const H& h, I&& scratch)
 {
-	I e{std::forward<I>(scratch)};
-	I l{e};
-	I t{e};
-	I w{e};
-
-	// Initialize
-	e.reset();
-	l.reset();
-	for (size_t j = 0; j < n; ++j) {
-		reset_bits(p[j]);
-	}
-
-	// Work from MSB to LSB
-	size_t d  = D0;
-	size_t ho = m * n;
-	for (auto i = static_cast<intptr_t>(m - 1); i >= 0; --i) {
-		// Get the Hilbert index bits
-		ho -= n;
-		get_bits<H, I>(h, n, ho, w);
-
-		// t = GrayCode(w)
-		t = w;
-		gray_code(t);
-
-		// Reverse the transform
-		// l = T^{-1}_{(e,d)}(t)
-		l = t;
-		transform_inv<I>(e, d, n, l);
-
-		// Distribute these bits
-		// to the coordinates.
-		set_location<P, I>(p, n, static_cast<size_t>(i), l);
-
-		// Update the entry point and direction.
-		update1<I>(l, t, w, n, e, d);
-	}
+  I e{std::forward<I>(scratch)};
+  I l{e};
+  I t{e};
+  I w{e};
+
+  // Initialize
+  e.reset();
+  l.reset();
+  for (size_t j = 0; j < n; ++j) {
+    reset_bits(p[j]);
+  }
+
+  // Work from MSB to LSB
+  size_t d  = D0;
+  size_t ho = m * n;
+  for (auto i = static_cast<intptr_t>(m - 1); i >= 0; --i) {
+    // Get the Hilbert index bits
+    ho -= n;
+    get_bits<H, I>(h, n, ho, w);
+
+    // t = GrayCode(w)
+    t = w;
+    gray_code(t);
+
+    // Reverse the transform
+    // l = T^{-1}_{(e,d)}(t)
+    l = t;
+    transform_inv<I>(e, d, n, l);
+
+    // Distribute these bits
+    // to the coordinates.
+    set_location<P, I>(p, n, static_cast<size_t>(i), l);
+
+    // Update the entry point and direction.
+    update1<I>(l, t, w, n, e, d);
+  }
 }
 
-template <class P, class HC, class I>
+template<class P, class HC, class I>
 inline void
 coords_to_compact_index(const P&            p,
                         const size_t* const ms,
@@ -303,43 +303,43 @@ coords_to_compact_index(const P&            p,
                         size_t              M = 0,
                         size_t              m = 0)
 {
-	// Get total precision and max precision if not supplied
-	if (M == 0 || m == 0) {
-		M = m = 0;
-		for (size_t i = 0; i < n; ++i) {
-			assert(num_bits(p[i]) >= ms[i]);
-			if (ms[i] > m) {
-				m = ms[i];
-			}
-			M += ms[i];
-		}
-	}
-
-	const size_t mn = m * n;
-
-	assert(num_bits(hc) >= M);
-
-	// If we could avoid allocation altogether (ie: have a
-	// fixed buffer allocated on the stack) then this increases
-	// speed by a bit (4% when n=4, m=20)
-	auto* const ds = new size_t[m];
-
-	if (mn > SmallBitVec::bits_per_rack) {
-		DynamicBitVec h(mn);
-		detail::coords_to_index<P, DynamicBitVec, I>(
-		  p, m, n, h, std::forward<I>(scratch), ds);
-		compact_index<DynamicBitVec, HC>(ms, ds, n, m, h, hc);
-	} else {
-		SmallBitVec h(mn);
-		detail::coords_to_index<P, SmallBitVec, I>(
-		  p, m, n, h, std::forward<I>(scratch), ds);
-		compact_index<SmallBitVec, HC>(ms, ds, n, m, h, hc);
-	}
-
-	delete[] ds;
+  // Get total precision and max precision if not supplied
+  if (M == 0 || m == 0) {
+    M = m = 0;
+    for (size_t i = 0; i < n; ++i) {
+      assert(num_bits(p[i]) >= ms[i]);
+      if (ms[i] > m) {
+        m = ms[i];
+      }
+      M += ms[i];
+    }
+  }
+
+  const size_t mn = m * n;
+
+  assert(num_bits(hc) >= M);
+
+  // If we could avoid allocation altogether (ie: have a
+  // fixed buffer allocated on the stack) then this increases
+  // speed by a bit (4% when n=4, m=20)
+  auto* const ds = new size_t[m];
+
+  if (mn > SmallBitVec::bits_per_rack) {
+    DynamicBitVec h(mn);
+    detail::coords_to_index<P, DynamicBitVec, I>(
+      p, m, n, h, std::forward<I>(scratch), ds);
+    compact_index<DynamicBitVec, HC>(ms, ds, n, m, h, hc);
+  } else {
+    SmallBitVec h(mn);
+    detail::coords_to_index<P, SmallBitVec, I>(
+      p, m, n, h, std::forward<I>(scratch), ds);
+    compact_index<SmallBitVec, HC>(ms, ds, n, m, h, hc);
+  }
+
+  delete[] ds;
 }
 
-template <class P, class HC, class I>
+template<class P, class HC, class I>
 inline void
 compact_index_to_coords(P&            p,
                         const size_t* ms,
@@ -349,109 +349,107 @@ compact_index_to_coords(P&            p,
                         size_t        M,
                         size_t        m)
 {
-	I e{std::forward<I>(scratch)};
-	I l{e};
-	I t{e};
-	I w{e};
-	I r{e};
-	I mask{e};
-	I ptrn{e};
-
-	// Get total precision and max precision
-	// if not supplied
-	if (M == 0 || m == 0) {
-		M = m = 0;
-		for (size_t i = 0; i < n; ++i) {
-			if (ms[i] > m) {
-				m = ms[i];
-			}
-			M += ms[i];
-		}
-	}
-
-	assert(num_bits(hc) >= M);
-	assert(num_bits(p[0]) >= m);
-
-	// Initialize
-	e.reset();
-	l.reset();
-	for (size_t j = 0; j < n; ++j) {
-		reset_bits(p[j]);
-	}
-
-	// Work from MSB to LSB
-	size_t d = D0;
-
-	for (auto i = static_cast<intptr_t>(m - 1); i >= 0; --i) {
-		// Get the mask and ptrn
-		size_t b = 0;
-		extract_mask<I>(ms, n, d, static_cast<size_t>(i), mask, b);
-		ptrn = e;
-		assert(ptrn.size() == n);
-		ptrn.rotr(d);
-
-		// Get the Hilbert index bits
-		M -= b;
-		r.reset(); // GetBits doesn't do this
-		get_bits<HC, I>(hc, b, M, r);
-
-		// w = GrayCodeRankInv(r)
-		// t = GrayCode(w)
-		gray_code_rank_inv<I>(mask, ptrn, r, n, b, t, w);
-
-		// Reverse the transform
-		// l = T^{-1}_{(e,d)}(t)
-		l = t;
-		transform_inv<I>(e, d, n, l);
-
-		// Distribute these bits
-		// to the coordinates.
-		set_location<P, I>(p, n, static_cast<size_t>(i), l);
-
-		// Update the entry point and direction.
-		update1<I>(l, t, w, n, e, d);
-	}
+  I e{std::forward<I>(scratch)};
+  I l{e};
+  I t{e};
+  I w{e};
+  I r{e};
+  I mask{e};
+  I ptrn{e};
+
+  // Get total precision and max precision
+  // if not supplied
+  if (M == 0 || m == 0) {
+    M = m = 0;
+    for (size_t i = 0; i < n; ++i) {
+      if (ms[i] > m) {
+        m = ms[i];
+      }
+      M += ms[i];
+    }
+  }
+
+  assert(num_bits(hc) >= M);
+  assert(num_bits(p[0]) >= m);
+
+  // Initialize
+  e.reset();
+  l.reset();
+  for (size_t j = 0; j < n; ++j) {
+    reset_bits(p[j]);
+  }
+
+  // Work from MSB to LSB
+  size_t d = D0;
+
+  for (auto i = static_cast<intptr_t>(m - 1); i >= 0; --i) {
+    // Get the mask and ptrn
+    size_t b = 0;
+    extract_mask<I>(ms, n, d, static_cast<size_t>(i), mask, b);
+    ptrn = e;
+    assert(ptrn.size() == n);
+    ptrn.rotr(d);
+
+    // Get the Hilbert index bits
+    M -= b;
+    r.reset(); // GetBits doesn't do this
+    get_bits<HC, I>(hc, b, M, r);
+
+    // w = GrayCodeRankInv(r)
+    // t = GrayCode(w)
+    gray_code_rank_inv<I>(mask, ptrn, r, n, b, t, w);
+
+    // Reverse the transform
+    // l = T^{-1}_{(e,d)}(t)
+    l = t;
+    transform_inv<I>(e, d, n, l);
+
+    // Distribute these bits
+    // to the coordinates.
+    set_location<P, I>(p, n, static_cast<size_t>(i), l);
+
+    // Update the entry point and direction.
+    update1<I>(l, t, w, n, e, d);
+  }
 }
 
 } // namespace detail
 
-template <class P, class H>
+template<class P, class H>
 inline void
 coords_to_index(const P& p, const size_t m, const size_t n, H& h)
 {
-	assert(detail::num_bits(h) >= n * m);
-	assert(detail::num_bits(p[0]) >= m);
-
-	if (n <= SmallBitVec::bits_per_rack) {
-		// Intermediate variables will fit in fixed width
-		detail::coords_to_index<P, H, SmallBitVec>(p, m, n, h, SmallBitVec(n));
-	} else {
-		// Otherwise, they must be DynamicBitVecs
-		detail::coords_to_index<P, H, DynamicBitVec>(
-		  p, m, n, h, DynamicBitVec(n));
-	}
+  assert(detail::num_bits(h) >= n * m);
+  assert(detail::num_bits(p[0]) >= m);
+
+  if (n <= SmallBitVec::bits_per_rack) {
+    // Intermediate variables will fit in fixed width
+    detail::coords_to_index<P, H, SmallBitVec>(p, m, n, h, SmallBitVec(n));
+  } else {
+    // Otherwise, they must be DynamicBitVecs
+    detail::coords_to_index<P, H, DynamicBitVec>(p, m, n, h, DynamicBitVec(n));
+  }
 }
 
-template <class P, class H>
+template<class P, class H>
 inline void
 index_to_coords(P& p, const size_t m, const size_t n, const H& h)
 {
-	assert(m > 0);
-	assert(n > 0);
-	assert(detail::num_bits(h) >= n * m);
-	assert(detail::num_bits(p[0]) >= m);
-
-	if (n <= SmallBitVec::bits_per_rack) {
-		// Intermediate variables will fit in fixed width
-		detail::index_to_coords<P, H, SmallBitVec>(p, m, n, h, SmallBitVec(n));
-	} else {
-		// Otherwise, they must be DynamicBitVecs
-		detail::index_to_coords<P, H, DynamicBitVec>(
-		  p, m, n, h, DynamicBitVec(n));
-	}
+  assert(m > 0);
+  assert(n > 0);
+  assert(detail::num_bits(h) >= n * m);
+  assert(detail::num_bits(p[0]) >= m);
+
+  if (n <= SmallBitVec::bits_per_rack) {
+    // Intermediate variables will fit in fixed width
+    detail::index_to_coords<P, H, SmallBitVec>(p, m, n, h, SmallBitVec(n));
+  } else {
+    // Otherwise, they must be DynamicBitVecs
+    detail::index_to_coords<P, H, DynamicBitVec>(p, m, n, h, DynamicBitVec(n));
+  }
 }
 
-template <class P, class HC>
+template<class P, class HC>
 inline void
 coords_to_compact_index(const P&            p,
                         const size_t* const ms,
@@ -460,20 +458,20 @@ coords_to_compact_index(const P&            p,
                         const size_t        M,
                         const size_t        m)
 {
-	assert(hc.size() >= std::accumulate(ms, ms + n, size_t(0)));
-
-	if (n <= SmallBitVec::bits_per_rack) {
-		// Intermediate variables will fit in fixed width
-		detail::coords_to_compact_index<P, HC, SmallBitVec>(
-		  p, ms, n, hc, SmallBitVec(n), M, m);
-	} else {
-		// Otherwise, they must be DynamicBitVecs
-		detail::coords_to_compact_index<P, HC, DynamicBitVec>(
-		  p, ms, n, hc, DynamicBitVec(n), M, m);
-	}
+  assert(hc.size() >= std::accumulate(ms, ms + n, size_t(0)));
+
+  if (n <= SmallBitVec::bits_per_rack) {
+    // Intermediate variables will fit in fixed width
+    detail::coords_to_compact_index<P, HC, SmallBitVec>(
+      p, ms, n, hc, SmallBitVec(n), M, m);
+  } else {
+    // Otherwise, they must be DynamicBitVecs
+    detail::coords_to_compact_index<P, HC, DynamicBitVec>(
+      p, ms, n, hc, DynamicBitVec(n), M, m);
+  }
 }
 
-template <class P, class HC>
+template<class P, class HC>
 inline void
 compact_index_to_coords(P&                  p,
                         const size_t* const ms,
@@ -482,19 +480,19 @@ compact_index_to_coords(P&                  p,
                         const size_t        M,
                         const size_t        m)
 {
-	assert(hc.size() >= std::accumulate(ms, ms + n, size_t(0)));
-
-	if (n <= SmallBitVec::bits_per_rack) {
-		// Intermediate variables will fit in fixed width
-		SmallBitVec scratch(n);
-		detail::compact_index_to_coords<P, HC, SmallBitVec>(
-		  p, ms, n, hc, std::move(scratch), M, m);
-	} else {
-		// Otherwise, they must be DynamicBitVecs
-		DynamicBitVec scratch(n);
-		detail::compact_index_to_coords<P, HC, DynamicBitVec>(
-		  p, ms, n, hc, std::move(scratch), M, m);
-	}
+  assert(hc.size() >= std::accumulate(ms, ms + n, size_t(0)));
+
+  if (n <= SmallBitVec::bits_per_rack) {
+    // Intermediate variables will fit in fixed width
+    SmallBitVec scratch(n);
+    detail::compact_index_to_coords<P, HC, SmallBitVec>(
+      p, ms, n, hc, std::move(scratch), M, m);
+  } else {
+    // Otherwise, they must be DynamicBitVecs
+    DynamicBitVec scratch(n);
+    detail::compact_index_to_coords<P, HC, DynamicBitVec>(
+      p, ms, n, hc, std::move(scratch), M, m);
+  }
 }
 
 } // namespace chilbert
diff --git a/include/chilbert/detail/BitVecIndex.hpp b/include/chilbert/detail/BitVecIndex.hpp
index c8d0469..c337c70 100644
--- a/include/chilbert/detail/BitVecIndex.hpp
+++ b/include/chilbert/detail/BitVecIndex.hpp
@@ -27,22 +27,21 @@ namespace chilbert {
 namespace detail {
 
 /// Index into a multi-rack bit vector
-template <class BitVec>
-struct BitVecIndex
-{
-	using Rack = typename BitVec::Rack;
-
-	static constexpr size_t bits_per_rack = sizeof(Rack) * CHAR_BIT;
-
-	explicit BitVecIndex(const size_t bits)
-	  : rack{bits / bits_per_rack}
-	  , bit{bits - rack * bits_per_rack}
-	{
-		assert(bit < bits_per_rack);
-	}
-
-	size_t rack;
-	size_t bit;
+template<class BitVec>
+struct BitVecIndex {
+  using Rack = typename BitVec::Rack;
+
+  static constexpr size_t bits_per_rack = sizeof(Rack) * CHAR_BIT;
+
+  explicit BitVecIndex(const size_t bits)
+    : rack{bits / bits_per_rack}
+    , bit{bits - rack * bits_per_rack}
+  {
+    assert(bit < bits_per_rack);
+  }
+
+  size_t rack;
+  size_t bit;
 };
 
 } // namespace detail
diff --git a/include/chilbert/detail/BitVecIterator.hpp b/include/chilbert/detail/BitVecIterator.hpp
index 2458da9..43af2ee 100644
--- a/include/chilbert/detail/BitVecIterator.hpp
+++ b/include/chilbert/detail/BitVecIterator.hpp
@@ -26,72 +26,69 @@
 namespace chilbert {
 namespace detail {
 
-template <class BitVec>
+template<class BitVec>
 class BitVecIteratorBase : public BitVecMask<typename BitVec::Rack>
 {
 public:
-	using Mask = typename BitVec::Mask;
+  using Mask = typename BitVec::Mask;
 
-	BitVecIteratorBase& operator++()
-	{
-		Mask::operator++();
-		return *this;
-	}
+  BitVecIteratorBase& operator++()
+  {
+    Mask::operator++();
+    return *this;
+  }
 
-	BitVecIteratorBase& operator--()
-	{
-		Mask::operator--();
-		return *this;
-	}
+  BitVecIteratorBase& operator--()
+  {
+    Mask::operator--();
+    return *this;
+  }
 
-	bool operator==(const BitVecIteratorBase& rhs) const
-	{
-		return m_vec == rhs.m_vec && Mask::operator==(rhs);
-	}
+  bool operator==(const BitVecIteratorBase& rhs) const
+  {
+    return m_vec == rhs.m_vec && Mask::operator==(rhs);
+  }
 
-	bool operator!=(const BitVecIteratorBase& rhs) const
-	{
-		return !operator==(rhs);
-	}
+  bool operator!=(const BitVecIteratorBase& rhs) const
+  {
+    return !operator==(rhs);
+  }
 
-	bool operator*() const { return m_vec->test(*this); }
+  bool operator*() const { return m_vec->test(*this); }
 
 protected:
-	BitVecIteratorBase(BitVec& vec, const size_t index)
-	  : Mask{index}
-	  , m_vec{&vec}
-	{
-	}
+  BitVecIteratorBase(BitVec& vec, const size_t index)
+    : Mask{index}
+    , m_vec{&vec}
+  {}
 
-	BitVec* m_vec;
+  BitVec* m_vec;
 };
 
-template <class BitVec>
+template<class BitVec>
 class BitVecIterator : public BitVecIteratorBase<BitVec>
 {
 public:
-	void set() { this->m_vec->set(*this); }
-	void reset() { this->m_vec->reset(*this); }
+  void set() { this->m_vec->set(*this); }
+  void reset() { this->m_vec->reset(*this); }
 
 private:
-	friend BitVec;
+  friend BitVec;
 
-	BitVecIterator(BitVec& vec, const size_t index)
-	  : BitVecIteratorBase<BitVec>{vec, index}
-	{
-	}
+  BitVecIterator(BitVec& vec, const size_t index)
+    : BitVecIteratorBase<BitVec>{vec, index}
+  {}
 };
 
-template <class BitVec>
+template<class BitVec>
 class ConstBitVecIterator : public BitVecIteratorBase<const BitVec>
 {
 private:
-	friend BitVec;
+  friend BitVec;
 
-	ConstBitVecIterator(const BitVec& vec, const size_t index)
-	  : BitVecIteratorBase<const BitVec>{vec, index}
-	{
-	}
+  ConstBitVecIterator(const BitVec& vec, const size_t index)
+    : BitVecIteratorBase<const BitVec>{vec, index}
+  {}
 };
 
 } // namespace detail
diff --git a/include/chilbert/detail/BitVecMask.hpp b/include/chilbert/detail/BitVecMask.hpp
index 5674114..855e4f3 100644
--- a/include/chilbert/detail/BitVecMask.hpp
+++ b/include/chilbert/detail/BitVecMask.hpp
@@ -30,42 +30,40 @@ namespace detail {
  * This enables fast iteration while setting or resetting bits since it is
  * internally stored as a mask which avoids repeated index math.
  */
-template <class Rack>
-struct BitVecMask
-{
-	static constexpr size_t bits_per_rack = sizeof(Rack) * CHAR_BIT;
+template<class Rack>
+struct BitVecMask {
+  static constexpr size_t bits_per_rack = sizeof(Rack) * CHAR_BIT;
 
-	explicit BitVecMask(const size_t index)
-	  : rack{index / bits_per_rack}
-	  , mask{Rack{1} << (index - rack * bits_per_rack)}
-	{
-	}
+  explicit BitVecMask(const size_t index)
+    : rack{index / bits_per_rack}
+    , mask{Rack{1} << (index - rack * bits_per_rack)}
+  {}
 
-	void operator++()
-	{
-		if ((mask <<= 1U) == 0U) {
-			mask = 1U;
-			++rack;
-		}
-	}
+  void operator++()
+  {
+    if ((mask <<= 1U) == 0U) {
+      mask = 1U;
+      ++rack;
+    }
+  }
 
-	void operator--()
-	{
-		if ((mask >>= 1U) == 0U) {
-			mask = Rack{1} << (bits_per_rack - 1U);
-			--rack;
-		}
-	}
+  void operator--()
+  {
+    if ((mask >>= 1U) == 0U) {
+      mask = Rack{1} << (bits_per_rack - 1U);
+      --rack;
+    }
+  }
 
-	bool operator==(const BitVecMask& rhs) const
-	{
-		return rack == rhs.rack && mask == rhs.mask;
-	}
+  bool operator==(const BitVecMask& rhs) const
+  {
+    return rack == rhs.rack && mask == rhs.mask;
+  }
 
-	bool operator!=(const BitVecMask& rhs) const { return !operator==(rhs); }
+  bool operator!=(const BitVecMask& rhs) const { return !operator==(rhs); }
 
-	size_t rack;
-	Rack   mask;
+  size_t rack;
+  Rack   mask;
 };
 
 } // namespace detail
diff --git a/include/chilbert/detail/MultiBitVec.hpp b/include/chilbert/detail/MultiBitVec.hpp
index 6db1065..e8ae08a 100644
--- a/include/chilbert/detail/MultiBitVec.hpp
+++ b/include/chilbert/detail/MultiBitVec.hpp
@@ -33,358 +33,358 @@
 namespace chilbert {
 namespace detail {
 
-template <class Derived>
+template<class Derived>
 class MultiBitVec
 {
 public:
-	using Rack = uintptr_t;
-	using Mask = BitVecMask<Rack>;
-
-	using iterator       = BitVecIterator<MultiBitVec<Derived>>;
-	using const_iterator = ConstBitVecIterator<MultiBitVec<Derived>>;
-
-	static constexpr size_t bits_per_rack = sizeof(Rack) * CHAR_BIT;
-
-	/// Return the value of the bit covered by `mask`
-	bool test(const Mask mask) const { return rack(mask.rack) & mask.mask; }
-
-	/// Return the value of the `index`th bit
-	bool test(const size_t index) const { return test(mask(index)); }
-
-	/// Set all bits to one
-	Derived& set()
-	{
-		if (size()) {
-			memset(data(), 0xFF, data_size());
-			self()->truncate();
-		}
-
-		return *self();
-	}
-
-	/// Set the bit covered by `mask` to 1
-	Derived& set(const Mask mask)
-	{
-		rack(mask.rack) |= mask.mask;
-		return *self();
-	}
-
-	/// Set the `index`th bit to 1
-	Derived& set(const size_t index) { return set(mask(index)); }
-
-	/// Set the bit covered by `mask` to `value`
-	Derived& set(const Mask mask, const bool value)
-	{
-		auto& r = rack(mask.rack);
-		r ^= (-Rack{value} ^ r) & mask.mask;
-		return *self();
-	}
-
-	/// Set the `index`th bit to `value`
-	Derived& set(const size_t index, const bool value)
-	{
-		return set(mask(index), value);
-	}
-
-	/// Set all bits to zero
-	Derived& reset()
-	{
-		memset(data(), 0, data_size());
-		return *self();
-	}
-
-	/// Reset the bit covered by `mask` to 0
-	Derived& reset(const Mask mask)
-	{
-		rack(mask.rack) &= ~mask.mask;
-		return *self();
-	}
-
-	/// Reset the `index`th bit to 0
-	Derived& reset(const size_t index) { return reset(mask(index)); }
-
-	/// Flip all bits (one's complement)
-	Derived& flip()
-	{
-		for (size_t i = 0; i < num_racks(); ++i) {
-			rack(i) = ~rack(i);
-		}
-		return *self();
-	}
-
-	/// Flip the value of the bit covered by `mask`
-	Derived& flip(const Mask mask)
-	{
-		rack(mask.rack) ^= mask.mask;
-		return *self();
-	}
-
-	/// Flip the value of the `index`th bit
-	Derived& flip(const size_t index) { return flip(mask(index)); }
-
-	/// Clear any bits in storage outside the valid range if necessary
-	void truncate()
-	{
-		if (const auto pad = num_racks() * bits_per_rack - size()) {
-			rack(num_racks() - 1) &= ~Rack{0} >> pad;
-		}
-	}
-
-	/// Right-rotate by `bits` positions
-	Derived& rotr(const size_t bits)
-	{
-		assert(bits <= size());
-		if (bits == 0 || bits == size()) {
-			return *self();
-		}
-
-		Derived t1(*self());
-		*self() >>= bits;
-		t1 <<= (size() - bits);
-		*self() |= t1;
-
-		truncate();
-		return *self();
-	}
-
-	/// Left-rotate by `bits` positions
-	Derived& rotl(const size_t bits)
-	{
-		assert(bits <= size());
-		if (bits == 0 || bits == size()) {
-			return *self();
-		}
-
-		Derived t1(*self());
-		*self() <<= bits;
-		t1 >>= (size() - bits);
-		*self() |= t1;
-
-		truncate();
-		return *self();
-	}
-
-	/// Return true iff all bits are zero
-	bool none() const
-	{
-		for (size_t i = 0; i < num_racks(); ++i) {
-			if (rack(i)) {
-				return false;
-			}
-		}
-		return true;
-	}
-
-	/// Return 1 + the index of the first set bit, or 0 if there are none
-	size_t find_first() const
-	{
-		for (size_t i = 0; i < num_racks(); ++i) {
-			const int j = chilbert::detail::find_first(rack(i));
-			if (j) {
-				return (i * bits_per_rack) + static_cast<size_t>(j);
-			}
-		}
-		return 0;
-	}
-
-	/// Return the number of set bits
-	size_t count() const
-	{
-		size_t c = 0;
-		for (size_t i = 0; i < num_racks(); ++i) {
-			c += static_cast<size_t>(pop_count(rack(i)));
-		}
-		return c;
-	}
-
-	/// Return a mask that covers the bit with index `i`
-	Mask mask(const size_t i = 0) const
-	{
-		assert(i <= size());
-		return Mask{i};
-	}
-
-	bool operator==(const Derived& vec) const
-	{
-		return (num_racks() == vec.num_racks() &&
-		        (num_racks() == 0 ||
-		         !memcmp(data(), vec.data(), num_racks() * sizeof(Rack))));
-	}
-
-	bool operator!=(const Derived& vec) const { return !(*this == vec); }
-
-	bool operator<(const Derived& vec) const
-	{
-		assert(size() == vec.size());
-
-		for (size_t ri = 0; ri < num_racks(); ++ri) {
-			const size_t i = num_racks() - ri - 1;
-			if (rack(i) < vec.rack(i)) {
-				return true;
-			}
-
-			if (rack(i) > vec.rack(i)) {
-				return false;
-			}
-		}
-		return false;
-	}
-
-	Derived& operator&=(const Derived& vec)
-	{
-		for (size_t i = 0; i < std::min(num_racks(), vec.num_racks()); ++i) {
-			rack(i) &= vec.rack(i);
-		}
-
-		return *self();
-	}
-
-	Derived& operator|=(const Derived& vec)
-	{
-		for (size_t i = 0; i < std::min(num_racks(), vec.num_racks()); ++i) {
-			rack(i) |= vec.rack(i);
-		}
-
-		return *self();
-	}
-
-	Derived& operator^=(const Derived& vec)
-	{
-		for (size_t i = 0; i < std::min(num_racks(), vec.num_racks()); ++i) {
-			rack(i) ^= vec.rack(i);
-		}
-
-		return *self();
-	}
-
-	Derived& operator<<=(const size_t bits)
-	{
-		if (bits == 0) {
-			return *self();
-		}
-
-		if (bits >= size()) {
-			reset();
-			return *self();
-		}
-
-		const Index index{bits};
-
-		if (index.bit == 0) {
-			// Simple rack-aligned shift
-			for (size_t i = num_racks() - 1; i >= index.rack; --i) {
-				rack(i) = rack(i - index.rack);
-			}
-		} else {
-			// Rack + bit offset shift
-			const size_t right_shift = bits_per_rack - index.bit;
-			for (size_t i = num_racks() - index.rack - 1; i > 0; --i) {
-				rack(i + index.rack) =
-				  (rack(i) << index.bit) | (rack(i - 1) >> right_shift);
-			}
-
-			rack(index.rack) = rack(0) << index.bit;
-		}
-
-		// Zero least significant racks
-		for (size_t i = 0; i < index.rack; ++i) {
-			rack(i) = 0;
-		}
-
-		return *self();
-	}
-
-	Derived& operator>>=(const size_t bits)
-	{
-		if (bits == 0) {
-			return *self();
-		}
-
-		if (bits >= size()) {
-			reset();
-			return *self();
-		}
-
-		const Index index{bits};
-
-		if (index.bit == 0) {
-			// Simple rack-aligned shift
-			for (size_t i = 0; i < num_racks() - index.rack; ++i) {
-				rack(i) = rack(i + index.rack);
-			}
-		} else {
-			// Rack + bit offset shift
-			const size_t last       = num_racks() - 1;
-			const size_t left_shift = bits_per_rack - index.bit;
-			for (size_t i = index.rack; i < last; ++i) {
-				rack(i - index.rack) =
-				  (rack(i) >> index.bit) | (rack(i + 1) << left_shift);
-			}
-
-			rack(last - index.rack) = rack(last) >> index.bit;
-		}
-
-		// Zero most significant racks
-		for (size_t i = num_racks() - index.rack; i < num_racks(); ++i) {
-			rack(i) = 0;
-		}
-
-		return *self();
-	}
-
-	auto begin(const size_t i = 0) { return iterator(*self(), i); }
-	auto end() { return iterator(*self(), size()); }
-	auto begin(const size_t i = 0) const { return const_iterator(*self(), i); }
-	auto end() const { return const_iterator(self(), size()); }
-
-	const Rack& rack(const size_t index) const { return self()->rack(index); }
-	Rack&       rack(const size_t index) { return self()->rack(index); }
-	Rack*       data() { return self()->data(); }
-	const Rack* data() const { return self()->data(); }
-	size_t      num_racks() const { return self()->num_racks(); }
-	size_t      size() const { return self()->size(); }
-	size_t      data_size() const { return self()->data_size(); }
+  using Rack = uintptr_t;
+  using Mask = BitVecMask<Rack>;
+
+  using iterator       = BitVecIterator<MultiBitVec<Derived>>;
+  using const_iterator = ConstBitVecIterator<MultiBitVec<Derived>>;
+
+  static constexpr size_t bits_per_rack = sizeof(Rack) * CHAR_BIT;
+
+  /// Return the value of the bit covered by `mask`
+  bool test(const Mask mask) const { return rack(mask.rack) & mask.mask; }
+
+  /// Return the value of the `index`th bit
+  bool test(const size_t index) const { return test(mask(index)); }
+
+  /// Set all bits to one
+  Derived& set()
+  {
+    if (size()) {
+      memset(data(), 0xFF, data_size());
+      self()->truncate();
+    }
+
+    return *self();
+  }
+
+  /// Set the bit covered by `mask` to 1
+  Derived& set(const Mask mask)
+  {
+    rack(mask.rack) |= mask.mask;
+    return *self();
+  }
+
+  /// Set the `index`th bit to 1
+  Derived& set(const size_t index) { return set(mask(index)); }
+
+  /// Set the bit covered by `mask` to `value`
+  Derived& set(const Mask mask, const bool value)
+  {
+    auto& r = rack(mask.rack);
+    r ^= (-Rack{value} ^ r) & mask.mask;
+    return *self();
+  }
+
+  /// Set the `index`th bit to `value`
+  Derived& set(const size_t index, const bool value)
+  {
+    return set(mask(index), value);
+  }
+
+  /// Set all bits to zero
+  Derived& reset()
+  {
+    memset(data(), 0, data_size());
+    return *self();
+  }
+
+  /// Reset the bit covered by `mask` to 0
+  Derived& reset(const Mask mask)
+  {
+    rack(mask.rack) &= ~mask.mask;
+    return *self();
+  }
+
+  /// Reset the `index`th bit to 0
+  Derived& reset(const size_t index) { return reset(mask(index)); }
+
+  /// Flip all bits (one's complement)
+  Derived& flip()
+  {
+    for (size_t i = 0; i < num_racks(); ++i) {
+      rack(i) = ~rack(i);
+    }
+    return *self();
+  }
+
+  /// Flip the value of the bit covered by `mask`
+  Derived& flip(const Mask mask)
+  {
+    rack(mask.rack) ^= mask.mask;
+    return *self();
+  }
+
+  /// Flip the value of the `index`th bit
+  Derived& flip(const size_t index) { return flip(mask(index)); }
+
+  /// Clear any bits in storage outside the valid range if necessary
+  void truncate()
+  {
+    if (const auto pad = num_racks() * bits_per_rack - size()) {
+      rack(num_racks() - 1) &= ~Rack{0} >> pad;
+    }
+  }
+
+  /// Right-rotate by `bits` positions
+  Derived& rotr(const size_t bits)
+  {
+    assert(bits <= size());
+    if (bits == 0 || bits == size()) {
+      return *self();
+    }
+
+    Derived t1(*self());
+    *self() >>= bits;
+    t1 <<= (size() - bits);
+    *self() |= t1;
+
+    truncate();
+    return *self();
+  }
+
+  /// Left-rotate by `bits` positions
+  Derived& rotl(const size_t bits)
+  {
+    assert(bits <= size());
+    if (bits == 0 || bits == size()) {
+      return *self();
+    }
+
+    Derived t1(*self());
+    *self() <<= bits;
+    t1 >>= (size() - bits);
+    *self() |= t1;
+
+    truncate();
+    return *self();
+  }
+
+  /// Return true iff all bits are zero
+  bool none() const
+  {
+    for (size_t i = 0; i < num_racks(); ++i) {
+      if (rack(i)) {
+        return false;
+      }
+    }
+    return true;
+  }
+
+  /// Return 1 + the index of the first set bit, or 0 if there are none
+  size_t find_first() const
+  {
+    for (size_t i = 0; i < num_racks(); ++i) {
+      const int j = chilbert::detail::find_first(rack(i));
+      if (j) {
+        return (i * bits_per_rack) + static_cast<size_t>(j);
+      }
+    }
+    return 0;
+  }
+
+  /// Return the number of set bits
+  size_t count() const
+  {
+    size_t c = 0;
+    for (size_t i = 0; i < num_racks(); ++i) {
+      c += static_cast<size_t>(pop_count(rack(i)));
+    }
+    return c;
+  }
+
+  /// Return a mask that covers the bit with index `i`
+  Mask mask(const size_t i = 0) const
+  {
+    assert(i <= size());
+    return Mask{i};
+  }
+
+  bool operator==(const Derived& vec) const
+  {
+    return (num_racks() == vec.num_racks() &&
+            (num_racks() == 0 ||
+             !memcmp(data(), vec.data(), num_racks() * sizeof(Rack))));
+  }
+
+  bool operator!=(const Derived& vec) const { return !(*this == vec); }
+
+  bool operator<(const Derived& vec) const
+  {
+    assert(size() == vec.size());
+
+    for (size_t ri = 0; ri < num_racks(); ++ri) {
+      const size_t i = num_racks() - ri - 1;
+      if (rack(i) < vec.rack(i)) {
+        return true;
+      }
+
+      if (rack(i) > vec.rack(i)) {
+        return false;
+      }
+    }
+    return false;
+  }
+
+  Derived& operator&=(const Derived& vec)
+  {
+    for (size_t i = 0; i < std::min(num_racks(), vec.num_racks()); ++i) {
+      rack(i) &= vec.rack(i);
+    }
+
+    return *self();
+  }
+
+  Derived& operator|=(const Derived& vec)
+  {
+    for (size_t i = 0; i < std::min(num_racks(), vec.num_racks()); ++i) {
+      rack(i) |= vec.rack(i);
+    }
+
+    return *self();
+  }
+
+  Derived& operator^=(const Derived& vec)
+  {
+    for (size_t i = 0; i < std::min(num_racks(), vec.num_racks()); ++i) {
+      rack(i) ^= vec.rack(i);
+    }
+
+    return *self();
+  }
+
+  Derived& operator<<=(const size_t bits)
+  {
+    if (bits == 0) {
+      return *self();
+    }
+
+    if (bits >= size()) {
+      reset();
+      return *self();
+    }
+
+    const Index index{bits};
+
+    if (index.bit == 0) {
+      // Simple rack-aligned shift
+      for (size_t i = num_racks() - 1; i >= index.rack; --i) {
+        rack(i) = rack(i - index.rack);
+      }
+    } else {
+      // Rack + bit offset shift
+      const size_t right_shift = bits_per_rack - index.bit;
+      for (size_t i = num_racks() - index.rack - 1; i > 0; --i) {
+        rack(i + index.rack) =
+          (rack(i) << index.bit) | (rack(i - 1) >> right_shift);
+      }
+
+      rack(index.rack) = rack(0) << index.bit;
+    }
+
+    // Zero least significant racks
+    for (size_t i = 0; i < index.rack; ++i) {
+      rack(i) = 0;
+    }
+
+    return *self();
+  }
+
+  Derived& operator>>=(const size_t bits)
+  {
+    if (bits == 0) {
+      return *self();
+    }
+
+    if (bits >= size()) {
+      reset();
+      return *self();
+    }
+
+    const Index index{bits};
+
+    if (index.bit == 0) {
+      // Simple rack-aligned shift
+      for (size_t i = 0; i < num_racks() - index.rack; ++i) {
+        rack(i) = rack(i + index.rack);
+      }
+    } else {
+      // Rack + bit offset shift
+      const size_t last       = num_racks() - 1;
+      const size_t left_shift = bits_per_rack - index.bit;
+      for (size_t i = index.rack; i < last; ++i) {
+        rack(i - index.rack) =
+          (rack(i) >> index.bit) | (rack(i + 1) << left_shift);
+      }
+
+      rack(last - index.rack) = rack(last) >> index.bit;
+    }
+
+    // Zero most significant racks
+    for (size_t i = num_racks() - index.rack; i < num_racks(); ++i) {
+      rack(i) = 0;
+    }
+
+    return *self();
+  }
+
+  auto begin(const size_t i = 0) { return iterator(*self(), i); }
+  auto end() { return iterator(*self(), size()); }
+  auto begin(const size_t i = 0) const { return const_iterator(*self(), i); }
+  auto end() const { return const_iterator(self(), size()); }
+
+  const Rack& rack(const size_t index) const { return self()->rack(index); }
+  Rack&       rack(const size_t index) { return self()->rack(index); }
+  Rack*       data() { return self()->data(); }
+  const Rack* data() const { return self()->data(); }
+  size_t      num_racks() const { return self()->num_racks(); }
+  size_t      size() const { return self()->size(); }
+  size_t      data_size() const { return self()->data_size(); }
 
 private:
-	using Index = detail::BitVecIndex<Derived>;
+  using Index = detail::BitVecIndex<Derived>;
 
-	Derived*       self() { return static_cast<Derived*>(this); }
-	const Derived* self() const { return static_cast<const Derived*>(this); }
+  Derived*       self() { return static_cast<Derived*>(this); }
+  const Derived* self() const { return static_cast<const Derived*>(this); }
 };
 
-template <class Derived>
+template<class Derived>
 void
 gray_code(MultiBitVec<Derived>& value)
 {
-	typename MultiBitVec<Derived>::Rack s = 0;
-
-	constexpr size_t left_shift = MultiBitVec<Derived>::bits_per_rack - 1;
-	for (size_t ri = 0; ri < value.num_racks(); ++ri) {
-		const size_t i    = value.num_racks() - ri - 1U;
-		auto&        rack = value.rack(i);
-		const auto   t    = rack & 1U;
-		gray_code(rack);
-		rack ^= (s << left_shift);
-		s = t;
-	}
+  typename MultiBitVec<Derived>::Rack s = 0;
+
+  constexpr size_t left_shift = MultiBitVec<Derived>::bits_per_rack - 1;
+  for (size_t ri = 0; ri < value.num_racks(); ++ri) {
+    const size_t i    = value.num_racks() - ri - 1U;
+    auto&        rack = value.rack(i);
+    const auto   t    = rack & 1U;
+    gray_code(rack);
+    rack ^= (s << left_shift);
+    s = t;
+  }
 }
 
-template <class Derived>
+template<class Derived>
 void
 gray_code_inv(MultiBitVec<Derived>& value)
 {
-	using Rack = typename MultiBitVec<Derived>::Rack;
-
-	constexpr std::array<Rack, 2> masks{Rack{0}, ~Rack{0}};
-	bool                          s = false;
-
-	for (size_t ri = 0; ri < value.num_racks(); ++ri) {
-		const size_t i    = value.num_racks() - ri - 1;
-		auto&        rack = value.rack(i);
-		gray_code_inv(rack);
-		rack ^= masks[s];
-		s = rack & 1U;
-	}
+  using Rack = typename MultiBitVec<Derived>::Rack;
+
+  constexpr std::array<Rack, 2> masks{Rack{0}, ~Rack{0}};
+  bool                          s = false;
+
+  for (size_t ri = 0; ri < value.num_racks(); ++ri) {
+    const size_t i    = value.num_racks() - ri - 1;
+    auto&        rack = value.rack(i);
+    gray_code_inv(rack);
+    rack ^= masks[s];
+    s = rack & 1U;
+  }
 }
 
 } // namespace detail
diff --git a/include/chilbert/detail/gray_code_rank.hpp b/include/chilbert/detail/gray_code_rank.hpp
index b6811ff..a326617 100644
--- a/include/chilbert/detail/gray_code_rank.hpp
+++ b/include/chilbert/detail/gray_code_rank.hpp
@@ -30,7 +30,7 @@ namespace detail {
 // a Compact Hilbert Index.  It compresses a previously
 // calculated index when provided the rotation
 // at each level of precision.
-template <class H, class HC>
+template<class H, class HC>
 inline void
 compact_index(const size_t* const ms,
               const size_t* const ds,
@@ -39,60 +39,60 @@ compact_index(const size_t* const ms,
               H&                  h,
               HC&                 hc)
 {
-	assert(h.size() >= n * m);
-	assert(hc.size() >= std::accumulate(ms, ms + n, size_t(0)));
-
-	reset_bits(hc);
-
-	auto hm  = h.mask(0);
-	auto hcm = hc.mask(0);
-
-	// Run through the levels of precision
-	for (size_t i = 0; i < m; ++i) {
-		// Run through the dimensions
-		size_t j = ds[i];
-		do {
-			// This dimension contributes a bit?
-			if (ms[j] > i) {
-				if (h.test(hm)) {
-					hc.set(hcm);
-				}
-				++hcm;
-			}
-
-			if (++j == n) {
-				j = 0;
-			}
-			++hm;
-		} while (j != ds[i]);
-	}
+  assert(h.size() >= n * m);
+  assert(hc.size() >= std::accumulate(ms, ms + n, size_t(0)));
+
+  reset_bits(hc);
+
+  auto hm  = h.mask(0);
+  auto hcm = hc.mask(0);
+
+  // Run through the levels of precision
+  for (size_t i = 0; i < m; ++i) {
+    // Run through the dimensions
+    size_t j = ds[i];
+    do {
+      // This dimension contributes a bit?
+      if (ms[j] > i) {
+        if (h.test(hm)) {
+          hc.set(hcm);
+        }
+        ++hcm;
+      }
+
+      if (++j == n) {
+        j = 0;
+      }
+      ++hm;
+    } while (j != ds[i]);
+  }
 }
 
-template <class I>
+template<class I>
 inline void
 gray_code_rank(const I& mask, const I& gi, const size_t n, I& r)
 {
-	assert(mask.size() == n);
-	assert(gi.size() == n);
-	assert(r.size() == n);
-
-	r.reset();
-
-	auto mi  = mask.begin();
-	auto gii = gi.begin();
-	auto ri  = r.begin();
-
-	for (size_t i = 0; i < n; ++i, ++mi, ++gii) {
-		if (*mi) {
-			if (*gii) {
-				ri.set();
-			}
-			++ri;
-		}
-	}
+  assert(mask.size() == n);
+  assert(gi.size() == n);
+  assert(r.size() == n);
+
+  r.reset();
+
+  auto mi  = mask.begin();
+  auto gii = gi.begin();
+  auto ri  = r.begin();
+
+  for (size_t i = 0; i < n; ++i, ++mi, ++gii) {
+    if (*mi) {
+      if (*gii) {
+        ri.set();
+      }
+      ++ri;
+    }
+  }
 }
 
-template <class I>
+template<class I>
 inline void
 gray_code_rank_inv(const I&     mask,
                    const I&     ptrn,
@@ -102,51 +102,51 @@ gray_code_rank_inv(const I&     mask,
                    I&           g,
                    I&           gi)
 {
-	assert(mask.size() == n);
-	assert(ptrn.size() == n);
-	assert(r.size() == n);
-	assert(g.size() == n);
-	assert(gi.size() == n);
-
-	g.reset();
-	gi.reset();
-
-	auto m  = mask.mask(n - 1);
-	auto ri = r.begin(b - 1);
-
-	typename I::Rack gi0 = 0;
-	typename I::Rack gi1 = 0;
-	typename I::Rack g0  = 0;
-
-	for (size_t i = 0; i < n; ++i) {
-		if (mask.test(m)) { // Unconstrained bit
-			gi1 = gi0;
-			gi0 = *ri;
-			g0  = gi0 ^ gi1;
-			if (gi0) {
-				gi.set(m);
-			}
-			if (g0) {
-				g.set(m);
-			}
-			--ri;
-		} else { // Constrained bit
-			g0  = (ptrn.test(m) > 0);
-			gi1 = gi0;
-			gi0 = g0 ^ gi1;
-			if (gi0) {
-				gi.set(m);
-			}
-			if (g0) {
-				g.set(m);
-			}
-		}
-
-		--m;
-	}
+  assert(mask.size() == n);
+  assert(ptrn.size() == n);
+  assert(r.size() == n);
+  assert(g.size() == n);
+  assert(gi.size() == n);
+
+  g.reset();
+  gi.reset();
+
+  auto m  = mask.mask(n - 1);
+  auto ri = r.begin(b - 1);
+
+  typename I::Rack gi0 = 0;
+  typename I::Rack gi1 = 0;
+  typename I::Rack g0  = 0;
+
+  for (size_t i = 0; i < n; ++i) {
+    if (mask.test(m)) { // Unconstrained bit
+      gi1 = gi0;
+      gi0 = *ri;
+      g0  = gi0 ^ gi1;
+      if (gi0) {
+        gi.set(m);
+      }
+      if (g0) {
+        g.set(m);
+      }
+      --ri;
+    } else { // Constrained bit
+      g0  = (ptrn.test(m) > 0);
+      gi1 = gi0;
+      gi0 = g0 ^ gi1;
+      if (gi0) {
+        gi.set(m);
+      }
+      if (g0) {
+        g.set(m);
+      }
+    }
+
+    --m;
+  }
 }
 
-template <class I>
+template<class I>
 inline void
 extract_mask(const size_t* const ms,
              const size_t        n,
@@ -155,25 +155,25 @@ extract_mask(const size_t* const ms,
              I&                  mask,
              size_t&             b)
 {
-	assert(d < n);
-	assert(mask.size() == n);
-
-	mask.reset();
-	b = 0;
-
-	auto   mi = mask.begin();
-	size_t j  = d; // #D j = (d==n-1) ? 0 : d+1;
-	do {
-		if (ms[j] > i) {
-			mi.set();
-			++b;
-		}
-
-		++mi;
-		if (++j == n) {
-			j = 0;
-		}
-	} while (j != d);
+  assert(d < n);
+  assert(mask.size() == n);
+
+  mask.reset();
+  b = 0;
+
+  auto   mi = mask.begin();
+  size_t j  = d; // #D j = (d==n-1) ? 0 : d+1;
+  do {
+    if (ms[j] > i) {
+      mi.set();
+      ++b;
+    }
+
+    ++mi;
+    if (++j == n) {
+      j = 0;
+    }
+  } while (j != d);
 }
 
 } // namespace detail
diff --git a/include/chilbert/detail/operations.hpp b/include/chilbert/detail/operations.hpp
index a196ba4..ed9dfac 100644
--- a/include/chilbert/detail/operations.hpp
+++ b/include/chilbert/detail/operations.hpp
@@ -32,135 +32,139 @@ namespace chilbert {
 namespace detail {
 
 /// Reset all bits in `field`
-template <typename T>
+template<typename T>
 std::enable_if_t<std::is_integral<T>::value>
 reset_bits(T& field)
 {
-	field = static_cast<T>(0);
+  field = static_cast<T>(0);
 }
 
 /// Reset all bits in `field`
-template <typename T>
+template<typename T>
 std::enable_if_t<is_bitvec_v<T>>
 reset_bits(T& field)
 {
-	field.reset();
+  field.reset();
 }
 
 /// Return the `index`th bit in `field`
-template <typename T>
+template<typename T>
 std::enable_if_t<std::is_integral<T>::value, bool>
 test_bit(const T& field, const size_t index)
 {
-	assert(size_t(index) < sizeof(field) * CHAR_BIT);
-	return field & (T{1} << index);
+  assert(size_t(index) < sizeof(field) * CHAR_BIT);
+  return field & (T{1} << index);
 }
 
 /// Return the `index`th bit in `field`
-template <typename T>
+template<typename T>
 std::enable_if_t<is_bitvec_v<T>, bool>
 test_bit(const T& field, const size_t index)
 {
-	return field.test(index);
+  return field.test(index);
 }
 
 /// Set the `index`th bit in `field`
-template <typename T>
+template<typename T>
 std::enable_if_t<std::is_integral<T>::value>
 set_bit(T& field, const size_t index)
 {
-	assert(size_t(index) < sizeof(field) * CHAR_BIT);
-	field |= (T{1} << index);
-	assert(test_bit(field, index));
+  assert(size_t(index) < sizeof(field) * CHAR_BIT);
+  field |= (T{1} << index);
+  assert(test_bit(field, index));
 }
 
 /// Set the `index`th bit in `field` to `value`
-template <typename T>
+template<typename T>
 std::enable_if_t<std::is_integral<T>::value>
 set_bit(T& field, const size_t index, const bool value)
 {
-	assert(size_t(index) < sizeof(field) * CHAR_BIT);
-	field ^= (-T{value} ^ field) & (T{1U} << index);
-	assert(test_bit(field, index) == value);
+  assert(size_t(index) < sizeof(field) * CHAR_BIT);
+  field ^= (-T{value} ^ field) & (T{1U} << index);
+  assert(test_bit(field, index) == value);
 }
 
 /// Set the `index`th bit in `field`
-template <typename T>
+template<typename T>
 std::enable_if_t<is_bitvec_v<T>>
 set_bit(T& field, const size_t index)
 {
-	field.set(index);
+  field.set(index);
 }
 
 /// Set the `index`th bit in `field` to `value`
-template <typename T>
+template<typename T>
 std::enable_if_t<is_bitvec_v<T>>
 set_bit(T& field, const size_t index, const bool value)
 {
-	field.set(index, value);
+  field.set(index, value);
 }
 
 /// Return 1 + the index of the least significant 1-bit of `field`, or zero
-template <typename T>
-inline int pop_count(const T& field);
+template<typename T>
+inline int
+pop_count(const T& field);
 
-template <>
+template<>
 inline int
 pop_count<unsigned long>(const unsigned long& field)
 {
-	return __builtin_popcountl(field);
+  return __builtin_popcountl(field);
 }
 
-template <>
+template<>
 inline int
 pop_count<unsigned long long>(const unsigned long long& field)
 {
-	return __builtin_popcountll(field);
+  return __builtin_popcountll(field);
 }
 
 /// Return 1 + the index of the least significant 1-bit of `field`, or zero
-template <typename T>
-inline int find_first(const T field);
+template<typename T>
+inline int
+find_first(const T field);
 
-template <>
+template<>
 inline int
 find_first<unsigned long>(const unsigned long field)
 {
-	return __builtin_ffsl(static_cast<long>(field));
+  return __builtin_ffsl(static_cast<long>(field));
 }
 
-template <>
+template<>
 inline int
 find_first<unsigned long long>(const unsigned long long field)
 {
-	return __builtin_ffsll(static_cast<long long>(field));
+  return __builtin_ffsll(static_cast<long long>(field));
 }
 
 /// Calculates the Gray Code of `value` in place
-template <typename T>
-std::enable_if_t<is_bitvec_v<T>> gray_code(T& value);
+template<typename T>
+std::enable_if_t<is_bitvec_v<T>>
+gray_code(T& value);
 
 /// Implementation of grayCode for any integral type
-template <typename T>
+template<typename T>
 std::enable_if_t<std::is_integral<T>::value>
 gray_code(T& value)
 {
-	value ^= (value >> 1U);
+  value ^= (value >> 1U);
 }
 
 /// Calculates the inverse Gray Code of `value` in place
-template <typename T>
-std::enable_if_t<is_bitvec_v<T>> gray_code_inv(T& value);
+template<typename T>
+std::enable_if_t<is_bitvec_v<T>>
+gray_code_inv(T& value);
 
 /// Implementation of gray_code_inv for any integral type
-template <typename T>
+template<typename T>
 std::enable_if_t<std::is_integral<T>::value>
 gray_code_inv(T& value)
 {
-	constexpr T shift_end = sizeof(T) * CHAR_BIT;
-	for (T shift = 1U; shift < shift_end; shift <<= 1U) {
-		value ^= (value >> shift);
-	}
+  constexpr T shift_end = sizeof(T) * CHAR_BIT;
+  for (T shift = 1U; shift < shift_end; shift <<= 1U) {
+    value ^= (value >> shift);
+  }
 }
 
 } // namespace detail
diff --git a/include/chilbert/detail/traits.hpp b/include/chilbert/detail/traits.hpp
index 3e3f4d2..e6a2023 100644
--- a/include/chilbert/detail/traits.hpp
+++ b/include/chilbert/detail/traits.hpp
@@ -23,14 +23,13 @@ namespace chilbert {
 namespace detail {
 
 /// Member `value` is true iff T is a chilbert bitset
-template <class T>
-struct is_bitvec
-{
-	static constexpr bool value = false;
+template<class T>
+struct is_bitvec {
+  static constexpr bool value = false;
 };
 
 /// True iff T is a chilbert bitset
-template <class T>
+template<class T>
 static constexpr bool is_bitvec_v = is_bitvec<T>::value;
 
 } // namespace detail
diff --git a/include/chilbert/operators.hpp b/include/chilbert/operators.hpp
index efad7f4..4a92e0d 100644
--- a/include/chilbert/operators.hpp
+++ b/include/chilbert/operators.hpp
@@ -28,67 +28,68 @@ namespace chilbert {
 
 using detail::is_bitvec_v;
 
-template <class T>
-std::enable_if_t<is_bitvec_v<T>, T> operator&(const T& lhs, const T& rhs)
+template<class T>
+std::enable_if_t<is_bitvec_v<T>, T>
+operator&(const T& lhs, const T& rhs)
 {
-	T r{lhs};
-	r &= rhs;
-	return r;
+  T r{lhs};
+  r &= rhs;
+  return r;
 }
 
-template <class T>
+template<class T>
 std::enable_if_t<is_bitvec_v<T>, T>
 operator|(const T& lhs, const T& rhs)
 {
-	T r{lhs};
-	r |= rhs;
-	return r;
+  T r{lhs};
+  r |= rhs;
+  return r;
 }
 
-template <class T>
+template<class T>
 std::enable_if_t<is_bitvec_v<T>, T>
 operator^(const T& lhs, const T& rhs)
 {
-	T r{lhs};
-	r ^= rhs;
-	return r;
+  T r{lhs};
+  r ^= rhs;
+  return r;
 }
 
-template <class T>
+template<class T>
 std::enable_if_t<is_bitvec_v<T>, T>
 operator~(const T& vec)
 {
-	T r{vec};
-	r.flip();
-	return r;
+  T r{vec};
+  r.flip();
+  return r;
 }
 
-template <class T>
+template<class T>
 std::enable_if_t<is_bitvec_v<T>, T>
 operator<<(const T& vec, const size_t bits)
 {
-	T r{vec};
-	r <<= bits;
-	return r;
+  T r{vec};
+  r <<= bits;
+  return r;
 }
 
-template <class T>
+template<class T>
 std::enable_if_t<is_bitvec_v<T>, T>
 operator>>(const T& vec, const size_t bits)
 {
-	T r{vec};
-	r >>= bits;
-	return r;
+  T r{vec};
+  r >>= bits;
+  return r;
 }
 
-template <class T>
+template<class T>
 inline std::enable_if_t<is_bitvec_v<T>, std::ostream>&
 operator<<(std::ostream& os, const T& vec)
 {
-	for (size_t i = 0; i < vec.size(); ++i) {
-		os << vec.test(vec.size() - i - 1);
-	}
-	return os;
+  for (size_t i = 0; i < vec.size(); ++i) {
+    os << vec.test(vec.size() - i - 1);
+  }
+  return os;
 }
 
 } // namespace chilbert
diff --git a/src/chilbert_obj.cpp b/src/chilbert_obj.cpp
index 0f31fbd..1fdd1ee 100644
--- a/src/chilbert_obj.cpp
+++ b/src/chilbert_obj.cpp
@@ -25,30 +25,30 @@
 int
 main(int argc, char** argv)
 {
-	if (argc != 2) {
-		fprintf(stderr, "Usage: %s NUM_POINTS\n", argv[0]);
-		return 1;
-	}
-
-	const unsigned long num_points = std::strtoul(argv[1], nullptr, 10);
-	if (num_points == 0 || num_points == ULONG_MAX) {
-		fprintf(stderr, "Usage: %s NUM_POINTS\n", argv[0]);
-		return 1;
-	}
-
-	// Vertices
-	for (uint64_t i = 0; i < num_points; ++i) {
-		std::array<uint32_t, 3> point;
-		chilbert::index_to_coords(point, 16, 3, i);
-		printf("v %u %u %u\n", point[0], point[1], point[2]);
-	}
-
-	// One polyline through all vertices
-	printf("\nl");
-	for (unsigned i = 0; i < num_points - 1; ++i) {
-		printf(" %u", i + 1);
-	}
-	printf("\n");
-
-	return 0;
+  if (argc != 2) {
+    fprintf(stderr, "Usage: %s NUM_POINTS\n", argv[0]);
+    return 1;
+  }
+
+  const unsigned long num_points = std::strtoul(argv[1], nullptr, 10);
+  if (num_points == 0 || num_points == ULONG_MAX) {
+    fprintf(stderr, "Usage: %s NUM_POINTS\n", argv[0]);
+    return 1;
+  }
+
+  // Vertices
+  for (uint64_t i = 0; i < num_points; ++i) {
+    std::array<uint32_t, 3> point;
+    chilbert::index_to_coords(point, 16, 3, i);
+    printf("v %u %u %u\n", point[0], point[1], point[2]);
+  }
+
+  // One polyline through all vertices
+  printf("\nl");
+  for (unsigned i = 0; i < num_points - 1; ++i) {
+    printf(" %u", i + 1);
+  }
+  printf("\n");
+
+  return 0;
 }
diff --git a/src/chilbert_svg.cpp b/src/chilbert_svg.cpp
index dd4a3d4..7033bec 100644
--- a/src/chilbert_svg.cpp
+++ b/src/chilbert_svg.cpp
@@ -26,38 +26,38 @@
 int
 main(int argc, char** argv)
 {
-	if (argc != 2) {
-		fprintf(stderr, "Usage: %s NUM_POINTS\n", argv[0]);
-		return 1;
-	}
-
-	const unsigned long num_points = std::strtoul(argv[1], nullptr, 10);
-	if (num_points == 0 || num_points == ULONG_MAX) {
-		fprintf(stderr, "Usage: %s NUM_POINTS\n", argv[0]);
-		return 1;
-	}
-
-	const uint32_t w =
-	  uint32_t(sqrt(1 << uint32_t(ceil(log2(double(num_points)))))) - 1;
-
-	// Header
-	printf("<svg xmlns='http://www.w3.org/2000/svg'"
-	       " version='1.2' baseProfile='tiny' viewBox='0 0 %u %u'>\n",
-	       w,
-	       w);
-	printf("<desc>Hilbert Curve</desc>\n");
-	printf("<polyline vector-effect='non-scaling-stroke' fill='none' "
-	       "stroke='black' stroke-width='1' points='");
-
-	// One polyline through all vertices
-	for (uint64_t i = 0; i <= num_points; ++i) {
-		std::array<uint32_t, 2> point;
-		chilbert::index_to_coords(point, 32, 2, i);
-		printf("%u,%u ", point[0], point[1]);
-	}
-
-	// Close off document
-	printf("' />\n</svg>\n");
-
-	return 0;
+  if (argc != 2) {
+    fprintf(stderr, "Usage: %s NUM_POINTS\n", argv[0]);
+    return 1;
+  }
+
+  const unsigned long num_points = std::strtoul(argv[1], nullptr, 10);
+  if (num_points == 0 || num_points == ULONG_MAX) {
+    fprintf(stderr, "Usage: %s NUM_POINTS\n", argv[0]);
+    return 1;
+  }
+
+  const uint32_t w =
+    uint32_t(sqrt(1 << uint32_t(ceil(log2(double(num_points)))))) - 1;
+
+  // Header
+  printf("<svg xmlns='http://www.w3.org/2000/svg'"
+         " version='1.2' baseProfile='tiny' viewBox='0 0 %u %u'>\n",
+         w,
+         w);
+  printf("<desc>Hilbert Curve</desc>\n");
+  printf("<polyline vector-effect='non-scaling-stroke' fill='none' "
+         "stroke='black' stroke-width='1' points='");
+
+  // One polyline through all vertices
+  for (uint64_t i = 0; i <= num_points; ++i) {
+    std::array<uint32_t, 2> point;
+    chilbert::index_to_coords(point, 32, 2, i);
+    printf("%u,%u ", point[0], point[1]);
+  }
+
+  // Close off document
+  printf("' />\n</svg>\n");
+
+  return 0;
 }
diff --git a/test/test_bitvec.cpp b/test/test_bitvec.cpp
index c123212..892744a 100644
--- a/test/test_bitvec.cpp
+++ b/test/test_bitvec.cpp
@@ -30,357 +30,357 @@
 #include <cstddef>
 #include <cstdlib>
 
-template <class T, size_t N>
+template<class T, size_t N>
 void
 test_and(Context& ctx)
 {
-	const T a = make_random_bitvec<T, N>(ctx);
-	const T b = make_random_bitvec<T, N>(ctx);
-	T       r = a;
-	assert((a & b) == (r &= b));
-
-	for (size_t i = 0; i < N; ++i) {
-		assert(r.test(i) == (a.test(i) && b.test(i)));
-	}
+  const T a = make_random_bitvec<T, N>(ctx);
+  const T b = make_random_bitvec<T, N>(ctx);
+  T       r = a;
+  assert((a & b) == (r &= b));
+
+  for (size_t i = 0; i < N; ++i) {
+    assert(r.test(i) == (a.test(i) && b.test(i)));
+  }
 }
 
-template <class T, size_t N>
+template<class T, size_t N>
 void
 test_or(Context& ctx)
 {
-	const T a = make_random_bitvec<T, N>(ctx);
-	const T b = make_random_bitvec<T, N>(ctx);
-	T       r = a;
-	assert((a | b) == (r |= b));
-
-	for (size_t i = 0; i < N; ++i) {
-		assert(r.test(i) == (a.test(i) || b.test(i)));
-	}
+  const T a = make_random_bitvec<T, N>(ctx);
+  const T b = make_random_bitvec<T, N>(ctx);
+  T       r = a;
+  assert((a | b) == (r |= b));
+
+  for (size_t i = 0; i < N; ++i) {
+    assert(r.test(i) == (a.test(i) || b.test(i)));
+  }
 }
 
-template <class T, size_t N>
+template<class T, size_t N>
 void
 test_xor(Context& ctx)
 {
-	const T a = make_random_bitvec<T, N>(ctx);
-	const T b = make_random_bitvec<T, N>(ctx);
-	T       r = a;
-	assert((a ^ b) == (r ^= b));
-
-	for (size_t i = 0; i < N; ++i) {
-		assert(r.test(i) == (a.test(i) != b.test(i)));
-	}
+  const T a = make_random_bitvec<T, N>(ctx);
+  const T b = make_random_bitvec<T, N>(ctx);
+  T       r = a;
+  assert((a ^ b) == (r ^= b));
+
+  for (size_t i = 0; i < N; ++i) {
+    assert(r.test(i) == (a.test(i) != b.test(i)));
+  }
 }
 
-template <class T, size_t N>
+template<class T, size_t N>
 void
 test_not(Context& ctx)
 {
-	const T v = make_random_bitvec<T, N>(ctx);
-	const T r = ~v;
+  const T v = make_random_bitvec<T, N>(ctx);
+  const T r = ~v;
 
-	for (size_t i = 0; i < N; ++i) {
-		assert(r.test(i) == !v.test(i));
-	}
+  for (size_t i = 0; i < N; ++i) {
+    assert(r.test(i) == !v.test(i));
+  }
 }
 
-template <class T, size_t N>
+template<class T, size_t N>
 void
 test_flip_one(Context&)
 {
-	T v = make_zero_bitvec<T, N>();
-	for (size_t i = 0; i < N; ++i) {
-		assert(v.none());
-		v.flip(i);
-		for (size_t j = 0; j < N; ++j) {
-			assert(v.test(j) == (j == i));
-		}
-
-		v.flip(i);
-		assert(v.none());
-	}
+  T v = make_zero_bitvec<T, N>();
+  for (size_t i = 0; i < N; ++i) {
+    assert(v.none());
+    v.flip(i);
+    for (size_t j = 0; j < N; ++j) {
+      assert(v.test(j) == (j == i));
+    }
+
+    v.flip(i);
+    assert(v.none());
+  }
 }
 
-template <class T, size_t N>
+template<class T, size_t N>
 void
 test_flip_all(Context& ctx)
 {
-	const T a = make_random_bitvec<T, N>(ctx);
-	T       r = a;
-	r.flip();
-	for (size_t i = 0; i < N; ++i) {
-		assert(r.test(i) == !a.test(i));
-	}
+  const T a = make_random_bitvec<T, N>(ctx);
+  T       r = a;
+  r.flip();
+  for (size_t i = 0; i < N; ++i) {
+    assert(r.test(i) == !a.test(i));
+  }
 }
 
-template <class T, size_t N>
+template<class T, size_t N>
 void
 test_none(Context&)
 {
-	T v = make_zero_bitvec<T, N>();
-	assert(v.none());
-	v.set();
-	assert(v.none() == (N == 0));
-	if (N > 1) {
-		v.reset(N / 2);
-		assert(!v.none());
-		v.reset();
-		v.set(N / 2);
-		assert(!v.none());
-	}
+  T v = make_zero_bitvec<T, N>();
+  assert(v.none());
+  v.set();
+  assert(v.none() == (N == 0));
+  if (N > 1) {
+    v.reset(N / 2);
+    assert(!v.none());
+    v.reset();
+    v.set(N / 2);
+    assert(!v.none());
+  }
 }
 
-template <class T, size_t N>
+template<class T, size_t N>
 void
 test_set_reset_one(Context&)
 {
-	T v = make_zero_bitvec<T, N>();
-	for (size_t i = 0; i < N; ++i) {
-		assert(v.none());
-		v.set(i);
-		for (size_t j = 0; j < N; ++j) {
-			assert(v.test(j) == (j == i));
-		}
-
-		v.reset(i);
-		assert(v.none());
-	}
+  T v = make_zero_bitvec<T, N>();
+  for (size_t i = 0; i < N; ++i) {
+    assert(v.none());
+    v.set(i);
+    for (size_t j = 0; j < N; ++j) {
+      assert(v.test(j) == (j == i));
+    }
+
+    v.reset(i);
+    assert(v.none());
+  }
 }
 
-template <class T, size_t N>
+template<class T, size_t N>
 void
 test_set_all(Context&)
 {
-	T v = make_zero_bitvec<T, N>();
-	v.set();
-	for (size_t i = 0; i < N; ++i) {
-		assert(v.test(i));
-	}
+  T v = make_zero_bitvec<T, N>();
+  v.set();
+  for (size_t i = 0; i < N; ++i) {
+    assert(v.test(i));
+  }
 }
 
-template <class T, size_t N>
+template<class T, size_t N>
 void
 test_reset_all(Context&)
 {
-	T v = make_zero_bitvec<T, N>();
-	v.set();
-	v.reset();
-	for (size_t i = 0; i < N; ++i) {
-		assert(!v.test(i));
-	}
+  T v = make_zero_bitvec<T, N>();
+  v.set();
+  v.reset();
+  for (size_t i = 0; i < N; ++i) {
+    assert(!v.test(i));
+  }
 }
 
-template <class T, size_t N>
+template<class T, size_t N>
 void
 test_left_shift(Context& ctx)
 {
-	for (size_t s = 0; s < N; ++s) {
-		const T v = make_random_bitvec<T, N>(ctx);
-		T       r = v;
-		assert((v << s) == (r <<= s));
-
-		for (size_t i = 0; i < s; ++i) {
-			assert(!r.test(i));
-		}
-
-		for (size_t i = s; i < N; ++i) {
-			assert(r.test(i) == v.test(i - s));
-		}
-	}
+  for (size_t s = 0; s < N; ++s) {
+    const T v = make_random_bitvec<T, N>(ctx);
+    T       r = v;
+    assert((v << s) == (r <<= s));
+
+    for (size_t i = 0; i < s; ++i) {
+      assert(!r.test(i));
+    }
+
+    for (size_t i = s; i < N; ++i) {
+      assert(r.test(i) == v.test(i - s));
+    }
+  }
 }
 
-template <class T, size_t N>
+template<class T, size_t N>
 void
 test_right_shift(Context& ctx)
 {
-	for (size_t s = 0; s < N; ++s) {
-		const T v = make_random_bitvec<T, N>(ctx);
-		T       r = v;
-		assert((v >> s) == (r >>= s));
-
-		for (size_t i = N - 1; i > N - s - 1; --i) {
-			assert(!r.test(i));
-		}
-
-		for (size_t i = 0; i <= N - s - 1; ++i) {
-			assert(r.test(i) == v.test(i + s));
-		}
-	}
+  for (size_t s = 0; s < N; ++s) {
+    const T v = make_random_bitvec<T, N>(ctx);
+    T       r = v;
+    assert((v >> s) == (r >>= s));
+
+    for (size_t i = N - 1; i > N - s - 1; --i) {
+      assert(!r.test(i));
+    }
+
+    for (size_t i = 0; i <= N - s - 1; ++i) {
+      assert(r.test(i) == v.test(i + s));
+    }
+  }
 }
 
-template <class T, size_t N>
+template<class T, size_t N>
 void
 test_left_rotate(Context& ctx)
 {
-	const T v = make_random_bitvec<T, N>(ctx);
-	for (size_t bits = 0; bits <= N; ++bits) {
-		T r = v;
-		r.rotl(bits);
-
-		if (N > 0) {
-			for (size_t i = 0; i < N; ++i) {
-				assert(r.test((i + bits) % N) == v.test(i));
-			}
-		}
-	}
+  const T v = make_random_bitvec<T, N>(ctx);
+  for (size_t bits = 0; bits <= N; ++bits) {
+    T r = v;
+    r.rotl(bits);
+
+    if (N > 0) {
+      for (size_t i = 0; i < N; ++i) {
+        assert(r.test((i + bits) % N) == v.test(i));
+      }
+    }
+  }
 }
 
-template <class T, size_t N>
+template<class T, size_t N>
 void
 test_right_rotate(Context& ctx)
 {
-	const T v = make_random_bitvec<T, N>(ctx);
-	for (size_t bits = 0; bits <= N; ++bits) {
-		T r = v;
-		r.rotr(bits);
-
-		if (N > 0) {
-			for (size_t i = 0; i < N; ++i) {
-				assert(r.test(i) == v.test((i + bits) % N));
-			}
-		}
-	}
+  const T v = make_random_bitvec<T, N>(ctx);
+  for (size_t bits = 0; bits <= N; ++bits) {
+    T r = v;
+    r.rotr(bits);
+
+    if (N > 0) {
+      for (size_t i = 0; i < N; ++i) {
+        assert(r.test(i) == v.test((i + bits) % N));
+      }
+    }
+  }
 }
 
-template <class T, size_t N>
+template<class T, size_t N>
 void
 test_find_first(Context&)
 {
-	T v = make_zero_bitvec<T, N>();
-	for (size_t i = 0; i < N; ++i) {
-		v.reset();
-		v.set(i);
-		for (size_t j = i + 1; j < N; ++j) {
-			v.set(j, rand() % 2);
-		}
-		assert(size_t(v.find_first()) == i + 1);
-	}
+  T v = make_zero_bitvec<T, N>();
+  for (size_t i = 0; i < N; ++i) {
+    v.reset();
+    v.set(i);
+    for (size_t j = i + 1; j < N; ++j) {
+      v.set(j, rand() % 2);
+    }
+    assert(size_t(v.find_first()) == i + 1);
+  }
 }
 
-template <class T, size_t N>
+template<class T, size_t N>
 void
 test_gray_code(Context& ctx)
 {
-	const T v = make_random_bitvec<T, N>(ctx);
-	T       r = v;
-	chilbert::detail::gray_code(r);
+  const T v = make_random_bitvec<T, N>(ctx);
+  T       r = v;
+  chilbert::detail::gray_code(r);
 
-	if (N > 0) {
-		assert(N == 1 || r == (v ^ (v >> 1)));
+  if (N > 0) {
+    assert(N == 1 || r == (v ^ (v >> 1)));
 
-		T s = r;
-		chilbert::detail::gray_code_inv(s);
-		assert(s == v);
-	}
+    T s = r;
+    chilbert::detail::gray_code_inv(s);
+    assert(s == v);
+  }
 }
 
-template <class T, size_t N>
+template<class T, size_t N>
 void
 test_comparison(Context&)
 {
-	T a = make_zero_bitvec<T, N>();
-	T b = make_zero_bitvec<T, N>();
+  T a = make_zero_bitvec<T, N>();
+  T b = make_zero_bitvec<T, N>();
 
-	for (size_t bit = 1; bit < N; ++bit) {
-		chilbert::detail::set_bit(a, bit, 1);
+  for (size_t bit = 1; bit < N; ++bit) {
+    chilbert::detail::set_bit(a, bit, 1);
 
-		for (size_t i = 0; i < bit; ++i) {
-			chilbert::detail::set_bit(a, i, rand() % 2 == 0);
-			chilbert::detail::set_bit(b, i, rand() % 2 == 0);
-		}
+    for (size_t i = 0; i < bit; ++i) {
+      chilbert::detail::set_bit(a, i, rand() % 2 == 0);
+      chilbert::detail::set_bit(b, i, rand() % 2 == 0);
+    }
 
-		assert(b < a);
-	}
+    assert(b < a);
+  }
 }
 
-template <class T, size_t N>
+template<class T, size_t N>
 void
 test_iteration(Context&)
 {
-	T      v     = make_zero_bitvec<T, N>();
-	size_t count = 0;
-	for (const auto bit : v) {
-		assert(!bit);
-		++count;
-	}
-	assert(count == N);
-
-	v.flip();
-	count = 0;
-	for (const auto bit : v) {
-		assert(bit);
-		++count;
-	}
-	assert(count == N);
+  T      v     = make_zero_bitvec<T, N>();
+  size_t count = 0;
+  for (const auto bit : v) {
+    assert(!bit);
+    ++count;
+  }
+  assert(count == N);
+
+  v.flip();
+  count = 0;
+  for (const auto bit : v) {
+    assert(bit);
+    ++count;
+  }
+  assert(count == N);
 }
 
-template <class T, size_t N>
+template<class T, size_t N>
 void
 test(Context& ctx)
 {
-	test_and<T, N>(ctx);
-	test_or<T, N>(ctx);
-	test_xor<T, N>(ctx);
-	test_not<T, N>(ctx);
-	test_flip_one<T, N>(ctx);
-	test_flip_all<T, N>(ctx);
-	test_none<T, N>(ctx);
-	test_set_reset_one<T, N>(ctx);
-	test_set_all<T, N>(ctx);
-	test_reset_all<T, N>(ctx);
-	test_left_shift<T, N>(ctx);
-	test_right_shift<T, N>(ctx);
-	test_left_rotate<T, N>(ctx);
-	test_right_rotate<T, N>(ctx);
-	test_find_first<T, N>(ctx);
-	test_gray_code<T, N>(ctx);
-	test_comparison<T, N>(ctx);
-	test_iteration<T, N>(ctx);
+  test_and<T, N>(ctx);
+  test_or<T, N>(ctx);
+  test_xor<T, N>(ctx);
+  test_not<T, N>(ctx);
+  test_flip_one<T, N>(ctx);
+  test_flip_all<T, N>(ctx);
+  test_none<T, N>(ctx);
+  test_set_reset_one<T, N>(ctx);
+  test_set_all<T, N>(ctx);
+  test_reset_all<T, N>(ctx);
+  test_left_shift<T, N>(ctx);
+  test_right_shift<T, N>(ctx);
+  test_left_rotate<T, N>(ctx);
+  test_right_rotate<T, N>(ctx);
+  test_find_first<T, N>(ctx);
+  test_gray_code<T, N>(ctx);
+  test_comparison<T, N>(ctx);
+  test_iteration<T, N>(ctx);
 }
 
 int
 main()
 {
-	Context ctx;
-
-	// test<chilbert::SmallBitVec, 0>(ctx);
-	test<chilbert::SmallBitVec, 1>(ctx);
-	test<chilbert::SmallBitVec, 31>(ctx);
-	test<chilbert::SmallBitVec, 32>(ctx);
-	test<chilbert::SmallBitVec, 33>(ctx);
-	test<chilbert::SmallBitVec, 63>(ctx);
-	test<chilbert::SmallBitVec, 64>(ctx);
-
-	test<chilbert::DynamicBitVec, 0>(ctx);
-	test<chilbert::DynamicBitVec, 1>(ctx);
-	test<chilbert::DynamicBitVec, 31>(ctx);
-	test<chilbert::DynamicBitVec, 32>(ctx);
-	test<chilbert::DynamicBitVec, 33>(ctx);
-	test<chilbert::DynamicBitVec, 63>(ctx);
-	test<chilbert::DynamicBitVec, 64>(ctx);
-	test<chilbert::DynamicBitVec, 65>(ctx);
-	test<chilbert::DynamicBitVec, 997>(ctx);
-
-	test<chilbert::StaticBitVec<0>, 0>(ctx);
-	test<chilbert::StaticBitVec<1>, 1>(ctx);
-	test<chilbert::StaticBitVec<31>, 31>(ctx);
-	test<chilbert::StaticBitVec<32>, 32>(ctx);
-	test<chilbert::StaticBitVec<33>, 33>(ctx);
-	test<chilbert::StaticBitVec<63>, 63>(ctx);
-	test<chilbert::StaticBitVec<64>, 64>(ctx);
-	test<chilbert::StaticBitVec<65>, 65>(ctx);
-	test<chilbert::StaticBitVec<997>, 997>(ctx);
-
-	test<chilbert::BoundedBitVec<0>, 0>(ctx);
-	test<chilbert::BoundedBitVec<1>, 1>(ctx);
-	test<chilbert::BoundedBitVec<31>, 31>(ctx);
-	test<chilbert::BoundedBitVec<32>, 32>(ctx);
-	test<chilbert::BoundedBitVec<64>, 33>(ctx);
-	test<chilbert::BoundedBitVec<64>, 63>(ctx);
-	test<chilbert::BoundedBitVec<64>, 64>(ctx);
-	test<chilbert::BoundedBitVec<128>, 65>(ctx);
-	test<chilbert::BoundedBitVec<997>, 997>(ctx);
-	test<chilbert::BoundedBitVec<2048>, 997>(ctx);
-
-	return 0;
+  Context ctx;
+
+  // test<chilbert::SmallBitVec, 0>(ctx);
+  test<chilbert::SmallBitVec, 1>(ctx);
+  test<chilbert::SmallBitVec, 31>(ctx);
+  test<chilbert::SmallBitVec, 32>(ctx);
+  test<chilbert::SmallBitVec, 33>(ctx);
+  test<chilbert::SmallBitVec, 63>(ctx);
+  test<chilbert::SmallBitVec, 64>(ctx);
+
+  test<chilbert::DynamicBitVec, 0>(ctx);
+  test<chilbert::DynamicBitVec, 1>(ctx);
+  test<chilbert::DynamicBitVec, 31>(ctx);
+  test<chilbert::DynamicBitVec, 32>(ctx);
+  test<chilbert::DynamicBitVec, 33>(ctx);
+  test<chilbert::DynamicBitVec, 63>(ctx);
+  test<chilbert::DynamicBitVec, 64>(ctx);
+  test<chilbert::DynamicBitVec, 65>(ctx);
+  test<chilbert::DynamicBitVec, 997>(ctx);
+
+  test<chilbert::StaticBitVec<0>, 0>(ctx);
+  test<chilbert::StaticBitVec<1>, 1>(ctx);
+  test<chilbert::StaticBitVec<31>, 31>(ctx);
+  test<chilbert::StaticBitVec<32>, 32>(ctx);
+  test<chilbert::StaticBitVec<33>, 33>(ctx);
+  test<chilbert::StaticBitVec<63>, 63>(ctx);
+  test<chilbert::StaticBitVec<64>, 64>(ctx);
+  test<chilbert::StaticBitVec<65>, 65>(ctx);
+  test<chilbert::StaticBitVec<997>, 997>(ctx);
+
+  test<chilbert::BoundedBitVec<0>, 0>(ctx);
+  test<chilbert::BoundedBitVec<1>, 1>(ctx);
+  test<chilbert::BoundedBitVec<31>, 31>(ctx);
+  test<chilbert::BoundedBitVec<32>, 32>(ctx);
+  test<chilbert::BoundedBitVec<64>, 33>(ctx);
+  test<chilbert::BoundedBitVec<64>, 63>(ctx);
+  test<chilbert::BoundedBitVec<64>, 64>(ctx);
+  test<chilbert::BoundedBitVec<128>, 65>(ctx);
+  test<chilbert::BoundedBitVec<997>, 997>(ctx);
+  test<chilbert::BoundedBitVec<2048>, 997>(ctx);
+
+  return 0;
 }
diff --git a/test/test_gray_code_rank.cpp b/test/test_gray_code_rank.cpp
index 3b6f9bb..906b26c 100644
--- a/test/test_gray_code_rank.cpp
+++ b/test/test_gray_code_rank.cpp
@@ -33,132 +33,132 @@
 #include <cstddef>
 #include <iterator>
 
-template <class T, size_t Max, size_t D>
+template<class T, size_t Max, size_t D>
 T
 get_mask(const std::array<size_t, D>& ms, const size_t d, const size_t step)
 {
-	T      mask = make_zero_bitvec<T, D>();
-	size_t b    = 0U;
-	chilbert::detail::extract_mask(ms.data(), D, d, step, mask, b);
+  T      mask = make_zero_bitvec<T, D>();
+  size_t b    = 0U;
+  chilbert::detail::extract_mask(ms.data(), D, d, step, mask, b);
 
-	assert(b == mask.count());
+  assert(b == mask.count());
 
-	return mask;
+  return mask;
 }
 
-template <class T, size_t Max, size_t D>
+template<class T, size_t Max, size_t D>
 void
 test_extract_mask(Context& ctx)
 {
-	for (size_t d = 0; d < D; ++d) {
-		const auto   ms    = make_random_precisions<Max, D>(ctx);
-		const auto   max_m = *std::max_element(std::begin(ms), std::end(ms));
-		const size_t step  = rand_between(ctx, 0, max_m);
-		const auto   mask  = get_mask<T, D>(ms, d, step);
-
-		for (size_t i = 0; i < D; ++i) {
-			assert(mask.test(i) == (ms[(d + i) % D] > step));
-		}
-	}
+  for (size_t d = 0; d < D; ++d) {
+    const auto   ms    = make_random_precisions<Max, D>(ctx);
+    const auto   max_m = *std::max_element(std::begin(ms), std::end(ms));
+    const size_t step  = rand_between(ctx, 0, max_m);
+    const auto   mask  = get_mask<T, D>(ms, d, step);
+
+    for (size_t i = 0; i < D; ++i) {
+      assert(mask.test(i) == (ms[(d + i) % D] > step));
+    }
+  }
 }
 
-template <class T, size_t Max, size_t D>
+template<class T, size_t Max, size_t D>
 void
 test_gray_code_rank(Context& ctx)
 {
-	for (size_t d = 0; d < D; ++d) {
-		// Generate random mask
-		const auto mask = make_random_bitvec<T, D>(ctx);
-
-		// Generate two random values and their gray codes
-		const auto a = make_random_bitvec<T, D>(ctx);
-		const auto b = make_random_bitvec<T, D>(ctx);
-
-		auto ga = a;
-		chilbert::detail::gray_code(ga);
-
-		auto gb = b;
-		chilbert::detail::gray_code(gb);
-
-		// Calculate gray code ranks
-		auto ra = make_zero_bitvec<T, D>();
-		chilbert::detail::gray_code_rank(mask, ga, D, ra);
-
-		auto rb = make_zero_bitvec<T, D>();
-		chilbert::detail::gray_code_rank(mask, gb, D, rb);
-
-		// Ensure ranks have at most mask.count() bits
-		auto max = make_zero_bitvec<T, D>();
-		chilbert::detail::set_bit(max, mask.count(), 1);
-		assert(ra < max);
-		assert(rb < max);
-
-		// Test fundamental property of gray code ranks
-		const auto mga = ga & mask;
-		const auto mgb = gb & mask;
-		assert((mga < mgb) == (ra < rb));
-
-		// Test inversion
-		const auto pat     = ~mask;
-		auto       ga_out  = make_zero_bitvec<T, D>();
-		auto       gag_out = make_zero_bitvec<T, D>();
-		chilbert::detail::gray_code_rank_inv(
-		  mask, pat, ra, D, mask.count(), gag_out, ga_out);
-		assert((ga_out & mask) == (ga & mask));
-
-		auto gag_check = ga_out;
-		chilbert::detail::gray_code(gag_check);
-		assert(gag_check == gag_out);
-	}
+  for (size_t d = 0; d < D; ++d) {
+    // Generate random mask
+    const auto mask = make_random_bitvec<T, D>(ctx);
+
+    // Generate two random values and their gray codes
+    const auto a = make_random_bitvec<T, D>(ctx);
+    const auto b = make_random_bitvec<T, D>(ctx);
+
+    auto ga = a;
+    chilbert::detail::gray_code(ga);
+
+    auto gb = b;
+    chilbert::detail::gray_code(gb);
+
+    // Calculate gray code ranks
+    auto ra = make_zero_bitvec<T, D>();
+    chilbert::detail::gray_code_rank(mask, ga, D, ra);
+
+    auto rb = make_zero_bitvec<T, D>();
+    chilbert::detail::gray_code_rank(mask, gb, D, rb);
+
+    // Ensure ranks have at most mask.count() bits
+    auto max = make_zero_bitvec<T, D>();
+    chilbert::detail::set_bit(max, mask.count(), 1);
+    assert(ra < max);
+    assert(rb < max);
+
+    // Test fundamental property of gray code ranks
+    const auto mga = ga & mask;
+    const auto mgb = gb & mask;
+    assert((mga < mgb) == (ra < rb));
+
+    // Test inversion
+    const auto pat     = ~mask;
+    auto       ga_out  = make_zero_bitvec<T, D>();
+    auto       gag_out = make_zero_bitvec<T, D>();
+    chilbert::detail::gray_code_rank_inv(
+      mask, pat, ra, D, mask.count(), gag_out, ga_out);
+    assert((ga_out & mask) == (ga & mask));
+
+    auto gag_check = ga_out;
+    chilbert::detail::gray_code(gag_check);
+    assert(gag_check == gag_out);
+  }
 }
 
-template <class T, size_t Max, size_t D>
+template<class T, size_t Max, size_t D>
 void
 test(Context& ctx)
 {
-	test_extract_mask<T, Max, D>(ctx);
-	test_gray_code_rank<T, Max, D>(ctx);
+  test_extract_mask<T, Max, D>(ctx);
+  test_gray_code_rank<T, Max, D>(ctx);
 }
 
 int
 main()
 {
-	Context ctx;
-
-	test<chilbert::SmallBitVec, 64, 1>(ctx);
-	test<chilbert::SmallBitVec, 64, 31>(ctx);
-	test<chilbert::SmallBitVec, 64, 32>(ctx);
-	test<chilbert::SmallBitVec, 64, 33>(ctx);
-	test<chilbert::SmallBitVec, 64, 60>(ctx);
-	test<chilbert::SmallBitVec, 64, 64>(ctx);
-
-	test<chilbert::DynamicBitVec, 64, 1>(ctx);
-	test<chilbert::DynamicBitVec, 64, 31>(ctx);
-	test<chilbert::DynamicBitVec, 64, 32>(ctx);
-	test<chilbert::DynamicBitVec, 64, 33>(ctx);
-	test<chilbert::DynamicBitVec, 64, 60>(ctx);
-	test<chilbert::DynamicBitVec, 64, 64>(ctx);
-	test<chilbert::DynamicBitVec, 96, 65>(ctx);
-	test<chilbert::DynamicBitVec, 1024, 997>(ctx);
-
-	test<chilbert::StaticBitVec<1>, 64, 1>(ctx);
-	test<chilbert::StaticBitVec<31>, 64, 31>(ctx);
-	test<chilbert::StaticBitVec<32>, 64, 32>(ctx);
-	test<chilbert::StaticBitVec<33>, 64, 33>(ctx);
-	test<chilbert::StaticBitVec<60>, 64, 60>(ctx);
-	test<chilbert::StaticBitVec<64>, 64, 64>(ctx);
-	test<chilbert::StaticBitVec<65>, 96, 65>(ctx);
-	test<chilbert::StaticBitVec<997>, 1024, 997>(ctx);
-
-	test<chilbert::BoundedBitVec<1>, 64, 1>(ctx);
-	test<chilbert::BoundedBitVec<31>, 64, 31>(ctx);
-	test<chilbert::BoundedBitVec<32>, 64, 32>(ctx);
-	test<chilbert::BoundedBitVec<64>, 64, 33>(ctx);
-	test<chilbert::BoundedBitVec<64>, 64, 60>(ctx);
-	test<chilbert::BoundedBitVec<64>, 64, 64>(ctx);
-	test<chilbert::BoundedBitVec<128>, 96, 65>(ctx);
-	test<chilbert::BoundedBitVec<997>, 1024, 997>(ctx);
-	test<chilbert::BoundedBitVec<2048>, 1024, 997>(ctx);
-
-	return 0;
+  Context ctx;
+
+  test<chilbert::SmallBitVec, 64, 1>(ctx);
+  test<chilbert::SmallBitVec, 64, 31>(ctx);
+  test<chilbert::SmallBitVec, 64, 32>(ctx);
+  test<chilbert::SmallBitVec, 64, 33>(ctx);
+  test<chilbert::SmallBitVec, 64, 60>(ctx);
+  test<chilbert::SmallBitVec, 64, 64>(ctx);
+
+  test<chilbert::DynamicBitVec, 64, 1>(ctx);
+  test<chilbert::DynamicBitVec, 64, 31>(ctx);
+  test<chilbert::DynamicBitVec, 64, 32>(ctx);
+  test<chilbert::DynamicBitVec, 64, 33>(ctx);
+  test<chilbert::DynamicBitVec, 64, 60>(ctx);
+  test<chilbert::DynamicBitVec, 64, 64>(ctx);
+  test<chilbert::DynamicBitVec, 96, 65>(ctx);
+  test<chilbert::DynamicBitVec, 1024, 997>(ctx);
+
+  test<chilbert::StaticBitVec<1>, 64, 1>(ctx);
+  test<chilbert::StaticBitVec<31>, 64, 31>(ctx);
+  test<chilbert::StaticBitVec<32>, 64, 32>(ctx);
+  test<chilbert::StaticBitVec<33>, 64, 33>(ctx);
+  test<chilbert::StaticBitVec<60>, 64, 60>(ctx);
+  test<chilbert::StaticBitVec<64>, 64, 64>(ctx);
+  test<chilbert::StaticBitVec<65>, 96, 65>(ctx);
+  test<chilbert::StaticBitVec<997>, 1024, 997>(ctx);
+
+  test<chilbert::BoundedBitVec<1>, 64, 1>(ctx);
+  test<chilbert::BoundedBitVec<31>, 64, 31>(ctx);
+  test<chilbert::BoundedBitVec<32>, 64, 32>(ctx);
+  test<chilbert::BoundedBitVec<64>, 64, 33>(ctx);
+  test<chilbert::BoundedBitVec<64>, 64, 60>(ctx);
+  test<chilbert::BoundedBitVec<64>, 64, 64>(ctx);
+  test<chilbert::BoundedBitVec<128>, 96, 65>(ctx);
+  test<chilbert::BoundedBitVec<997>, 1024, 997>(ctx);
+  test<chilbert::BoundedBitVec<2048>, 1024, 997>(ctx);
+
+  return 0;
 }
diff --git a/test/test_hilbert.cpp b/test/test_hilbert.cpp
index a0a64c0..1ed6f8c 100644
--- a/test/test_hilbert.cpp
+++ b/test/test_hilbert.cpp
@@ -35,191 +35,191 @@
 #include <cstdint>
 
 /// Return a `D`-dimensional point within `ms` per-dimension precision
-template <size_t D>
+template<size_t D>
 std::array<uint64_t, D>
 make_random_point(Context& ctx, const std::array<size_t, D>& ms)
 {
-	std::array<uint64_t, D> p{};
-	for (size_t i = 0; i < D; ++i) {
-		p[i] = rand_between(ctx, 0, (1UL << ms[i]) - 1);
-	}
-	return p;
+  std::array<uint64_t, D> p{};
+  for (size_t i = 0; i < D; ++i) {
+    p[i] = rand_between(ctx, 0, (1UL << ms[i]) - 1);
+  }
+  return p;
 }
 
 /// Return the squared distance from point `a` to point `b`
-template <class T, size_t D>
+template<class T, size_t D>
 T
 squared_distance(const std::array<T, D>& a, const std::array<T, D>& b)
 {
-	T sdist = 0;
-	for (size_t i = 0; i < D; ++i) {
-		const T diff = a[i] > b[i] ? a[i] - b[i] : b[i] - a[i];
-		sdist += diff * diff;
-	}
-	return sdist;
+  T sdist = 0;
+  for (size_t i = 0; i < D; ++i) {
+    const T diff = a[i] > b[i] ? a[i] - b[i] : b[i] - a[i];
+    sdist += diff * diff;
+  }
+  return sdist;
 }
 
 /// Convert bit vector `vec` to a big integer
-template <class T>
+template<class T>
 mpz_class
 to_big_int(const T& vec)
 {
-	using Rack = typename T::Rack;
-
-	mpz_t ia;
-	mpz_init(ia);
-	mpz_import(ia, vec.num_racks(), -1, sizeof(Rack), 0, 0, vec.data());
-	const mpz_class num(ia);
-	mpz_clear(ia);
-	return num;
+  using Rack = typename T::Rack;
+
+  mpz_t ia;
+  mpz_init(ia);
+  mpz_import(ia, vec.num_racks(), -1, sizeof(Rack), 0, 0, vec.data());
+  const mpz_class num(ia);
+  mpz_clear(ia);
+  return num;
 }
 
 /// Convert big integer `num` to a bit vector
-template <class T, size_t M>
+template<class T, size_t M>
 T
 from_big_int(const mpz_class& num)
 {
-	using Rack = typename T::Rack;
+  using Rack = typename T::Rack;
 
-	T      vec   = make_zero_bitvec<T, M>();
-	size_t count = 0;
-	mpz_export(vec.data(), &count, -1, sizeof(Rack), 0, 0, num.get_mpz_t());
-	assert(count <= static_cast<size_t>(vec.num_racks()));
-	return vec;
+  T      vec   = make_zero_bitvec<T, M>();
+  size_t count = 0;
+  mpz_export(vec.data(), &count, -1, sizeof(Rack), 0, 0, num.get_mpz_t());
+  assert(count <= static_cast<size_t>(vec.num_racks()));
+  return vec;
 }
 
-template <class H, size_t M, size_t D>
+template<class H, size_t M, size_t D>
 void
 test_standard(Context& ctx)
 {
-	static_assert(M < sizeof(typename H::Rack) * CHAR_BIT, "");
+  static_assert(M < sizeof(typename H::Rack) * CHAR_BIT, "");
 
-	// Generate random point and its hilbert index
-	const auto pa = make_random_point<M, D>(ctx);
+  // Generate random point and its hilbert index
+  const auto pa = make_random_point<M, D>(ctx);
 
-	H ha = make_zero_bitvec<H, D * M>();
-	assert(ha.size() >= D * M);
-	chilbert::coords_to_index(pa, M, D, ha);
+  H ha = make_zero_bitvec<H, D * M>();
+  assert(ha.size() >= D * M);
+  chilbert::coords_to_index(pa, M, D, ha);
 
-	{
-		// Ensure unmapping results in the original point
-		auto pa_out = make_random_point<M, D>(ctx);
-		chilbert::index_to_coords(pa_out, M, D, ha);
-		assert(pa_out == pa);
-	}
+  {
+    // Ensure unmapping results in the original point
+    auto pa_out = make_random_point<M, D>(ctx);
+    chilbert::index_to_coords(pa_out, M, D, ha);
+    assert(pa_out == pa);
+  }
 
-	// Convert hilbert indices to a big integer for manipulation/comparison
-	const auto ia = to_big_int(ha);
+  // Convert hilbert indices to a big integer for manipulation/comparison
+  const auto ia = to_big_int(ha);
 
-	// Generate the next hilbert index
-	const auto ib = ia + 1;
-	const auto hb = from_big_int<H, D * M>(ib);
+  // Generate the next hilbert index
+  const auto ib = ia + 1;
+  const auto hb = from_big_int<H, D * M>(ib);
 
-	// Unmap next hilbert index to a point
-	auto pb = make_random_point<M, D>(ctx);
-	chilbert::index_to_coords(pb, M, D, hb);
+  // Unmap next hilbert index to a point
+  auto pb = make_random_point<M, D>(ctx);
+  chilbert::index_to_coords(pb, M, D, hb);
 
-	// Ensure next point is 1 unit of distance away from first
-	assert(squared_distance(pa, pb) == 1);
+  // Ensure next point is 1 unit of distance away from first
+  assert(squared_distance(pa, pb) == 1);
 }
 
-template <class T, size_t M, size_t D>
+template<class T, size_t M, size_t D>
 void
 test_compact(Context& ctx)
 {
-	static_assert(M < sizeof(typename T::Rack) * CHAR_BIT, "");
+  static_assert(M < sizeof(typename T::Rack) * CHAR_BIT, "");
 
-	// Generate random point and its hilbert index
-	const auto ms = make_random_precisions<D * M, D>(ctx);
-	const auto pa = make_random_point<D>(ctx, ms);
+  // Generate random point and its hilbert index
+  const auto ms = make_random_precisions<D * M, D>(ctx);
+  const auto pa = make_random_point<D>(ctx, ms);
 
-	T ha = make_zero_bitvec<T, D * M>();
-	assert(ha.size() >= D * M);
-	chilbert::coords_to_compact_index(pa, ms.data(), D, ha);
+  T ha = make_zero_bitvec<T, D * M>();
+  assert(ha.size() >= D * M);
+  chilbert::coords_to_compact_index(pa, ms.data(), D, ha);
 
-	{
-		// Ensure unmapping results in the original point
-		auto pa_out = make_random_point<M, D>(ctx);
-		chilbert::compact_index_to_coords(pa_out, ms.data(), D, ha);
-		assert(pa_out == pa);
-	}
+  {
+    // Ensure unmapping results in the original point
+    auto pa_out = make_random_point<M, D>(ctx);
+    chilbert::compact_index_to_coords(pa_out, ms.data(), D, ha);
+    assert(pa_out == pa);
+  }
 
-	// Convert hilbert indices to a big integer for manipulation/comparison
-	const auto ia = to_big_int(ha);
+  // Convert hilbert indices to a big integer for manipulation/comparison
+  const auto ia = to_big_int(ha);
 
-	// Generate the next hilbert index
-	const auto ib = ia + 1;
-	const auto hb = from_big_int<T, D * M>(ib);
+  // Generate the next hilbert index
+  const auto ib = ia + 1;
+  const auto hb = from_big_int<T, D * M>(ib);
 
-	// Unmap next hilbert index to a point
-	auto pb = make_random_point<M, D>(ctx);
-	chilbert::compact_index_to_coords(pb, ms.data(), D, hb);
+  // Unmap next hilbert index to a point
+  auto pb = make_random_point<M, D>(ctx);
+  chilbert::compact_index_to_coords(pb, ms.data(), D, hb);
 
-	// Ensure next point is 1 unit of distance away from first
-	assert(squared_distance(pa, pb) == 1);
+  // Ensure next point is 1 unit of distance away from first
+  assert(squared_distance(pa, pb) == 1);
 }
 
 int
 main()
 {
-	Context ctx;
-
-	test_standard<chilbert::SmallBitVec, 4, 2>(ctx);
-	test_standard<chilbert::SmallBitVec, 32, 2>(ctx);
-	test_standard<chilbert::SmallBitVec, 16, 4>(ctx);
-	test_standard<chilbert::SmallBitVec, 8, 8>(ctx);
-	test_standard<chilbert::SmallBitVec, 4, 16>(ctx);
-	test_standard<chilbert::SmallBitVec, 2, 32>(ctx);
-	test_standard<chilbert::SmallBitVec, 1, 64>(ctx);
-
-	test_standard<chilbert::DynamicBitVec, 4, 65>(ctx);
-	test_standard<chilbert::DynamicBitVec, 32, 64>(ctx);
-	test_standard<chilbert::DynamicBitVec, 63, 128>(ctx);
-
-	test_standard<chilbert::StaticBitVec<4 * 2>, 4, 2>(ctx);
-	test_standard<chilbert::StaticBitVec<32 * 2>, 32, 2>(ctx);
-	test_standard<chilbert::StaticBitVec<16 * 4>, 16, 4>(ctx);
-	test_standard<chilbert::StaticBitVec<8 * 8>, 8, 8>(ctx);
-	test_standard<chilbert::StaticBitVec<4 * 16>, 4, 16>(ctx);
-	test_standard<chilbert::StaticBitVec<2 * 32>, 2, 32>(ctx);
-	test_standard<chilbert::StaticBitVec<1 * 64>, 1, 64>(ctx);
-	test_standard<chilbert::StaticBitVec<4 * 65>, 4, 65>(ctx);
-	test_standard<chilbert::StaticBitVec<32 * 64>, 32, 64>(ctx);
-	test_standard<chilbert::StaticBitVec<63 * 128>, 63, 128>(ctx);
-
-	test_standard<chilbert::BoundedBitVec<4 * 2>, 4, 2>(ctx);
-	test_standard<chilbert::BoundedBitVec<32 * 2>, 32, 2>(ctx);
-	test_standard<chilbert::BoundedBitVec<16 * 4>, 16, 4>(ctx);
-	test_standard<chilbert::BoundedBitVec<8 * 8>, 8, 8>(ctx);
-	test_standard<chilbert::BoundedBitVec<4 * 16>, 4, 16>(ctx);
-	test_standard<chilbert::BoundedBitVec<2 * 32>, 2, 32>(ctx);
-	test_standard<chilbert::BoundedBitVec<1 * 64>, 1, 64>(ctx);
-	test_standard<chilbert::BoundedBitVec<4 * 128>, 4, 65>(ctx);
-	test_standard<chilbert::BoundedBitVec<32 * 128>, 32, 64>(ctx);
-	test_standard<chilbert::BoundedBitVec<63 * 128>, 63, 128>(ctx);
-
-	test_compact<chilbert::SmallBitVec, 4, 2>(ctx);
-	test_compact<chilbert::SmallBitVec, 32, 2>(ctx);
-	test_compact<chilbert::SmallBitVec, 16, 4>(ctx);
-	test_compact<chilbert::SmallBitVec, 8, 8>(ctx);
-	test_compact<chilbert::SmallBitVec, 4, 16>(ctx);
-	test_compact<chilbert::SmallBitVec, 2, 32>(ctx);
-	test_compact<chilbert::SmallBitVec, 1, 64>(ctx);
-
-	test_compact<chilbert::DynamicBitVec, 4, 65>(ctx);
-	test_compact<chilbert::DynamicBitVec, 32, 64>(ctx);
-	test_compact<chilbert::DynamicBitVec, 63, 128>(ctx);
-
-	test_compact<chilbert::StaticBitVec<4 * 2>, 4, 2>(ctx);
-	test_compact<chilbert::StaticBitVec<32 * 2>, 32, 2>(ctx);
-	test_compact<chilbert::StaticBitVec<16 * 4>, 16, 4>(ctx);
-	test_compact<chilbert::StaticBitVec<8 * 8>, 8, 8>(ctx);
-	test_compact<chilbert::StaticBitVec<4 * 16>, 4, 16>(ctx);
-	test_compact<chilbert::StaticBitVec<2 * 32>, 2, 32>(ctx);
-	test_compact<chilbert::StaticBitVec<1 * 64>, 1, 64>(ctx);
-	test_compact<chilbert::StaticBitVec<4 * 65>, 4, 65>(ctx);
-	test_compact<chilbert::StaticBitVec<32 * 64>, 32, 64>(ctx);
-	test_compact<chilbert::StaticBitVec<63 * 128>, 63, 128>(ctx);
-
-	return 0;
+  Context ctx;
+
+  test_standard<chilbert::SmallBitVec, 4, 2>(ctx);
+  test_standard<chilbert::SmallBitVec, 32, 2>(ctx);
+  test_standard<chilbert::SmallBitVec, 16, 4>(ctx);
+  test_standard<chilbert::SmallBitVec, 8, 8>(ctx);
+  test_standard<chilbert::SmallBitVec, 4, 16>(ctx);
+  test_standard<chilbert::SmallBitVec, 2, 32>(ctx);
+  test_standard<chilbert::SmallBitVec, 1, 64>(ctx);
+
+  test_standard<chilbert::DynamicBitVec, 4, 65>(ctx);
+  test_standard<chilbert::DynamicBitVec, 32, 64>(ctx);
+  test_standard<chilbert::DynamicBitVec, 63, 128>(ctx);
+
+  test_standard<chilbert::StaticBitVec<4 * 2>, 4, 2>(ctx);
+  test_standard<chilbert::StaticBitVec<32 * 2>, 32, 2>(ctx);
+  test_standard<chilbert::StaticBitVec<16 * 4>, 16, 4>(ctx);
+  test_standard<chilbert::StaticBitVec<8 * 8>, 8, 8>(ctx);
+  test_standard<chilbert::StaticBitVec<4 * 16>, 4, 16>(ctx);
+  test_standard<chilbert::StaticBitVec<2 * 32>, 2, 32>(ctx);
+  test_standard<chilbert::StaticBitVec<1 * 64>, 1, 64>(ctx);
+  test_standard<chilbert::StaticBitVec<4 * 65>, 4, 65>(ctx);
+  test_standard<chilbert::StaticBitVec<32 * 64>, 32, 64>(ctx);
+  test_standard<chilbert::StaticBitVec<63 * 128>, 63, 128>(ctx);
+
+  test_standard<chilbert::BoundedBitVec<4 * 2>, 4, 2>(ctx);
+  test_standard<chilbert::BoundedBitVec<32 * 2>, 32, 2>(ctx);
+  test_standard<chilbert::BoundedBitVec<16 * 4>, 16, 4>(ctx);
+  test_standard<chilbert::BoundedBitVec<8 * 8>, 8, 8>(ctx);
+  test_standard<chilbert::BoundedBitVec<4 * 16>, 4, 16>(ctx);
+  test_standard<chilbert::BoundedBitVec<2 * 32>, 2, 32>(ctx);
+  test_standard<chilbert::BoundedBitVec<1 * 64>, 1, 64>(ctx);
+  test_standard<chilbert::BoundedBitVec<4 * 128>, 4, 65>(ctx);
+  test_standard<chilbert::BoundedBitVec<32 * 128>, 32, 64>(ctx);
+  test_standard<chilbert::BoundedBitVec<63 * 128>, 63, 128>(ctx);
+
+  test_compact<chilbert::SmallBitVec, 4, 2>(ctx);
+  test_compact<chilbert::SmallBitVec, 32, 2>(ctx);
+  test_compact<chilbert::SmallBitVec, 16, 4>(ctx);
+  test_compact<chilbert::SmallBitVec, 8, 8>(ctx);
+  test_compact<chilbert::SmallBitVec, 4, 16>(ctx);
+  test_compact<chilbert::SmallBitVec, 2, 32>(ctx);
+  test_compact<chilbert::SmallBitVec, 1, 64>(ctx);
+
+  test_compact<chilbert::DynamicBitVec, 4, 65>(ctx);
+  test_compact<chilbert::DynamicBitVec, 32, 64>(ctx);
+  test_compact<chilbert::DynamicBitVec, 63, 128>(ctx);
+
+  test_compact<chilbert::StaticBitVec<4 * 2>, 4, 2>(ctx);
+  test_compact<chilbert::StaticBitVec<32 * 2>, 32, 2>(ctx);
+  test_compact<chilbert::StaticBitVec<16 * 4>, 16, 4>(ctx);
+  test_compact<chilbert::StaticBitVec<8 * 8>, 8, 8>(ctx);
+  test_compact<chilbert::StaticBitVec<4 * 16>, 4, 16>(ctx);
+  test_compact<chilbert::StaticBitVec<2 * 32>, 2, 32>(ctx);
+  test_compact<chilbert::StaticBitVec<1 * 64>, 1, 64>(ctx);
+  test_compact<chilbert::StaticBitVec<4 * 65>, 4, 65>(ctx);
+  test_compact<chilbert::StaticBitVec<32 * 64>, 32, 64>(ctx);
+  test_compact<chilbert::StaticBitVec<63 * 128>, 63, 128>(ctx);
+
+  return 0;
 }
diff --git a/test/test_utils.hpp b/test/test_utils.hpp
index 65bd1b1..c3b6edb 100644
--- a/test/test_utils.hpp
+++ b/test/test_utils.hpp
@@ -26,70 +26,69 @@
 #include <random>
 
 /// Test context
-struct Context
-{
-	std::random_device                    rng;
-	std::uniform_int_distribution<size_t> dist{0, SIZE_MAX};
+struct Context {
+  std::random_device                    rng;
+  std::uniform_int_distribution<size_t> dist{0, SIZE_MAX};
 };
 
 /// Return a bit vector of type T with N zero bits
-template <class T, size_t N>
+template<class T, size_t N>
 T
 make_zero_bitvec()
 {
-	T v(N);
-	v.reset();
-	return v;
+  T v(N);
+  v.reset();
+  return v;
 }
 
 /// Return a bit vector of type T with N random bits
-template <class T, size_t N>
+template<class T, size_t N>
 T
 make_random_bitvec(Context& ctx)
 {
-	T v(N);
-	for (size_t i = 0; i < N; ++i) {
-		v.set(i, ctx.dist(ctx.rng) & 1U);
-	}
-	return v;
+  T v(N);
+  for (size_t i = 0; i < N; ++i) {
+    v.set(i, ctx.dist(ctx.rng) & 1U);
+  }
+  return v;
 }
 
 /// Return a random number in [min, max)
 static inline size_t
 rand_between(Context& ctx, const size_t min, const size_t max)
 {
-	assert(max >= min);
-	const size_t r = (max == min) ? min : ctx.dist(ctx.rng) % (max - min) + min;
-	assert(r >= min && r < max);
-	return r;
+  assert(max >= min);
+  const size_t r = (max == min) ? min : ctx.dist(ctx.rng) % (max - min) + min;
+  assert(r >= min && r < max);
+  return r;
 }
 
 /// Return an array of D precisions each at most Max and with sum at most N
-template <size_t N, size_t D, size_t Max = 64>
+template<size_t N, size_t D, size_t Max = 64>
 std::array<size_t, D>
 make_random_precisions(Context& ctx)
 {
-	std::array<size_t, D> ms{};
+  std::array<size_t, D> ms{};
 
-	size_t bits_left = N;
-	for (size_t i = 0; i < D; ++i) {
-		ms[i] = rand_between(ctx, 1, std::min(Max, bits_left / (D - i) + 1));
-		bits_left -= ms[i];
-	}
+  size_t bits_left = N;
+  for (size_t i = 0; i < D; ++i) {
+    ms[i] = rand_between(ctx, 1, std::min(Max, bits_left / (D - i) + 1));
+    bits_left -= ms[i];
+  }
 
-	return ms;
+  return ms;
 }
 
 /// Return a `D`-dimensional point with `M` bits of precision per dimension
-template <size_t M, size_t D>
+template<size_t M, size_t D>
 std::array<uint64_t, D>
 make_random_point(Context& ctx)
 {
-	std::array<uint64_t, D> p{};
-	for (size_t i = 0; i < D; ++i) {
-		p[i] = rand_between(ctx, 0, (1UL << M) - 1);
-	}
-	return p;
+  std::array<uint64_t, D> p{};
+  for (size_t i = 0; i < D; ++i) {
+    p[i] = rand_between(ctx, 0, (1UL << M) - 1);
+  }
+  return p;
 }
 
 #endif