diff --git a/2022/14.cpp b/2022/14.cpp
index 77f63bd..def2cc7 100644
--- a/2022/14.cpp
+++ b/2022/14.cpp
@@ -1,12 +1,11 @@
 #include <version>
 
-#include <concepts>
 #include <cstdint>
 #include <iostream>
-#include <set>
 #include <sstream>
 #include <stack>
 #include <tuple>
+#include <unordered_set>
 #include <vector>
 
 #include <boost/phoenix.hpp>
@@ -62,9 +61,9 @@ auto Parse(std::istream & in) -> Input
   return result;
 }
 
-auto ExpandLines(Input const& lines) -> std::set<aocpp::Coord>
+auto ExpandLines(Input const& lines) -> std::unordered_set<aocpp::Coord>
 {
-  auto result = std::set<aocpp::Coord>{};
+  auto result = std::unordered_set<aocpp::Coord>{};
   for (auto const& line : lines) {
     for (auto const i : boost::irange(std::size_t{1}, line.size())) {
       auto const& a = line[i-1];
@@ -99,29 +98,35 @@ auto FindBottom(Input const& input) -> std::int64_t
 
 enum class Action { Visit, Mark };
   
-auto Part1(std::int64_t bottom, std::set<aocpp::Coord> world) -> std::size_t
+auto Part1(std::int64_t bottom, std::unordered_set<aocpp::Coord> world) -> std::size_t
 {
   auto const starting_size = world.size();
   auto work = std::stack<std::pair<Action, aocpp::Coord>>{};
   work.emplace(Action::Visit, TOP);
 
   while (!work.empty()) {
-    auto const [action, here] = work.top();
-    work.pop();
-    
+    auto & [action, here] = work.top();
     switch (action) {
       case Action::Visit:
         if (!world.contains(here)) {
+          
+          // Particle is falling off the world, stop simulation
           if (here.y == bottom) {
             return world.size() - starting_size;
           }
-          work.emplace(Action::Mark, here);
-          work.emplace(Action::Visit, aocpp::Right(aocpp::Down(here)));
-          work.emplace(Action::Visit, aocpp::Left(aocpp::Down(here)));
-          work.emplace(Action::Visit, aocpp::Down(here));
+          
+          action = Action::Mark;
+          auto const down = aocpp::Down(here);
+          // action and here are invalidated after this line
+          work.emplace(Action::Visit, aocpp::Right(down));
+          work.emplace(Action::Visit, aocpp::Left(down));
+          work.emplace(Action::Visit, down);
+        } else {
+          work.pop();
         }
         break;
       case Action::Mark:
+        work.pop();
         world.insert(here);
         break;
     }
@@ -130,27 +135,30 @@ auto Part1(std::int64_t bottom, std::set<aocpp::Coord> world) -> std::size_t
   throw std::runtime_error{"No solution for part 1"};
 }
 
-auto Part2(std::int64_t bottom, std::set<aocpp::Coord> world) -> std::size_t
+auto Part2(std::int64_t bottom, std::unordered_set<aocpp::Coord> world) -> std::size_t
 {
   auto const starting_size = world.size();
   auto work = std::stack<std::pair<Action, aocpp::Coord>>{};
   work.emplace(Action::Visit, TOP);
 
   while (!work.empty()) {
-    auto const [action, here] = work.top();
-    work.pop();
+    auto & [action, here] = work.top();
     
     switch (action) {
       case Action::Visit:
         if (here.y - 2 != bottom && !world.contains(here)) {
-          work.emplace(Action::Mark, here);
-          work.emplace(Action::Visit, aocpp::Right(aocpp::Down(here)));
-          work.emplace(Action::Visit, aocpp::Left(aocpp::Down(here)));
-          work.emplace(Action::Visit, aocpp::Down(here));
+          action = Action::Mark;
+          auto const down = aocpp::Down(here);
+          work.emplace(Action::Visit, aocpp::Right(down));
+          work.emplace(Action::Visit, aocpp::Left(down));
+          work.emplace(Action::Visit, down);
+        } else {
+          work.pop();
         }
         break;
       case Action::Mark:
         world.insert(here);
+        work.pop();
         break;
     }
   }
diff --git a/2022/15.cpp b/2022/15.cpp
index dd4c288..48b6fae 100644
--- a/2022/15.cpp
+++ b/2022/15.cpp
@@ -11,8 +11,8 @@
 
 #include <doctest.h>
 
-#include <aocpp/Startup.hpp>
 #include <aocpp/Coord.hpp>
+#include <aocpp/Startup.hpp>
 
 namespace {
 
diff --git a/2022/21.cpp b/2022/21.cpp
new file mode 100644
index 0000000..70a584a
--- /dev/null
+++ b/2022/21.cpp
@@ -0,0 +1,191 @@
+#include <cstdint>
+#include <iostream>
+#include <sstream>
+#include <string>
+#include <unordered_map>
+#include <variant>
+#include <vector>
+
+#include <boost/phoenix.hpp>
+#include <boost/range/irange.hpp>
+#include <boost/spirit/include/qi.hpp>
+
+#include <doctest.h>
+
+#include <aocpp/Overloaded.hpp>
+#include <aocpp/Startup.hpp>
+
+namespace {
+
+namespace phx = boost::phoenix;
+namespace qi = boost::spirit::qi;
+
+enum class Op { Add, Sub, Mul, Div };
+
+struct Expr {
+  std::string lhs;
+  std::string rhs;
+  Op op;
+};
+
+struct Entry {
+  std::string lvalue;
+  std::variant<std::int64_t, Expr> rvalue;
+};
+
+using Input = std::vector<Entry>;
+
+// Input file grammar
+template <typename It>
+class Grammar : public qi::grammar<It, Input()> {
+  qi::rule<It, std::string()> variable;
+  qi::rule<It, Op()> op;
+  qi::rule<It, Expr()> expr;
+  qi::rule<It, std::variant<std::int64_t, Expr>()> rhs;
+  qi::rule<It, Entry()> line;
+  qi::rule<It, Input()> input;
+
+public:
+  Grammar() : Grammar::base_type{input} {
+    using namespace qi::labels;
+
+    variable %= qi::as_string[+qi::alpha];
+    op = qi::string("+") [ _val = Op::Add ]
+       | qi::string("-") [ _val = Op::Sub ]
+       | qi::string("*") [ _val = Op::Mul ]
+       | qi::string("/") [ _val = Op::Div ];
+    expr = variable [phx::bind(&Expr::lhs, _val) = _1]
+        >> " "
+        >> op [phx::bind(&Expr::op, _val) = _1]
+        >> " "
+        >> variable [phx::bind(&Expr::rhs, _val) = _1];
+    rhs = expr [ _val = _1 ] | qi::long_long [ _val = _1 ];
+    line = variable [ phx::bind(&Entry::lvalue, _val) = _1]
+        >> ": "
+        >> rhs [ phx::bind(&Entry::rvalue, _val) = _1]
+        >> "\n";
+    input = *line;
+  }
+};
+
+/// Parse the complete input stream according to the rules defined in 'Grammar'.
+/// Throws: std::runtime_error on failed parse
+auto Parse(std::istream & in) -> Input
+{
+  auto result = Input{};
+  auto const content = std::string{std::istreambuf_iterator{in}, {}};
+  auto b = content.begin(); // updated on successful parse
+  auto const e = content.end();
+
+  if (!qi::parse(b, e, Grammar<decltype(b)>{}, result) || b != e) {
+    throw std::runtime_error{"Bad input file: " + content};
+  }
+
+  return result;
+}
+
+auto Eval(
+  std::unordered_map<std::string, std::int64_t> & values,
+  std::unordered_map<std::string, Expr> const& exprs,
+  std::string const& var)
+  -> std::int64_t
+{
+  auto it = values.find(var);
+  if (it != values.end()) {
+    return it->second;
+  }
+  auto const& [lhs,rhs,op] = exprs.at(var);
+  auto const l = Eval(values, exprs, lhs);
+  auto const r = Eval(values, exprs, rhs);
+  std::int64_t o;
+  switch (op) {
+    case Op::Add: o = l + r; break;
+    case Op::Sub: o = l - r; break;
+    case Op::Mul: o = l * r; break;
+    case Op::Div: o = l / r; break;
+  }
+  values.emplace(var, o);
+  return o;
+}
+
+auto Part1(Input const& input) -> std::int64_t
+{
+  std::unordered_map<std::string, std::int64_t> values;
+  std::unordered_map<std::string, Expr> exprs;
+  for (auto const& entry : input) {
+    std::visit(overloaded {
+      [&](Expr const& expr) { exprs.emplace(entry.lvalue, expr); },
+      [&](std::int64_t val) { values.emplace(entry.lvalue, val); }
+    }, entry.rvalue);
+  }
+  return Eval(values, exprs, "root");
+}
+
+auto Part2(Input const& input) -> std::int64_t
+{
+  std::unordered_map<std::string, std::int64_t> values;
+  std::unordered_map<std::string, Expr> exprs;
+  for (auto const& entry : input) {
+    std::visit(overloaded {
+      [&](Expr const& expr) { exprs.emplace(entry.lvalue, expr); },
+      [&](std::int64_t val) { values.emplace(entry.lvalue, val); }
+    }, entry.rvalue);
+  }
+
+  exprs.at("root").op = Op::Sub;
+
+  auto eval = [&](std::int64_t humn) -> std::int64_t {
+    auto values_ = values;
+    values_.at("humn") = humn;
+    return Eval(values_, exprs, "root");
+  };
+
+  auto x0 = 10;
+  auto x1 = 20;
+  auto fx0 = eval(x0);
+  auto fx1 = eval(x1);
+  std::int64_t a = 0, b = 1;
+  while (x0 != x1 || a%b != 0) {
+    a = eval(x1) * (x1 - x0);
+    b = eval(x1) - eval(x0);
+    auto const x2 = x1 - a / b;
+    x0 = x1;
+    x1 = x2;
+  }
+
+  return x0;
+}
+
+} // namespace
+
+TEST_SUITE("2022-21 examples") {
+  TEST_CASE("documented example") {
+    auto in = std::istringstream{
+      "root: pppw + sjmn\n"
+      "dbpl: 5\n"
+      "cczh: sllz + lgvd\n"
+      "zczc: 2\n"
+      "ptdq: humn - dvpt\n"
+      "dvpt: 3\n"
+      "lfqf: 4\n"
+      "humn: 5\n"
+      "ljgn: 2\n"
+      "sjmn: drzm * dbpl\n"
+      "sllz: 4\n"
+      "pppw: cczh / lfqf\n"
+      "lgvd: ljgn * ptdq\n"
+      "drzm: hmdt - zczc\n"
+      "hmdt: 32\n"
+    };
+    auto const input = Parse(in);
+    CHECK(152 == Part1(input));
+    CHECK(301 == Part2(input));
+  }
+}
+
+auto Main(std::istream & in, std::ostream & out) -> void
+{
+  auto const input = Parse(in);
+  out << "Part 1: " << Part1(input) << std::endl;
+  out << "Part 2: " << Part2(input) << std::endl;
+}
diff --git a/2022/CMakeLists.txt b/2022/CMakeLists.txt
index e14c5ab..6641659 100644
--- a/2022/CMakeLists.txt
+++ b/2022/CMakeLists.txt
@@ -46,5 +46,8 @@ target_link_libraries(2022_18 aocpp)
 add_executable(2022_20 20.cpp)
 target_link_libraries(2022_20 aocpp Boost::headers)
 
+add_executable(2022_21 21.cpp)
+target_link_libraries(2022_21 aocpp Boost::headers)
+
 add_executable(2022_25 25.cpp)
 target_link_libraries(2022_25 aocpp)