diff --git a/.clang-tidy b/.clang-tidy
new file mode 100644
index 0000000000000000000000000000000000000000..c03f3a3f034cfa48cb255e5a3582356f26e543f9
--- /dev/null
+++ b/.clang-tidy
@@ -0,0 +1,22 @@
+Checks: "clang-analyzer-*,
+         cppcoreguidelines-*,
+         misc-*,
+         modernize-*,
+         performance-*,
+         portability-*,
+         readability-*,
+         -cppcoreguidelines-avoid-c-arrays,
+         -cppcoreguidelines-avoid-magic-numbers,
+         -cppcoreguidelines-macro-usage,
+         -cppcoreguidelines-pro-bounds-array-to-pointer-decay,
+         -cppcoreguidelines-pro-bounds-pointer-arithmetic,
+         -cppcoreguidelines-pro-type-union-access,
+         -misc-non-private-member-variables-in-classes,
+         -misc-no-recursion,
+         -modernize-avoid-c-arrays,
+         -modernize-use-trailing-return-type,
+         -readability-function-cognitive-complexity,
+         -readability-magic-numbers,
+         -readability-named-parameter"
+WarningsAsErrors: "*"
+HeaderFilterRegex: "src/"
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 3f8304e375e1442932b0cf79e2855505c74690d3..5d087f8644c61f7e14d89fe3cde7ffc17f12cb87 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,19 +1,24 @@
 cmake_minimum_required(VERSION 3.8)
-
-project(
-	"B-ASIC"
+project("B-ASIC"
 	VERSION 1.0.0
 	DESCRIPTION "Better ASIC Toolbox for Python 3"
-	LANGUAGES C CXX
-)
+	LANGUAGES C CXX)
 
-# Find dependencies.
-find_package(fmt REQUIRED)
-find_package(pybind11 CONFIG REQUIRED)
+option(ASIC_USE_FETCHCONTENT "Automatically download dependencies" ON)
+option(ASIC_USE_CLANG_TIDY "Use clang-tidy for static analysis" OFF)
+option(ASIC_BUILDING_PYTHON_DISTRIBUTION "Don't copy compiled binaries to project directory" OFF)
 
 set(LIBRARY_NAME "b_asic") # Name of the python library directory.
 set(TARGET_NAME "_${LIBRARY_NAME}") # Name of this extension module.
 
+# Find dependencies.
+if(ASIC_USE_FETCHCONTENT)
+	add_subdirectory(dependencies)
+else()
+	find_package(fmt REQUIRED)
+	find_package(pybind11 CONFIG REQUIRED)
+endif()
+
 # Set output directory for compiled binaries.
 if(NOT CMAKE_LIBRARY_OUTPUT_DIRECTORY)
 	include(GNUInstallDirs)
@@ -32,9 +37,7 @@ set(CMAKE_RUNTIME_OUTPUT_DIRECTORY_DEBUG "${CMAKE_LIBRARY_OUTPUT_DIRECTORY}")
 set(CMAKE_RUNTIME_OUTPUT_DIRECTORY_RELEASE "${CMAKE_LIBRARY_OUTPUT_DIRECTORY}")
 
 # Add files to be compiled into Python module.
-pybind11_add_module(
-	"${TARGET_NAME}"
-
+pybind11_add_module("${TARGET_NAME}"
 	# Main files.
 	"${CMAKE_CURRENT_SOURCE_DIR}/src/main.cpp"
 	"${CMAKE_CURRENT_SOURCE_DIR}/src/simulation.cpp"
@@ -53,49 +56,44 @@ pybind11_add_module(
 )
 
 # Include headers.
-target_include_directories(
-	"${TARGET_NAME}"
-    PRIVATE
-        "${CMAKE_CURRENT_SOURCE_DIR}/src"
-)
+target_include_directories("${TARGET_NAME}" PRIVATE "${CMAKE_CURRENT_SOURCE_DIR}/src")
 
 # Use C++17.
-target_compile_features(
-	"${TARGET_NAME}"
-	PRIVATE
-		cxx_std_17
-)
+target_compile_features("${TARGET_NAME}" PRIVATE cxx_std_17)
 
 # Set compiler-specific options using generator expressions.
-target_compile_options(
-	"${TARGET_NAME}"
-	PRIVATE
-		$<$<OR:$<CXX_COMPILER_ID:GNU>,$<CXX_COMPILER_ID:Clang>>:
-			-W -Wall -Wextra -Werror -Wno-psabi
-			$<$<CONFIG:Debug>:-g>
-			$<$<NOT:$<CONFIG:Debug>>:-O3>
-		>
-		$<$<CXX_COMPILER_ID:MSVC>:
-			/W3 /WX /permissive- /utf-8
-			$<$<CONFIG:Debug>:/Od>
-			$<$<NOT:$<CONFIG:Debug>>:/Ot>
-		>
-)
+target_compile_options("${TARGET_NAME}" PRIVATE
+	$<$<CXX_COMPILER_ID:GNU>:   -std=c++17  -Wall -Wextra   -Wpedantic      -Werror -Wno-psabi      $<$<CONFIG:Debug>:-Og -g>   $<$<CONFIG:Release>:-O3>    $<$<CONFIG:MinSizeRel>:-Os> $<$<CONFIG:RelWithDebInfo>:-O3 -g>>
+	$<$<CXX_COMPILER_ID:Clang>: -std=c++17  -Wall -Wextra   -Wpedantic      -Werror                 $<$<CONFIG:Debug>:-Og -g>   $<$<CONFIG:Release>:-O3>    $<$<CONFIG:MinSizeRel>:-Os> $<$<CONFIG:RelWithDebInfo>:-O3 -g>>
+	$<$<CXX_COMPILER_ID:MSVC>:  /std:c++17  /W3             /permissive-    /WX     /wd4996 /utf-8  $<$<CONFIG:Debug>:/Od>      $<$<CONFIG:Release>:/Ot>    $<$<CONFIG:MinSizeRel>:/Os> $<$<CONFIG:RelWithDebInfo>:/Ot /Od>>)
 
 # Add libraries. Note: pybind11 is already added in pybind11_add_module.
-target_link_libraries(
-	"${TARGET_NAME}"
-	PRIVATE
+if(ASIC_USE_FETCHCONTENT)
+	target_link_libraries("${TARGET_NAME}" PRIVATE
+		dependency_fmt)
+else()
+	target_link_libraries("${TARGET_NAME}" PRIVATE
 		$<TARGET_NAME_IF_EXISTS:fmt::fmt-header-only>
-		$<$<NOT:$<TARGET_EXISTS:fmt::fmt-header-only>>:fmt::fmt>
-)
+		$<$<NOT:$<TARGET_EXISTS:fmt::fmt-header-only>>:fmt::fmt>)
+endif()
+
+# Set up clang-tidy.
+if(ASIC_USE_CLANG_TIDY)
+	find_program(CLANG_TIDY NAMES clang-tidy REQUIRED)
+	set_property(TARGET "${TARGET_NAME}" PROPERTY CXX_CLANG_TIDY ${CLANG_TIDY})
+	if(MSVC)
+		set_target_properties("${TARGET_NAME}" PROPERTIES
+			VS_GLOBAL_RunCodeAnalysis true
+			VS_GLOBAL_EnableMicrosoftCodeAnalysis false
+			VS_GLOBAL_EnableClangTidyCodeAnalysis true
+			VS_GLOBAL_ClangTidyToolExe "${CLANG_TIDY}")
+	endif()
+endif()
 
 # Copy binaries to project folder for debugging during development.
 if(NOT ASIC_BUILDING_PYTHON_DISTRIBUTION)
-	add_custom_target(
-		copy_binaries ALL
+	add_custom_target(copy_binaries ALL
 		COMMAND ${CMAKE_COMMAND} -E copy "$<TARGET_FILE:${TARGET_NAME}>" "${CMAKE_CURRENT_LIST_DIR}"
 		COMMENT "Copying binaries to ${CMAKE_CURRENT_LIST_DIR}"
-		DEPENDS "${TARGET_NAME}"
-	)
+		DEPENDS "${TARGET_NAME}")
 endif()
diff --git a/MANIFEST.in b/MANIFEST.in
index 96b265fd04394abbdbf033028b46c98619817869..36ba9fca42e0017fece97fcf50bd793cd372eda8 100644
--- a/MANIFEST.in
+++ b/MANIFEST.in
@@ -2,4 +2,4 @@ include README.md
 include LICENSE
 include CMakeLists.txt
 include b_asic/GUI/operation_icons/*
-recursive-include src *.cpp *.h
+recursive-include src *.cpp *.hpp
diff --git a/README.md b/README.md
index f0743ea53ec24c5b1b62fb7d845b0008f8f45a9d..36bfe17d7ba25dd3794816db13aee0c19dfbd999 100644
--- a/README.md
+++ b/README.md
@@ -1,45 +1,55 @@
 <img src="logo.png" width="278" height="100">
 
 # B-ASIC - Better ASIC Toolbox
+
 B-ASIC is an ASIC toolbox for Python 3 that simplifies circuit design and optimization.
 
 ## Development
+
 How to build and debug the library during development.
 
 ### Prerequisites
+
 The following packages are required in order to build the library:
-* cmake 3.8+
-* gcc 7+/clang 7+/msvc 16+
-* fmtlib
-* pybind11 2.3.0+
-* python 3.6+
-* Python:
-  * graphviz
-  * matplotlib
-  * numpy
-  * pybind11
-  * pyside2
-  * qtpy
-  * scipy
-  * setuptools
+
+-   cmake 3.8+
+-   gcc 7+/clang 7+/msvc 16+
+-   fmtlib
+-   pybind11 2.3.0+
+-   python 3.6+
+-   Python:
+    -   graphviz
+    -   matplotlib
+    -   numpy
+    -   pybind11
+    -   pyside2
+    -   qtpy
+    -   scipy
+    -   setuptools
 
 To build a binary distribution, the following additional packages are required:
-* Python:
-  * wheel
+
+-   Python:
+    -   wheel
 
 To run the test suite, the following additional packages are required:
-* Python:
-  * pytest
-  * pytest-cov (for testing with coverage)
-  
+
+-   Python:
+    -   pytest
+    -   pytest-cov (for testing with coverage)
+
 To generate the documentation, the following additional packages are required:
-* doxygen
+
+-   doxygen
 
 ### Using CMake directly
+
 How to build using CMake.
 
 #### Configuring
+
 In `B-ASIC`:
+
 ```
 mkdir build
 cd build
@@ -47,53 +57,73 @@ cmake ..
 ```
 
 #### Building (Debug)
+
 In `B-ASIC/build`:
+
 ```
 cmake --build .
 ```
+
 The output gets written to `B-ASIC/build/lib`.
 
 #### Building (Release)
+
 In `B-ASIC/build`:
+
 ```
 cmake --build . --config Release
 ```
+
 The output gets written to `B-ASIC/build/lib`.
 
 ### Using setuptools to create a package
+
 How to create a package using setuptools that can be installed using pip.
 
 #### Setup (Binary distribution)
+
 In `B-ASIC`:
+
 ```
 python3 setup.py bdist_wheel
 ```
+
 The output gets written to `B-ASIC/dist/b_asic-<version>-<python_tag>-<abi_tag>-<platform_tag>.whl`.
 
 #### Setup (Source distribution)
+
 In `B-ASIC`:
+
 ```
 python3 setup.py sdist
 ```
+
 The output gets written to `B-ASIC/dist/b-asic-<version>.tar.gz`.
 
 #### Installation (Binary distribution)
+
 In `B-ASIC/dist`:
+
 ```
 pip install b_asic-<version>-<python_tag>-<abi_tag>-<platform_tag>.whl
 ```
 
 #### Installation (Source distribution)
+
 In `B-ASIC/dist`:
+
 ```
 pip install b-asic-<version>.tar.gz
 ```
 
 ### Running tests
+
 How to run the tests using pytest in a virtual environment.
 
 #### Linux/OS X
+
 In `B-ASIC`:
+
 ```
 python3 -m venv env
 source env/bin/activate
@@ -102,7 +132,9 @@ pytest
 ```
 
 #### Windows
+
 In `B-ASIC` (as admin):
+
 ```
 python3 -m venv env
 .\env\Scripts\activate.bat
@@ -111,26 +143,33 @@ pytest
 ```
 
 #### Test with coverage
+
 ```
 pytest --cov=b_asic --cov-report html test
 ```
 
 ### Generating documentation
+
 In `B-ASIC`:
+
 ```
 doxygen
 ```
+
 The output gets written to `B-ASIC/doc`.
 
 ## Usage
+
 How to build and use the library as a user.
 
 ### Installation
+
 ```
 pip install b_asic
 ```
 
 ### Importing
+
 ```
 python3
 >>> import b_asic as asic
@@ -138,5 +177,6 @@ python3
 ```
 
 ## License
+
 B-ASIC is distributed under the MIT license.
 See the included LICENSE file for more information.
diff --git a/dependencies/CMakeLists.txt b/dependencies/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..229a492f5ce8c0f1e232fffa3a06e334d057d686
--- /dev/null
+++ b/dependencies/CMakeLists.txt
@@ -0,0 +1,4 @@
+include(FetchContent)
+
+add_subdirectory(fmt)
+add_subdirectory(pybind11)
diff --git a/dependencies/fmt/CMakeLists.txt b/dependencies/fmt/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..fda9e099476f8f7c29878a61253623ca3e1b8699
--- /dev/null
+++ b/dependencies/fmt/CMakeLists.txt
@@ -0,0 +1,11 @@
+message(STATUS "Fetching fmt...")
+
+FetchContent_Declare(fmt
+	GIT_REPOSITORY https://github.com/fmtlib/fmt
+	GIT_TAG d141cdbeb0fb422a3fb7173b285fd38e0d1772dc # 8.0.1
+)
+FetchContent_MakeAvailable(fmt)
+
+add_library(dependency_fmt INTERFACE)
+target_include_directories(dependency_fmt SYSTEM INTERFACE $<TARGET_PROPERTY:fmt-header-only,INTERFACE_INCLUDE_DIRECTORIES>)
+target_link_libraries(dependency_fmt INTERFACE fmt-header-only)
diff --git a/dependencies/pybind11/CMakeLists.txt b/dependencies/pybind11/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..dea520bc8c89a1793c6b609c1d075a286ed37289
--- /dev/null
+++ b/dependencies/pybind11/CMakeLists.txt
@@ -0,0 +1,7 @@
+message(STATUS "Fetching pybind11...")
+
+FetchContent_Declare(pybind11
+	GIT_REPOSITORY https://github.com/pybind/pybind11
+	GIT_TAG f7b499615e14d70ab098a20deb0cdb3889998a1a # 2.8.1
+)
+FetchContent_MakeAvailable(pybind11)
diff --git a/legacy/README.md b/legacy/README.md
index 746d1efdb57c04f4e52bb2ce7f0a7def99c744ed..a450266c51265e81844e9c15353e99d90c798d8e 100644
--- a/legacy/README.md
+++ b/legacy/README.md
@@ -8,4 +8,4 @@ or to be used as a refererence for future development.
 This folder contains a C++ implementation of the Simulation class designed
 using Object-Oriented Programming, as opposed to the current version that uses
 Data-Oriented Design. They are functionally identical, but use different
-styles of programming and have different performance characteristics.
\ No newline at end of file
+styles of programming and have different performance characteristics.
diff --git a/legacy/simulation_oop/core_operations.hpp b/legacy/simulation_oop/core_operations.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..16b05746d436f3e59c0766879c0cf723d9ff459f
--- /dev/null
+++ b/legacy/simulation_oop/core_operations.hpp
@@ -0,0 +1,235 @@
+#ifndef ASIC_SIMULATION_CORE_OPERATIONS_HPP
+#define ASIC_SIMULATION_CORE_OPERATIONS_HPP
+
+#include "../debug.hpp"
+#include "../number.hpp"
+#include "operation.hpp"
+
+#include <algorithm>
+#include <cmath>
+#include <cstddef>
+#include <stdexcept>
+#include <utility>
+
+namespace asic {
+
+class constant_operation final : public abstract_operation {
+public:
+	constant_operation(result_key key, number value)
+		: abstract_operation(std::move(key))
+		, m_value(value) {}
+
+	[[nodiscard]] std::size_t output_count() const noexcept final {
+		return 1;
+	}
+
+private:
+	[[nodiscard]] number evaluate_output_impl(std::size_t, evaluation_context const&) const final {
+		ASIC_DEBUG_MSG("Evaluating constant.");
+		return m_value;
+	}
+
+	number m_value;
+};
+
+class addition_operation final : public binary_operation {
+public:
+	explicit addition_operation(result_key key)
+		: binary_operation(std::move(key)) {}
+
+	[[nodiscard]] std::size_t output_count() const noexcept final {
+		return 1;
+	}
+
+private:
+	[[nodiscard]] number evaluate_output_impl(std::size_t, evaluation_context const& context) const final {
+		ASIC_DEBUG_MSG("Evaluating addition.");
+		return this->evaluate_lhs(context) + this->evaluate_rhs(context);
+	}
+};
+
+class subtraction_operation final : public binary_operation {
+public:
+	explicit subtraction_operation(result_key key)
+		: binary_operation(std::move(key)) {}
+
+	[[nodiscard]] std::size_t output_count() const noexcept final {
+		return 1;
+	}
+
+private:
+	[[nodiscard]] number evaluate_output_impl(std::size_t, evaluation_context const& context) const final {
+		ASIC_DEBUG_MSG("Evaluating subtraction.");
+		return this->evaluate_lhs(context) - this->evaluate_rhs(context);
+	}
+};
+
+class multiplication_operation final : public binary_operation {
+public:
+	explicit multiplication_operation(result_key key)
+		: binary_operation(std::move(key)) {}
+
+	[[nodiscard]] std::size_t output_count() const noexcept final {
+		return 1;
+	}
+
+private:
+	[[nodiscard]] number evaluate_output_impl(std::size_t, evaluation_context const& context) const final {
+		ASIC_DEBUG_MSG("Evaluating multiplication.");
+		return this->evaluate_lhs(context) * this->evaluate_rhs(context);
+	}
+};
+
+class division_operation final : public binary_operation {
+public:
+	explicit division_operation(result_key key)
+		: binary_operation(std::move(key)) {}
+
+	[[nodiscard]] std::size_t output_count() const noexcept final {
+		return 1;
+	}
+
+private:
+	[[nodiscard]] number evaluate_output_impl(std::size_t, evaluation_context const& context) const final {
+		ASIC_DEBUG_MSG("Evaluating division.");
+		return this->evaluate_lhs(context) / this->evaluate_rhs(context);
+	}
+};
+
+class min_operation final : public binary_operation {
+public:
+	explicit min_operation(result_key key)
+		: binary_operation(std::move(key)) {}
+
+	[[nodiscard]] std::size_t output_count() const noexcept final {
+		return 1;
+	}
+
+private:
+	[[nodiscard]] number evaluate_output_impl(std::size_t, evaluation_context const& context) const final {
+		ASIC_DEBUG_MSG("Evaluating min.");
+		auto const lhs = this->evaluate_lhs(context);
+		if (lhs.imag() != 0) {
+			throw std::runtime_error{"Min does not support complex numbers."};
+		}
+		auto const rhs = this->evaluate_rhs(context);
+		if (rhs.imag() != 0) {
+			throw std::runtime_error{"Min does not support complex numbers."};
+		}
+		return std::min(lhs.real(), rhs.real());
+	}
+};
+
+class max_operation final : public binary_operation {
+public:
+	explicit max_operation(result_key key)
+		: binary_operation(std::move(key)) {}
+
+	[[nodiscard]] std::size_t output_count() const noexcept final {
+		return 1;
+	}
+
+private:
+	[[nodiscard]] number evaluate_output_impl(std::size_t, evaluation_context const& context) const final {
+		ASIC_DEBUG_MSG("Evaluating max.");
+		auto const lhs = this->evaluate_lhs(context);
+		if (lhs.imag() != 0) {
+			throw std::runtime_error{"Max does not support complex numbers."};
+		}
+		auto const rhs = this->evaluate_rhs(context);
+		if (rhs.imag() != 0) {
+			throw std::runtime_error{"Max does not support complex numbers."};
+		}
+		return std::max(lhs.real(), rhs.real());
+	}
+};
+
+class square_root_operation final : public unary_operation {
+public:
+	explicit square_root_operation(result_key key)
+		: unary_operation(std::move(key)) {}
+
+	[[nodiscard]] std::size_t output_count() const noexcept final {
+		return 1;
+	}
+
+private:
+	[[nodiscard]] number evaluate_output_impl(std::size_t, evaluation_context const& context) const final {
+		ASIC_DEBUG_MSG("Evaluating sqrt.");
+		return std::sqrt(this->evaluate_input(context));
+	}
+};
+
+class complex_conjugate_operation final : public unary_operation {
+public:
+	explicit complex_conjugate_operation(result_key key)
+		: unary_operation(std::move(key)) {}
+
+	[[nodiscard]] std::size_t output_count() const noexcept final {
+		return 1;
+	}
+
+private:
+	[[nodiscard]] number evaluate_output_impl(std::size_t, evaluation_context const& context) const final {
+		ASIC_DEBUG_MSG("Evaluating conj.");
+		return std::conj(this->evaluate_input(context));
+	}
+};
+
+class absolute_operation final : public unary_operation {
+public:
+	explicit absolute_operation(result_key key)
+		: unary_operation(std::move(key)) {}
+
+	[[nodiscard]] std::size_t output_count() const noexcept final {
+		return 1;
+	}
+
+private:
+	[[nodiscard]] number evaluate_output_impl(std::size_t, evaluation_context const& context) const final {
+		ASIC_DEBUG_MSG("Evaluating abs.");
+		return std::abs(this->evaluate_input(context));
+	}
+};
+
+class constant_multiplication_operation final : public unary_operation {
+public:
+	constant_multiplication_operation(result_key key, number value)
+		: unary_operation(std::move(key))
+		, m_value(value) {}
+
+	[[nodiscard]] std::size_t output_count() const noexcept final {
+		return 1;
+	}
+
+private:
+	[[nodiscard]] number evaluate_output_impl(std::size_t, evaluation_context const& context) const final {
+		ASIC_DEBUG_MSG("Evaluating cmul.");
+		return this->evaluate_input(context) * m_value;
+	}
+
+	number m_value;
+};
+
+class butterfly_operation final : public binary_operation {
+public:
+	explicit butterfly_operation(result_key key)
+		: binary_operation(std::move(key)) {}
+
+	[[nodiscard]] std::size_t output_count() const noexcept final {
+		return 2;
+	}
+
+private:
+	[[nodiscard]] number evaluate_output_impl(std::size_t index, evaluation_context const& context) const final {
+		ASIC_DEBUG_MSG("Evaluating bfly.");
+		if (index == 0) {
+			return this->evaluate_lhs(context) + this->evaluate_rhs(context);
+		}
+		return this->evaluate_lhs(context) - this->evaluate_rhs(context);
+	}
+};
+
+} // namespace asic
+
+#endif // ASIC_SIMULATION_CORE_OPERATIONS_HPP
diff --git a/legacy/simulation_oop/custom_operation.cpp b/legacy/simulation_oop/custom_operation.cpp
index 9153eb5b651f3a481747f3b652f790ce581e70dc..03fbb4b7c8b8790c0004b885736609d9be4d6667 100644
--- a/legacy/simulation_oop/custom_operation.cpp
+++ b/legacy/simulation_oop/custom_operation.cpp
@@ -1,9 +1,8 @@
-#include "custom_operation.h"
+#include "custom_operation.hpp"
 
+#define NOMINMAX
 #include <pybind11/stl.h>
 
-namespace py = pybind11;
-
 namespace asic {
 
 custom_operation::custom_operation(result_key key, pybind11::object evaluate_output, pybind11::object truncate_input,
@@ -27,4 +26,4 @@ number custom_operation::truncate_input(std::size_t index, number value, std::si
 	return m_truncate_input(index, value, bits).cast<number>();
 }
 
-} // namespace asic
\ No newline at end of file
+} // namespace asic
diff --git a/legacy/simulation_oop/custom_operation.hpp b/legacy/simulation_oop/custom_operation.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..12a2a107b7b08c9e6752e61759cb0092937e1eba
--- /dev/null
+++ b/legacy/simulation_oop/custom_operation.hpp
@@ -0,0 +1,36 @@
+#ifndef ASIC_SIMULATION_CUSTOM_OPERATION_HPP
+#define ASIC_SIMULATION_CUSTOM_OPERATION_HPP
+
+#include "../algorithm.hpp"
+#include "../debug.hpp"
+#include "../number.hpp"
+#include "operation.hpp"
+
+#define NOMINMAX
+#include <cstddef>
+#include <fmt/format.h>
+#include <functional>
+#include <pybind11/pybind11.h>
+#include <stdexcept>
+#include <utility>
+
+namespace asic {
+
+class custom_operation final : public nary_operation {
+public:
+	custom_operation(result_key key, pybind11::object evaluate_output, pybind11::object truncate_input, std::size_t output_count);
+
+	[[nodiscard]] std::size_t output_count() const noexcept final;
+
+private:
+	[[nodiscard]] number evaluate_output_impl(std::size_t index, evaluation_context const& context) const final;
+	[[nodiscard]] number truncate_input(std::size_t index, number value, std::size_t bits) const final;
+
+	pybind11::object m_evaluate_output;
+	pybind11::object m_truncate_input;
+	std::size_t m_output_count;
+};
+
+} // namespace asic
+
+#endif // ASIC_SIMULATION_CUSTOM_OPERATION_HPP
diff --git a/legacy/simulation_oop/operation.cpp b/legacy/simulation_oop/operation.cpp
index a9738a0a05287f6ab2d430d4c73560a4c6bd57c5..57a187f8485ceb3f27ee2aceee1fede7c9b4e615 100644
--- a/legacy/simulation_oop/operation.cpp
+++ b/legacy/simulation_oop/operation.cpp
@@ -1,7 +1,8 @@
-#include "operation.h"
+#include "operation.hpp"
 
-#include "../debug.h"
+#include "../debug.hpp"
 
+#define NOMINMAX
 #include <pybind11/pybind11.h>
 
 namespace py = pybind11;
@@ -11,7 +12,7 @@ namespace asic {
 signal_source::signal_source(std::shared_ptr<const operation> op, std::size_t index, std::optional<std::size_t> bits)
 	: m_operation(std::move(op))
 	, m_index(index)
-	, m_bits(std::move(bits)) {}
+	, m_bits(bits) {}
 
 signal_source::operator bool() const noexcept {
 	return static_cast<bool>(m_operation);
@@ -100,7 +101,7 @@ signal_source const& unary_operation::input() const noexcept {
 number unary_operation::evaluate_input(evaluation_context const& context) const {
 	auto const value = m_in.evaluate_output(context);
 	auto const bits = context.bits_override.value_or(m_in.bits().value_or(0));
-	return (context.truncate && bits) ? this->truncate_input(0, value, bits) : value;
+	return (context.truncate && bits != 0) ? this->truncate_input(0, value, bits) : value;
 }
 
 binary_operation::binary_operation(result_key key)
@@ -122,13 +123,13 @@ signal_source const& binary_operation::rhs() const noexcept {
 number binary_operation::evaluate_lhs(evaluation_context const& context) const {
 	auto const value = m_lhs.evaluate_output(context);
 	auto const bits = context.bits_override.value_or(m_lhs.bits().value_or(0));
-	return (context.truncate && bits) ? this->truncate_input(0, value, bits) : value;
+	return (context.truncate && bits != 0) ? this->truncate_input(0, value, bits) : value;
 }
 
 number binary_operation::evaluate_rhs(evaluation_context const& context) const {
 	auto const value = m_rhs.evaluate_output(context);
 	auto const bits = context.bits_override.value_or(m_rhs.bits().value_or(0));
-	return (context.truncate && bits) ? this->truncate_input(0, value, bits) : value;
+	return (context.truncate && bits != 0) ? this->truncate_input(0, value, bits) : value;
 }
 
 nary_operation::nary_operation(result_key key)
@@ -148,9 +149,9 @@ std::vector<number> nary_operation::evaluate_inputs(evaluation_context const& co
 	for (auto const& input : m_inputs) {
 		auto const value = input.evaluate_output(context);
 		auto const bits = context.bits_override.value_or(input.bits().value_or(0));
-		values.push_back((context.truncate && bits) ? this->truncate_input(0, value, bits) : value);
+		values.push_back((context.truncate && bits != 0) ? this->truncate_input(0, value, bits) : value);
 	}
 	return values;
 }
 
-} // namespace asic
\ No newline at end of file
+} // namespace asic
diff --git a/legacy/simulation_oop/operation.hpp b/legacy/simulation_oop/operation.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..086666102a14e32895834c0616af1728e48dc400
--- /dev/null
+++ b/legacy/simulation_oop/operation.hpp
@@ -0,0 +1,134 @@
+#ifndef ASIC_SIMULATION_OPERATION_HPP
+#define ASIC_SIMULATION_OPERATION_HPP
+
+#include "../number.hpp"
+#include "../span.hpp"
+
+#include <cstddef>
+#include <cstdint>
+#include <fmt/format.h>
+#include <memory>
+#include <optional>
+#include <stdexcept>
+#include <string>
+#include <unordered_map>
+#include <utility>
+#include <vector>
+
+namespace asic {
+
+class operation;
+class signal_source;
+
+using result_key = std::string;
+using result_map = std::unordered_map<result_key, std::optional<number>>;
+using delay_map = std::unordered_map<result_key, number>;
+using delay_queue = std::vector<std::pair<result_key, signal_source const*>>;
+
+struct evaluation_context final {
+	result_map* results = nullptr;
+	delay_map* delays = nullptr;
+	delay_queue* deferred_delays = nullptr;
+	std::optional<std::size_t> bits_override{};
+	bool truncate = false;
+};
+
+class signal_source final {
+public:
+	signal_source() noexcept = default;
+	signal_source(std::shared_ptr<const operation> op, std::size_t index, std::optional<std::size_t> bits);
+
+	[[nodiscard]] explicit operator bool() const noexcept;
+
+	[[nodiscard]] std::optional<number> current_output(delay_map const& delays) const;
+	[[nodiscard]] number evaluate_output(evaluation_context const& context) const;
+
+	[[nodiscard]] std::optional<std::size_t> bits() const noexcept;
+
+private:
+	std::shared_ptr<const operation> m_operation{};
+	std::size_t m_index = 0;
+	std::optional<std::size_t> m_bits{};
+};
+
+class operation { // NOLINT(cppcoreguidelines-special-member-functions)
+public:
+	operation() noexcept = default;
+	virtual ~operation() = default;
+
+	[[nodiscard]] virtual std::size_t output_count() const noexcept = 0;
+	[[nodiscard]] virtual std::optional<number> current_output(std::size_t index, delay_map const& delays) const = 0;
+	[[nodiscard]] virtual number evaluate_output(std::size_t index, evaluation_context const& context) const = 0;
+};
+
+class abstract_operation : public operation { // NOLINT(cppcoreguidelines-special-member-functions)
+public:
+	explicit abstract_operation(result_key key);
+	~abstract_operation() override = default;
+
+	[[nodiscard]] std::optional<number> current_output(std::size_t, delay_map const&) const override;
+	[[nodiscard]] number evaluate_output(std::size_t index, evaluation_context const& context) const override;
+
+protected:
+	[[nodiscard]] virtual number evaluate_output_impl(std::size_t index, evaluation_context const& context) const = 0;
+	[[nodiscard]] virtual number truncate_input(std::size_t index, number value, std::size_t bits) const;
+
+	[[nodiscard]] result_key const& key_base() const;
+	[[nodiscard]] result_key key_of_output(std::size_t index) const;
+
+private:
+	result_key m_key;
+};
+
+class unary_operation : public abstract_operation { // NOLINT(cppcoreguidelines-special-member-functions)
+public:
+	explicit unary_operation(result_key key);
+	~unary_operation() override = default;
+
+	void connect(signal_source in);
+
+protected:
+	[[nodiscard]] bool connected() const noexcept;
+	[[nodiscard]] signal_source const& input() const noexcept;
+	[[nodiscard]] number evaluate_input(evaluation_context const& context) const;
+
+private:
+	signal_source m_in;
+};
+
+class binary_operation : public abstract_operation { // NOLINT(cppcoreguidelines-special-member-functions)
+public:
+	explicit binary_operation(result_key key);
+	~binary_operation() override = default;
+
+	void connect(signal_source lhs, signal_source rhs);
+
+protected:
+	[[nodiscard]] signal_source const& lhs() const noexcept;
+	[[nodiscard]] signal_source const& rhs() const noexcept;
+	[[nodiscard]] number evaluate_lhs(evaluation_context const& context) const;
+	[[nodiscard]] number evaluate_rhs(evaluation_context const& context) const;
+
+private:
+	signal_source m_lhs;
+	signal_source m_rhs;
+};
+
+class nary_operation : public abstract_operation { // NOLINT(cppcoreguidelines-special-member-functions)
+public:
+	explicit nary_operation(result_key key);
+	~nary_operation() override = default;
+
+	void connect(std::vector<signal_source> inputs);
+
+protected:
+	[[nodiscard]] span<signal_source const> inputs() const noexcept;
+	[[nodiscard]] std::vector<number> evaluate_inputs(evaluation_context const& context) const;
+
+private:
+	std::vector<signal_source> m_inputs{};
+};
+
+} // namespace asic
+
+#endif // ASIC_SIMULATION_OPERATION_HPP
diff --git a/legacy/simulation_oop/signal_flow_graph.cpp b/legacy/simulation_oop/signal_flow_graph.cpp
index 4c3763c81e8f97f22caf42a47d88c1186b9a874d..5248de741cff3bc7a6e03fc3e479998820c42707 100644
--- a/legacy/simulation_oop/signal_flow_graph.cpp
+++ b/legacy/simulation_oop/signal_flow_graph.cpp
@@ -1,6 +1,6 @@
-#include "signal_flow_graph.h"
+#include "signal_flow_graph.hpp"
 
-#include "../debug.h"
+#include "../debug.hpp"
 
 namespace py = pybind11;
 
@@ -67,8 +67,8 @@ std::shared_ptr<custom_operation> signal_flow_graph_operation::add_custom_operat
 																					std::string_view prefix, result_key key) {
 	auto const input_count = op.attr("input_count").cast<std::size_t>();
 	auto const output_count = op.attr("output_count").cast<std::size_t>();
-	auto const new_op = add_operation<custom_operation>(
-		op, added, key, op.attr("evaluate_output"), op.attr("truncate_input"), output_count);
+	auto new_op = add_operation<custom_operation>(
+		op, added, std::move(key), op.attr("evaluate_output"), op.attr("truncate_input"), output_count);
 	auto inputs = std::vector<signal_source>{};
 	inputs.reserve(input_count);
 	for (auto const i : range(input_count)) {
@@ -141,4 +141,4 @@ std::shared_ptr<operation> signal_flow_graph_operation::make_operation(pybind11:
 	return add_custom_operation(op, added, prefix, std::move(key));
 }
 
-} // namespace asic
\ No newline at end of file
+} // namespace asic
diff --git a/legacy/simulation_oop/signal_flow_graph.hpp b/legacy/simulation_oop/signal_flow_graph.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..c6593af0c9df9a125a477144d08cf546e1b42370
--- /dev/null
+++ b/legacy/simulation_oop/signal_flow_graph.hpp
@@ -0,0 +1,83 @@
+#ifndef ASIC_SIMULATION_SIGNAL_FLOW_GRAPH_HPP
+#define ASIC_SIMULATION_SIGNAL_FLOW_GRAPH_HPP
+
+#include "../algorithm.hpp"
+#include "../debug.hpp"
+#include "../number.hpp"
+#include "core_operations.hpp"
+#include "custom_operation.hpp"
+#include "operation.hpp"
+#include "special_operations.hpp"
+
+#define NOMINMAX
+#include <Python.h>
+#include <cstddef>
+#include <fmt/format.h>
+#include <functional>
+#include <memory>
+#include <pybind11/pybind11.h>
+#include <stdexcept>
+#include <string_view>
+#include <unordered_map>
+#include <utility>
+#include <vector>
+
+namespace asic {
+
+class signal_flow_graph_operation final : public abstract_operation {
+public:
+	using added_operation_cache = std::unordered_map<PyObject const*, std::shared_ptr<operation>>;
+
+	signal_flow_graph_operation(result_key key);
+
+	void create(pybind11::handle sfg, added_operation_cache& added);
+
+	[[nodiscard]] std::vector<std::shared_ptr<input_operation>> const& inputs() noexcept;
+	[[nodiscard]] std::size_t output_count() const noexcept final;
+
+	[[nodiscard]] number evaluate_output(std::size_t index, evaluation_context const& context) const final;
+
+private:
+	[[nodiscard]] number evaluate_output_impl(std::size_t index, evaluation_context const& context) const final;
+
+	[[nodiscard]] static signal_source make_source(pybind11::handle op, std::size_t input_index, added_operation_cache& added,
+												   std::string_view prefix);
+
+	template <typename Operation, typename... Args>
+	[[nodiscard]] static std::shared_ptr<Operation> add_operation(pybind11::handle op, added_operation_cache& added, Args&&... args) {
+		return std::static_pointer_cast<Operation>(
+			added.try_emplace(op.ptr(), std::make_shared<Operation>(std::forward<Args>(args)...)).first->second);
+	}
+
+	template <typename Operation, typename... Args>
+	[[nodiscard]] static std::shared_ptr<Operation> add_unary_operation(pybind11::handle op, added_operation_cache& added,
+																		std::string_view prefix, Args&&... args) {
+		auto new_op = add_operation<Operation>(op, added, std::forward<Args>(args)...);
+		new_op->connect(make_source(op, 0, added, prefix));
+		return new_op;
+	}
+
+	template <typename Operation, typename... Args>
+	[[nodiscard]] static std::shared_ptr<Operation> add_binary_operation(pybind11::handle op, added_operation_cache& added,
+																		 std::string_view prefix, Args&&... args) {
+		auto new_op = add_operation<Operation>(op, added, std::forward<Args>(args)...);
+		new_op->connect(make_source(op, 0, added, prefix), make_source(op, 1, added, prefix));
+		return new_op;
+	}
+
+	[[nodiscard]] static std::shared_ptr<operation> add_signal_flow_graph_operation(pybind11::handle sfg, added_operation_cache& added,
+																					std::string_view prefix, result_key key);
+
+	[[nodiscard]] static std::shared_ptr<custom_operation> add_custom_operation(pybind11::handle op, added_operation_cache& added,
+																				std::string_view prefix, result_key key);
+
+	[[nodiscard]] static std::shared_ptr<operation> make_operation(pybind11::handle op, added_operation_cache& added,
+																   std::string_view prefix);
+
+	std::vector<output_operation> m_output_operations{};
+	std::vector<std::shared_ptr<input_operation>> m_input_operations{};
+};
+
+} // namespace asic
+
+#endif // ASIC_SIMULATION_SIGNAL_FLOW_GRAPH_HPP
diff --git a/legacy/simulation_oop/simulation.cpp b/legacy/simulation_oop/simulation.cpp
index 4d6ff83337a6dade0a0d476e5942ee3b14178195..7d044c6a8a76b5d845f47342b2a82ec9a241a5b4 100644
--- a/legacy/simulation_oop/simulation.cpp
+++ b/legacy/simulation_oop/simulation.cpp
@@ -1,15 +1,13 @@
-#define NOMINMAX
-#include "simulation.h"
+#include "simulation.hpp"
 
-#include "../debug.h"
+#include "../debug.hpp"
 
 namespace py = pybind11;
 
 namespace asic {
 
-simulation::simulation(pybind11::handle sfg, std::optional<std::vector<std::optional<input_provider_t>>> input_providers)
-	: m_sfg("")
-	, m_input_functions(sfg.attr("input_count").cast<std::size_t>(), [](iteration_t) -> number { return number{}; }) {
+simulation::simulation(pybind11::handle sfg, std::optional<std::vector<std::optional<input_provider_type>>> input_providers)
+	: m_input_functions(sfg.attr("input_count").cast<std::size_t>(), [](iteration_type) -> number { return number{}; }) {
 	if (input_providers) {
 		this->set_inputs(std::move(*input_providers));
 	}
@@ -17,29 +15,30 @@ simulation::simulation(pybind11::handle sfg, std::optional<std::vector<std::opti
 	m_sfg.create(sfg, added);
 }
 
-void simulation::set_input(std::size_t index, input_provider_t input_provider) {
+void simulation::set_input(std::size_t index, input_provider_type input_provider) {
 	if (index >= m_input_functions.size()) {
 		throw py::index_error{fmt::format("Input index out of range (expected 0-{}, got {})", m_input_functions.size() - 1, index)};
 	}
-	if (auto* const callable = std::get_if<input_function_t>(&input_provider)) {
+	if (auto* const callable = std::get_if<input_function_type>(&input_provider)) {
 		m_input_functions[index] = std::move(*callable);
 	} else if (auto* const numeric = std::get_if<number>(&input_provider)) {
-		m_input_functions[index] = [value = *numeric](iteration_t) -> number {
+		m_input_functions[index] = [value = *numeric](iteration_type) -> number {
 			return value;
 		};
 	} else if (auto* const list = std::get_if<std::vector<number>>(&input_provider)) {
 		if (!m_input_length) {
-			m_input_length = static_cast<iteration_t>(list->size());
-		} else if (*m_input_length != static_cast<iteration_t>(list->size())) {
+			m_input_length = static_cast<iteration_type>(list->size());
+		} else if (*m_input_length != static_cast<iteration_type>(list->size())) {
 			throw py::value_error{fmt::format("Inconsistent input length for simulation (was {}, got {})", *m_input_length, list->size())};
 		}
-		m_input_functions[index] = [values = std::move(*list)](iteration_t n) -> number {
+		m_input_functions[index] = [values = std::move(*list)](iteration_type n) -> number {
 			return values.at(n);
 		};
 	}
 }
 
-void simulation::set_inputs(std::vector<std::optional<input_provider_t>> input_providers) {
+void simulation::set_inputs(
+	std::vector<std::optional<input_provider_type>> input_providers) { // NOLINT(performance-unnecessary-value-param)
 	if (input_providers.size() != m_input_functions.size()) {
 		throw py::value_error{fmt::format(
 			"Wrong number of inputs supplied to simulation (expected {}, got {})", m_input_functions.size(), input_providers.size())};
@@ -55,7 +54,7 @@ std::vector<number> simulation::step(bool save_results, std::optional<std::size_
 	return this->run_for(1, save_results, bits_override, truncate);
 }
 
-std::vector<number> simulation::run_until(iteration_t iteration, bool save_results, std::optional<std::size_t> bits_override,
+std::vector<number> simulation::run_until(iteration_type iteration, bool save_results, std::optional<std::size_t> bits_override,
 										  bool truncate) {
 	auto result = std::vector<number>{};
 	while (m_iteration < iteration) {
@@ -100,9 +99,9 @@ std::vector<number> simulation::run_until(iteration_t iteration, bool save_resul
 	return result;
 }
 
-std::vector<number> simulation::run_for(iteration_t iterations, bool save_results, std::optional<std::size_t> bits_override,
+std::vector<number> simulation::run_for(iteration_type iterations, bool save_results, std::optional<std::size_t> bits_override,
 										bool truncate) {
-	if (iterations > std::numeric_limits<iteration_t>::max() - m_iteration) {
+	if (iterations > std::numeric_limits<iteration_type>::max() - m_iteration) {
 		throw py::value_error("Simulation iteration type overflow!");
 	}
 	return this->run_until(m_iteration + iterations, save_results, bits_override, truncate);
@@ -115,7 +114,7 @@ std::vector<number> simulation::run(bool save_results, std::optional<std::size_t
 	throw py::index_error{"Tried to run unlimited simulation"};
 }
 
-iteration_t simulation::iteration() const noexcept {
+iteration_type simulation::iteration() const noexcept {
 	return m_iteration;
 }
 
diff --git a/legacy/simulation_oop/simulation.hpp b/legacy/simulation_oop/simulation.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..e10feb990c90ae314576a227d284c5ca2fc4f8ce
--- /dev/null
+++ b/legacy/simulation_oop/simulation.hpp
@@ -0,0 +1,67 @@
+#ifndef ASIC_SIMULATION_OOP_HPP
+#define ASIC_SIMULATION_OOP_HPP
+
+#include "../number.hpp"
+#include "core_operations.hpp"
+#include "custom_operation.hpp"
+#include "operation.hpp"
+#include "signal_flow_graph.hpp"
+#include "special_operations.hpp"
+
+#define NOMINMAX
+#include <cstddef>
+#include <cstdint>
+#include <fmt/format.h>
+#include <functional>
+#include <limits>
+#include <memory>
+#include <optional>
+#include <pybind11/functional.h>
+#include <pybind11/numpy.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include <string_view>
+#include <unordered_map>
+#include <utility>
+#include <variant>
+#include <vector>
+
+namespace asic {
+
+using iteration_type = std::uint32_t;
+using result_array_map = std::unordered_map<std::string, std::vector<number>>;
+using input_function_type = std::function<number(iteration_type)>;
+using input_provider_type = std::variant<number, std::vector<number>, input_function_type>;
+
+class simulation final {
+public:
+	simulation(pybind11::handle sfg, std::optional<std::vector<std::optional<input_provider_type>>> input_providers = std::nullopt);
+
+	void set_input(std::size_t index, input_provider_type input_provider);
+	void set_inputs(std::vector<std::optional<input_provider_type>> input_providers);
+
+	[[nodiscard]] std::vector<number> step(bool save_results, std::optional<std::size_t> bits_override, bool truncate);
+	[[nodiscard]] std::vector<number> run_until(iteration_type iteration, bool save_results, std::optional<std::size_t> bits_override,
+												bool truncate);
+	[[nodiscard]] std::vector<number> run_for(iteration_type iterations, bool save_results, std::optional<std::size_t> bits_override,
+											  bool truncate);
+	[[nodiscard]] std::vector<number> run(bool save_results, std::optional<std::size_t> bits_override, bool truncate);
+
+	[[nodiscard]] iteration_type iteration() const noexcept;
+	[[nodiscard]] pybind11::dict results() const noexcept;
+
+	void clear_results() noexcept;
+	void clear_state() noexcept;
+
+private:
+	signal_flow_graph_operation m_sfg{""};
+	result_array_map m_results{};
+	delay_map m_delays{};
+	iteration_type m_iteration = 0;
+	std::optional<iteration_type> m_input_length{};
+	std::vector<input_function_type> m_input_functions;
+};
+
+} // namespace asic
+
+#endif // ASIC_SIMULATION_OOP_HPP
diff --git a/legacy/simulation_oop/special_operations.cpp b/legacy/simulation_oop/special_operations.cpp
index 1f7a6519a90ba08224585e36093694becf495c4d..0a7100d954049979f59a773124fda6034fd892a3 100644
--- a/legacy/simulation_oop/special_operations.cpp
+++ b/legacy/simulation_oop/special_operations.cpp
@@ -1,6 +1,6 @@
-#include "special_operations.h"
+#include "special_operations.hpp"
 
-#include "../debug.h"
+#include "../debug.hpp"
 
 namespace asic {
 
diff --git a/legacy/simulation_oop/special_operations.hpp b/legacy/simulation_oop/special_operations.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..8c16c01a0e51cd956ee3c556548029c5b17d4be3
--- /dev/null
+++ b/legacy/simulation_oop/special_operations.hpp
@@ -0,0 +1,55 @@
+#ifndef ASIC_SIMULATION_SPECIAL_OPERATIONS_HPP
+#define ASIC_SIMULATION_SPECIAL_OPERATIONS_HPP
+
+#include "../debug.hpp"
+#include "../number.hpp"
+#include "operation.hpp"
+
+#include <cassert>
+#include <cstddef>
+#include <utility>
+
+namespace asic {
+
+class input_operation final : public unary_operation {
+public:
+	explicit input_operation(result_key key);
+
+	[[nodiscard]] std::size_t output_count() const noexcept final;
+	[[nodiscard]] number value() const noexcept;
+	void value(number value) noexcept;
+
+private:
+	[[nodiscard]] number evaluate_output_impl(std::size_t index, evaluation_context const& context) const final;
+
+	number m_value{};
+};
+
+class output_operation final : public unary_operation {
+public:
+	explicit output_operation(result_key key);
+
+	[[nodiscard]] std::size_t output_count() const noexcept final;
+
+private:
+	[[nodiscard]] number evaluate_output_impl(std::size_t index, evaluation_context const& context) const final;
+};
+
+class delay_operation final : public unary_operation {
+public:
+	delay_operation(result_key key, number initial_value);
+
+	[[nodiscard]] std::size_t output_count() const noexcept final;
+
+	[[nodiscard]] std::optional<number> current_output(std::size_t index, delay_map const& delays) const final;
+	[[nodiscard]] number evaluate_output(std::size_t index, evaluation_context const& context) const final;
+
+private:
+	[[nodiscard]] number evaluate_output_impl(std::size_t index, evaluation_context const& context) const final;
+
+	number m_initial_value;
+};
+
+} // namespace asic
+
+#endif // ASIC_SIMULATION_SPECIAL_OPERATIONS_HPP
diff --git a/setup.py b/setup.py
index 41376eb294fdee980533063ac4e207e67750b20b..c412d60e6645f8f224d44863bdf5cb53853d1d48 100644
--- a/setup.py
+++ b/setup.py
@@ -29,7 +29,7 @@ class CMakeBuild(build_ext):
             os.path.dirname(self.get_ext_fullpath(ext.name)))
         cmake_configure_argv = [
             CMAKE_EXE, ext.sourcedir,
-            "-DASIC_BUILDING_PYTHON_DISTRIBUTION=true",
+            "-DASIC_BUILDING_PYTHON_DISTRIBUTION=ON",
             "-DCMAKE_BUILD_TYPE=" + cmake_build_type,
             "-DCMAKE_LIBRARY_OUTPUT_DIRECTORY=" + cmake_output_dir,
             "-DPYTHON_EXECUTABLE=" + sys.executable,
@@ -80,7 +80,7 @@ setuptools.setup(
         "pybind11>=2.3.0",
         "pyside2",
         "qtpy",
-        "graphviz",
+        "graphviz<=0.17",
         "matplotlib",
         "scipy"
     ],
diff --git a/src/algorithm.hpp b/src/algorithm.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..4a609160d3bf9e9e0b065e105dd7ada3f759e8b9
--- /dev/null
+++ b/src/algorithm.hpp
@@ -0,0 +1,350 @@
+#ifndef ASIC_ALGORITHM_HPP
+#define ASIC_ALGORITHM_HPP
+
+#include <cstddef>
+#include <iterator>
+#include <memory>
+#include <type_traits>
+#include <utility>
+
+namespace asic {
+namespace detail {
+
+template <typename Reference>
+class arrow_proxy final {
+public:
+	template <typename Ref>
+	constexpr explicit arrow_proxy(Ref&& r)
+		: m_r(std::forward<Ref>(r)) {}
+
+	Reference* operator->() {
+		return std::addressof(m_r);
+	}
+
+private:
+	Reference m_r;
+};
+
+template <typename T>
+class range_view final {
+public:
+	class iterator final {
+	public:
+		using difference_type = std::ptrdiff_t;
+		using value_type = T const;
+		using reference = value_type&;
+		using pointer = value_type*;
+		using iterator_category = std::random_access_iterator_tag;
+
+		constexpr iterator() noexcept = default;
+		constexpr explicit iterator(T value) noexcept
+			: m_value(value) {}
+
+		[[nodiscard]] constexpr bool operator==(iterator const& other) const noexcept {
+			return m_value == other.m_value;
+		}
+
+		[[nodiscard]] constexpr bool operator!=(iterator const& other) const noexcept {
+			return m_value != other.m_value;
+		}
+
+		[[nodiscard]] constexpr bool operator<(iterator const& other) const noexcept {
+			return m_value < other.m_value;
+		}
+
+		[[nodiscard]] constexpr bool operator>(iterator const& other) const noexcept {
+			return m_value > other.m_value;
+		}
+
+		[[nodiscard]] constexpr bool operator<=(iterator const& other) const noexcept {
+			return m_value <= other.m_value;
+		}
+
+		[[nodiscard]] constexpr bool operator>=(iterator const& other) const noexcept {
+			return m_value >= other.m_value;
+		}
+
+		[[nodiscard]] constexpr reference operator*() const noexcept {
+			return m_value;
+		}
+
+		[[nodiscard]] constexpr pointer operator->() const noexcept {
+			return std::addressof(**this);
+		}
+
+		constexpr iterator& operator++() noexcept {
+			++m_value;
+			return *this;
+		}
+
+		constexpr iterator operator++(int) noexcept {
+			return iterator{m_value++};
+		}
+
+		constexpr iterator& operator--() noexcept {
+			--m_value;
+			return *this;
+		}
+
+		constexpr iterator operator--(int) noexcept {
+			return iterator{m_value--};
+		}
+
+		constexpr iterator& operator+=(difference_type n) noexcept {
+			m_value += n;
+			return *this;
+		}
+
+		constexpr iterator& operator-=(difference_type n) noexcept {
+			m_value -= n;
+			return *this;
+		}
+
+		[[nodiscard]] constexpr T operator[](difference_type n) noexcept {
+			return m_value + static_cast<T>(n);
+		}
+
+		[[nodiscard]] constexpr friend iterator operator+(iterator const& lhs, difference_type rhs) noexcept {
+			return iterator{lhs.m_value + rhs};
+		}
+
+		[[nodiscard]] constexpr friend iterator operator+(difference_type lhs, iterator const& rhs) noexcept {
+			return iterator{lhs + rhs.m_value};
+		}
+
+		[[nodiscard]] constexpr friend iterator operator-(iterator const& lhs, difference_type rhs) noexcept {
+			return iterator{lhs.m_value - rhs};
+		}
+
+		[[nodiscard]] constexpr friend difference_type operator-(iterator const& lhs, iterator const& rhs) noexcept {
+			return static_cast<difference_type>(lhs.m_value - rhs.m_value);
+		}
+
+	private:
+		T m_value{};
+	};
+
+	using sentinel = iterator;
+
+	template <typename First, typename Last>
+	constexpr range_view(First&& first, Last&& last) noexcept
+		: m_begin(std::forward<First>(first))
+		, m_end(std::forward<Last>(last)) {}
+
+	[[nodiscard]] constexpr iterator begin() const noexcept {
+		return m_begin;
+	}
+
+	[[nodiscard]] constexpr sentinel end() const noexcept {
+		return m_end;
+	}
+
+private:
+	iterator m_begin;
+	sentinel m_end;
+};
+
+template <typename Range, typename Iterator, typename Sentinel>
+class enumerate_view final {
+public:
+	using sentinel = Sentinel;
+
+	class iterator final {
+	public:
+		using difference_type = typename std::iterator_traits<Iterator>::difference_type;
+		using value_type = typename std::iterator_traits<Iterator>::value_type;
+		using reference = std::pair<std::size_t const&, decltype(*std::declval<Iterator const>())>;
+		using pointer = arrow_proxy<reference>;
+		using iterator_category =
+			std::common_type_t<typename std::iterator_traits<Iterator>::iterator_category, std::bidirectional_iterator_tag>;
+
+		constexpr iterator() = default;
+
+		constexpr iterator(Iterator it, std::size_t index)
+			: m_it(std::move(it))
+			, m_index(index) {}
+
+		[[nodiscard]] constexpr bool operator==(iterator const& other) const {
+			return m_it == other.m_it;
+		}
+
+		[[nodiscard]] constexpr bool operator!=(iterator const& other) const {
+			return m_it != other.m_it;
+		}
+
+		[[nodiscard]] constexpr bool operator==(sentinel const& other) const {
+			return m_it == other;
+		}
+
+		[[nodiscard]] constexpr bool operator!=(sentinel const& other) const {
+			return m_it != other;
+		}
+
+		[[nodiscard]] constexpr reference operator*() const {
+			return reference{m_index, *m_it};
+		}
+
+		[[nodiscard]] constexpr pointer operator->() const {
+			return pointer{**this};
+		}
+
+		constexpr iterator& operator++() {
+			++m_it;
+			++m_index;
+			return *this;
+		}
+
+		constexpr iterator operator++(int) {
+			return iterator{m_it++, m_index++};
+		}
+
+		constexpr iterator& operator--() {
+			--m_it;
+			--m_index;
+			return *this;
+		}
+
+		constexpr iterator operator--(int) {
+			return iterator{m_it--, m_index--};
+		}
+
+	private:
+		Iterator m_it;
+		std::size_t m_index = 0;
+	};
+
+	template <typename R>
+	constexpr explicit enumerate_view(R&& range)
+		: m_range(std::forward<R>(range)) {}
+
+	[[nodiscard]] constexpr iterator begin() const {
+		return iterator{std::begin(m_range), 0};
+	}
+
+	[[nodiscard]] constexpr sentinel end() const {
+		return std::end(m_range);
+	}
+
+private:
+	Range m_range;
+};
+
+template <typename Range1, typename Range2, typename Iterator1, typename Iterator2, typename Sentinel1, typename Sentinel2>
+class zip_view final {
+public:
+	using sentinel = std::pair<Sentinel1, Sentinel2>;
+
+	class iterator final {
+	public:
+		using difference_type = std::common_type_t<typename std::iterator_traits<Iterator1>::difference_type,
+												   typename std::iterator_traits<Iterator2>::difference_type>;
+		using value_type =
+			std::pair<typename std::iterator_traits<Iterator1>::value_type, typename std::iterator_traits<Iterator2>::value_type>;
+		using reference = std::pair<decltype(*std::declval<Iterator1 const>()), decltype(*std::declval<Iterator2 const>())>;
+		using pointer = arrow_proxy<reference>;
+		using iterator_category =
+			std::common_type_t<typename std::iterator_traits<Iterator1>::iterator_category,
+							   typename std::iterator_traits<Iterator2>::iterator_category, std::bidirectional_iterator_tag>;
+
+		constexpr iterator() = default;
+
+		constexpr iterator(Iterator1 it1, Iterator2 it2)
+			: m_it1(std::move(it1))
+			, m_it2(std::move(it2)) {}
+
+		[[nodiscard]] constexpr bool operator==(iterator const& other) const {
+			return m_it1 == other.m_it1 && m_it2 == other.m_it2;
+		}
+
+		[[nodiscard]] constexpr bool operator!=(iterator const& other) const {
+			return !(*this == other);
+		}
+
+		[[nodiscard]] constexpr bool operator==(sentinel const& other) const {
+			return m_it1 == other.first || m_it2 == other.second;
+		}
+
+		[[nodiscard]] constexpr bool operator!=(sentinel const& other) const {
+			return !(*this == other);
+		}
+
+		[[nodiscard]] constexpr reference operator*() const {
+			return reference{*m_it1, *m_it2};
+		}
+
+		[[nodiscard]] constexpr pointer operator->() const {
+			return pointer{**this};
+		}
+
+		constexpr iterator& operator++() {
+			++m_it1;
+			++m_it2;
+			return *this;
+		}
+
+		constexpr iterator operator++(int) {
+			return iterator{m_it1++, m_it2++};
+		}
+
+		constexpr iterator& operator--() {
+			--m_it1;
+			--m_it2;
+			return *this;
+		}
+
+		constexpr iterator operator--(int) {
+			return iterator{m_it1--, m_it2--};
+		}
+
+	private:
+		Iterator1 m_it1;
+		Iterator2 m_it2;
+	};
+
+	template <typename R1, typename R2>
+	constexpr zip_view(R1&& range1, R2&& range2)
+		: m_range1(std::forward<R1>(range1))
+		, m_range2(std::forward<R2>(range2)) {}
+
+	[[nodiscard]] constexpr iterator begin() const {
+		return iterator{std::begin(m_range1), std::begin(m_range2)};
+	}
+
+	[[nodiscard]] constexpr sentinel end() const {
+		return sentinel{std::end(m_range1), std::end(m_range2)};
+	}
+
+private:
+	Range1 m_range1;
+	Range2 m_range2;
+};
+
+} // namespace detail
+
+template <typename First, typename Last, typename T = std::remove_cv_t<std::remove_reference_t<First>>>
+[[nodiscard]] constexpr auto range(First&& first, Last&& last) {
+	return detail::range_view<T>{std::forward<First>(first), std::forward<Last>(last)};
+}
+
+template <typename Last, typename T = std::remove_cv_t<std::remove_reference_t<Last>>>
+[[nodiscard]] constexpr auto range(Last&& last) {
+	return detail::range_view<T>{T{}, std::forward<Last>(last)};
+}
+
+template <typename Range, typename Iterator = decltype(std::begin(std::declval<Range>())),
+		  typename Sentinel = decltype(std::end(std::declval<Range>()))>
+[[nodiscard]] constexpr auto enumerate(Range&& range) {
+	return detail::enumerate_view<Range, Iterator, Sentinel>{std::forward<Range>(range)};
+}
+
+template <typename Range1, typename Range2, typename Iterator1 = decltype(std::begin(std::declval<Range1>())),
+		  typename Iterator2 = decltype(std::begin(std::declval<Range2>())),
+		  typename Sentinel1 = decltype(std::end(std::declval<Range1>())), typename Sentinel2 = decltype(std::end(std::declval<Range2>()))>
+[[nodiscard]] constexpr auto zip(Range1&& range1, Range2&& range2) {
+	return detail::zip_view<Range1, Range2, Iterator1, Iterator2, Sentinel1, Sentinel2>{std::forward<Range1>(range1),
+																						std::forward<Range2>(range2)};
+}
+
+} // namespace asic
+
+#endif // ASIC_ALGORITHM_HPP
diff --git a/src/debug.hpp b/src/debug.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..e1253790b21aedfb5c9d0205139dd74cf388f8fc
--- /dev/null
+++ b/src/debug.hpp
@@ -0,0 +1,94 @@
+#ifndef ASIC_DEBUG_HPP
+#define ASIC_DEBUG_HPP
+
+#ifndef NDEBUG
+
+#ifndef ASIC_ENABLE_DEBUG_LOGGING
+#define ASIC_ENABLE_DEBUG_LOGGING 1
+#endif
+
+#ifndef ASIC_ENABLE_ASSERTS
+#define ASIC_ENABLE_ASSERTS 1
+#endif
+
+#else
+
+#ifndef ASIC_ENABLE_DEBUG_LOGGING
+#define ASIC_ENABLE_DEBUG_LOGGING 0
+#endif
+
+#ifndef ASIC_ENABLE_ASSERTS
+#define ASIC_ENABLE_ASSERTS 0
+#endif
+
+#endif // NDEBUG
+
+#if ASIC_ENABLE_DEBUG_LOGGING
+#include <filesystem>
+#include <fmt/format.h>
+#include <fstream>
+#include <ostream>
+#include <string_view>
+#include <utility>
+#endif // ASIC_ENABLE_DEBUG_LOGGING
+
+#if ASIC_ENABLE_ASSERTS
+#include <cstdio>
+#include <cstdlib>
+#include <filesystem>
+#include <fmt/format.h>
+#include <string_view>
+#endif // ASIC_ENABLE_ASSERTS
+
+namespace asic {
+
+constexpr auto debug_log_filename = "_b_asic_debug_log.txt";
+
+namespace detail {
+
+#if ASIC_ENABLE_DEBUG_LOGGING
+inline void log_debug_msg_string(std::string_view file, int line, std::string_view string) {
+	static auto log_file = std::ofstream{debug_log_filename, std::ios::trunc};
+	log_file << fmt::format("{:<40}: {}", fmt::format("{}:{}", std::filesystem::path{file}.filename().generic_string(), line), string)
+			 << std::endl;
+}
+
+template <typename Format, typename... Args>
+inline void log_debug_msg(std::string_view file, int line, Format&& format, Args&&... args) {
+	log_debug_msg_string(file, line, fmt::format(std::forward<Format>(format), std::forward<Args>(args)...));
+}
+#endif // ASIC_ENABLE_DEBUG_LOGGING
+
+#if ASIC_ENABLE_ASSERTS
+inline void fail_assert(std::string_view file, int line, std::string_view condition_string) {
+#if ASIC_ENABLE_DEBUG_LOGGING
+	log_debug_msg(file, line, "Assertion failed: {}", condition_string);
+#endif // ASIC_ENABLE_DEBUG_LOGGING
+	fmt::print(stderr, "{}:{}: Assertion failed: {}\n", std::filesystem::path{file}.filename().generic_string(), line, condition_string);
+	std::abort();
+}
+
+template <typename BoolConvertible>
+inline void check_assert(std::string_view file, int line, std::string_view condition_string, BoolConvertible&& condition) {
+	if (!static_cast<bool>(condition)) {
+		fail_assert(file, line, condition_string);
+	}
+}
+#endif // ASIC_ENABLE_ASSERTS
+
+} // namespace detail
+} // namespace asic
+
+#if ASIC_ENABLE_DEBUG_LOGGING
+#define ASIC_DEBUG_MSG(...) (asic::detail::log_debug_msg(__FILE__, __LINE__, __VA_ARGS__))
+#else
+#define ASIC_DEBUG_MSG(...) ((void)0)
+#endif // ASIC_ENABLE_DEBUG_LOGGING
+
+#if ASIC_ENABLE_ASSERTS
+#define ASIC_ASSERT(condition) (asic::detail::check_assert(__FILE__, __LINE__, #condition, (condition)))
+#else
+#define ASIC_ASSERT(condition) ((void)0)
+#endif
+
+#endif // ASIC_DEBUG_HPP
diff --git a/src/main.cpp b/src/main.cpp
index f5c4be532aa47468592e2e2f008308d1724e41b8..8ac36c73b6261d6a7a0daaa9760ed037e9d5cb5a 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -1,9 +1,9 @@
-#include "simulation.h"
-#include <pybind11/pybind11.h>
-
-namespace py = pybind11;
-
-PYBIND11_MODULE(_b_asic, module) {
-	module.doc() = "Better ASIC Toolbox Extension Module.";
-	asic::define_simulation_class(module);
-}
\ No newline at end of file
+#include "simulation.hpp"
+
+#define NOMINMAX
+#include <pybind11/pybind11.h>
+
+PYBIND11_MODULE(_b_asic, module) { // NOLINT
+	module.doc() = "Better ASIC Toolbox Extension Module.";
+	asic::define_simulation_class(module);
+}
diff --git a/src/number.hpp b/src/number.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..869a93f17965108d295022d0ddb2025ce7164ff8
--- /dev/null
+++ b/src/number.hpp
@@ -0,0 +1,14 @@
+#ifndef ASIC_NUMBER_HPP
+#define ASIC_NUMBER_HPP
+
+#define NOMINMAX
+#include <complex>
+#include <pybind11/complex.h>
+
+namespace asic {
+
+using number = std::complex<double>;
+
+} // namespace asic
+
+#endif // ASIC_NUMBER_HPP
diff --git a/src/simulation.cpp b/src/simulation.cpp
index 33280f604be77614f7eadf1ad40d868177dd6e95..2f9ccd3df64dd223f695af62a60a550822525b2a 100644
--- a/src/simulation.cpp
+++ b/src/simulation.cpp
@@ -1,5 +1,6 @@
-#include "simulation.h"
-#include "simulation/simulation.h"
+#include "simulation.hpp"
+
+#include "simulation/simulation.hpp"
 
 namespace py = pybind11;
 
@@ -12,7 +13,7 @@ void define_simulation_class(pybind11::module& module) {
 			py::arg("sfg"),
 			"SFG Constructor.")
 
-		.def(py::init<py::handle, std::optional<std::vector<std::optional<input_provider_t>>>>(),
+		.def(py::init<py::handle, std::optional<std::vector<std::optional<input_provider_type>>>>(),
 			py::arg("sfg"), py::arg("input_providers"),
 			"SFG Constructor.")
 
@@ -58,4 +59,4 @@ void define_simulation_class(pybind11::module& module) {
 	// clang-format on
 }
 
-} // namespace asic
\ No newline at end of file
+} // namespace asic
diff --git a/src/simulation.hpp b/src/simulation.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..5e7ab381e62b35a69aaf3f056e2808182c2e351d
--- /dev/null
+++ b/src/simulation.hpp
@@ -0,0 +1,13 @@
+#ifndef ASIC_SIMULATION_HPP
+#define ASIC_SIMULATION_HPP
+
+#define NOMINMAX
+#include <pybind11/pybind11.h>
+
+namespace asic {
+
+void define_simulation_class(pybind11::module& module);
+
+} // namespace asic
+
+#endif // ASIC_SIMULATION_HPP
diff --git a/src/simulation/compile.cpp b/src/simulation/compile.cpp
index 7fa7ac6721c9af0f5b9decdf8c3f0dca653e47ec..4b0e7fc813efaea3fba291fad0135f753ff9403b 100644
--- a/src/simulation/compile.cpp
+++ b/src/simulation/compile.cpp
@@ -1,11 +1,11 @@
-#define NOMINMAX
-#include "compile.h"
+#include "compile.hpp"
 
-#include "../algorithm.h"
-#include "../debug.h"
-#include "../span.h"
-#include "format_code.h"
+#include "../algorithm.hpp"
+#include "../debug.hpp"
+#include "../span.hpp"
+#include "format_code.hpp"
 
+#define NOMINMAX
 #include <Python.h>
 #include <fmt/format.h>
 #include <limits>
@@ -25,7 +25,7 @@ namespace asic {
 }
 
 [[nodiscard]] static result_key key_of_output(py::handle op, std::size_t output_index, std::string_view prefix) {
-	auto const base = key_base(op, prefix);
+	auto base = key_base(op, prefix);
 	if (base.empty()) {
 		return fmt::to_string(output_index);
 	}
@@ -70,10 +70,10 @@ private:
 	using sfg_info_stack = std::vector<sfg_info>;
 	using delay_queue = std::vector<std::tuple<std::size_t, py::handle, std::string, sfg_info_stack>>;
 	using added_output_cache = std::unordered_set<PyObject const*>;
-	using added_result_cache = std::unordered_map<PyObject const*, result_index_t>;
+	using added_result_cache = std::unordered_map<PyObject const*, result_index_type>;
 	using added_custom_operation_cache = std::unordered_map<PyObject const*, std::size_t>;
 
-	static constexpr auto no_result_index = std::numeric_limits<result_index_t>::max();
+	static constexpr auto no_result_index = std::numeric_limits<result_index_type>::max();
 
 	void initialize_code(std::size_t input_count, std::size_t output_count) {
 		m_code.required_stack_size = 0;
@@ -101,7 +101,7 @@ private:
 	void resolve_invalid_result_indices() {
 		for (auto& instruction : m_code.instructions) {
 			if (instruction.result_index == no_result_index) {
-				instruction.result_index = static_cast<result_index_t>(m_code.result_keys.size());
+				instruction.result_index = static_cast<result_index_type>(m_code.result_keys.size());
 			}
 		}
 	}
@@ -128,7 +128,7 @@ private:
 		return new_sfg_stack;
 	}
 
-	instruction& add_instruction(instruction_type type, result_index_t result_index, std::ptrdiff_t stack_diff) {
+	instruction& add_instruction(instruction_type type, result_index_type result_index, std::ptrdiff_t stack_diff) {
 		m_stack_depth += stack_diff;
 		if (m_stack_depth < 0) {
 			throw py::value_error{"Detected input/output count mismatch in simulation SFG"};
@@ -142,8 +142,9 @@ private:
 		return instruction;
 	}
 
-	[[nodiscard]] std::optional<result_index_t> begin_operation_output(py::handle op, std::size_t output_index, std::string_view prefix) {
-		auto const pointer = op.attr("outputs")[py::int_{output_index}].ptr();
+	[[nodiscard]] std::optional<result_index_type> begin_operation_output(py::handle op, std::size_t output_index,
+																		  std::string_view prefix) {
+		auto* const pointer = op.attr("outputs")[py::int_{output_index}].ptr();
 		if (m_incomplete_outputs.count(pointer) != 0) {
 			// Make sure the output doesn't depend on its own value, unless it's a delay operation.
 			if (op.attr("type_name")().cast<std::string_view>() != "t") {
@@ -151,11 +152,11 @@ private:
 			}
 		}
 		// Try to add a new result.
-		auto const [it, inserted] = m_added_results.try_emplace(pointer, static_cast<result_index_t>(m_code.result_keys.size()));
+		auto const [it, inserted] = m_added_results.try_emplace(pointer, static_cast<result_index_type>(m_code.result_keys.size()));
 		if (inserted) {
-			if (m_code.result_keys.size() >= static_cast<std::size_t>(std::numeric_limits<result_index_t>::max())) {
+			if (m_code.result_keys.size() >= static_cast<std::size_t>(std::numeric_limits<result_index_type>::max())) {
 				throw py::value_error{fmt::format("Simulation SFG requires too many outputs to be stored (limit: {})",
-												  std::numeric_limits<result_index_t>::max())};
+												  std::numeric_limits<result_index_type>::max())};
 			}
 			m_code.result_keys.push_back(key_of_output(op, output_index, prefix));
 			m_incomplete_outputs.insert(pointer);
@@ -168,7 +169,7 @@ private:
 	}
 
 	void end_operation_output(py::handle op, std::size_t output_index) {
-		auto const pointer = op.attr("outputs")[py::int_{output_index}].ptr();
+		auto* const pointer = op.attr("outputs")[py::int_{output_index}].ptr();
 		[[maybe_unused]] auto const erased = m_incomplete_outputs.erase(pointer);
 		ASIC_ASSERT(erased == 1);
 	}
@@ -184,7 +185,7 @@ private:
 		return it->second;
 	}
 
-	[[nodiscard]] std::size_t add_delay_info(number initial_value, result_index_t result_index) {
+	[[nodiscard]] std::size_t add_delay_info(number initial_value, result_index_type result_index) {
 		auto const delay_index = m_code.delays.size();
 		auto& delay = m_code.delays.emplace_back();
 		delay.initial_value = initial_value;
@@ -209,14 +210,14 @@ private:
 		}
 	}
 
-	void add_unary_operation_output(py::handle op, result_index_t result_index, std::string_view prefix, sfg_info_stack const& sfg_stack,
+	void add_unary_operation_output(py::handle op, result_index_type result_index, std::string_view prefix, sfg_info_stack const& sfg_stack,
 									delay_queue& deferred_delays, instruction_type type) {
 		this->add_source(op, 0, prefix, sfg_stack, deferred_delays);
 		this->add_instruction(type, result_index, 0);
 	}
 
-	void add_binary_operation_output(py::handle op, result_index_t result_index, std::string_view prefix, sfg_info_stack const& sfg_stack,
-									 delay_queue& deferred_delays, instruction_type type) {
+	void add_binary_operation_output(py::handle op, result_index_type result_index, std::string_view prefix,
+									 sfg_info_stack const& sfg_stack, delay_queue& deferred_delays, instruction_type type) {
 		this->add_source(op, 0, prefix, sfg_stack, deferred_delays);
 		this->add_source(op, 1, prefix, sfg_stack, deferred_delays);
 		this->add_instruction(type, result_index, -1);
@@ -299,10 +300,10 @@ private:
 		}
 	}
 
-	simulation_code m_code;
-	added_output_cache m_incomplete_outputs;
-	added_result_cache m_added_results;
-	added_custom_operation_cache m_added_custom_operations;
+	simulation_code m_code{};
+	added_output_cache m_incomplete_outputs{};
+	added_result_cache m_added_results{};
+	added_custom_operation_cache m_added_custom_operations{};
 	std::ptrdiff_t m_stack_depth = 0;
 };
 
@@ -310,4 +311,4 @@ simulation_code compile_simulation(pybind11::handle sfg) {
 	return compiler{}.compile(sfg);
 }
 
-} // namespace asic
\ No newline at end of file
+} // namespace asic
diff --git a/src/simulation/compile.hpp b/src/simulation/compile.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..ee7fa3f18c90592428250da25a7259963ddaf2a3
--- /dev/null
+++ b/src/simulation/compile.hpp
@@ -0,0 +1,62 @@
+#ifndef ASIC_SIMULATION_COMPILE_HPP
+#define ASIC_SIMULATION_COMPILE_HPP
+
+#include "instruction.hpp"
+
+#define NOMINMAX
+#include <cstddef>
+#include <pybind11/pybind11.h>
+#include <string>
+#include <vector>
+
+namespace asic {
+
+using result_key = std::string;
+
+struct simulation_code final {
+	struct custom_operation final {
+		// Python function used to evaluate the custom operation.
+		pybind11::object evaluate_output{};
+		// Number of inputs that the custom operation takes.
+		std::size_t input_count = 0;
+		// Number of outputs that the custom operation gives.
+		std::size_t output_count = 0;
+	};
+
+	struct custom_source final {
+		// Index into custom_operations where the custom_operation corresponding to this custom_source is located.
+		std::size_t custom_operation_index = 0;
+		// Output index of the custom_operation that this source gets it value from.
+		std::size_t output_index = 0;
+	};
+
+	struct delay_info final {
+		// Initial value to set at the start of the simulation.
+		number initial_value{};
+		// The result index where the current value should be stored at the start of each iteration.
+		result_index_type result_index = 0;
+	};
+
+	// Instructions to execute for one full iteration of the simulation.
+	std::vector<instruction> instructions{};
+	// Custom operations used by the simulation.
+	std::vector<custom_operation> custom_operations{};
+	// Signal sources that use custom operations.
+	std::vector<custom_source> custom_sources{};
+	// Info about the delay operations used in the simulation.
+	std::vector<delay_info> delays{};
+	// Keys for each result produced by the simulation. The index of the key matches the index of the result in the simulation state.
+	std::vector<result_key> result_keys{};
+	// Number of values expected as input to the simulation.
+	std::size_t input_count = 0;
+	// Number of values given as output from the simulation. This will be the number of values left on the stack after a full iteration of the simulation has been run.
+	std::size_t output_count = 0;
+	// Maximum number of values that need to be able to fit on the stack in order to run a full iteration of the simulation.
+	std::size_t required_stack_size = 0;
+};
+
+[[nodiscard]] simulation_code compile_simulation(pybind11::handle sfg);
+
+} // namespace asic
+
+#endif // ASIC_SIMULATION_COMPILE_HPP
diff --git a/src/simulation/format_code.hpp b/src/simulation/format_code.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..6b286d0e954289065882c8f795d9bdeb4dea3ed2
--- /dev/null
+++ b/src/simulation/format_code.hpp
@@ -0,0 +1,129 @@
+#ifndef ASIC_SIMULATION_FORMAT_CODE_HPP
+#define ASIC_SIMULATION_FORMAT_CODE_HPP
+
+#include "../algorithm.hpp"
+#include "../debug.hpp"
+#include "../number.hpp"
+#include "compile.hpp"
+#include "instruction.hpp"
+
+#include <fmt/format.h>
+#include <string>
+
+namespace asic {
+
+[[nodiscard]] inline std::string format_number(number const& value) {
+	if (value.imag() == 0) {
+		return fmt::to_string(value.real());
+	}
+	if (value.real() == 0) {
+		return fmt::format("{}j", value.imag());
+	}
+	if (value.imag() < 0) {
+		return fmt::format("{}-{}j", value.real(), -value.imag());
+	}
+	return fmt::format("{}+{}j", value.real(), value.imag());
+}
+
+[[nodiscard]] inline std::string format_compiled_simulation_code_result_keys(simulation_code const& code) {
+	auto result = std::string{};
+	for (auto const& [i, result_key] : enumerate(code.result_keys)) {
+		result += fmt::format("{:>2}: \"{}\"\n", i, result_key);
+	}
+	return result;
+}
+
+[[nodiscard]] inline std::string format_compiled_simulation_code_delays(simulation_code const& code) {
+	auto result = std::string{};
+	for (auto const& [i, delay] : enumerate(code.delays)) {
+		ASIC_ASSERT(delay.result_index < code.result_keys.size());
+		result += fmt::format("{:>2}: Initial value: {}, Result: {}: \"{}\"\n",
+							  i,
+							  format_number(delay.initial_value),
+							  delay.result_index,
+							  code.result_keys[delay.result_index]);
+	}
+	return result;
+}
+
+[[nodiscard]] inline std::string format_compiled_simulation_code_instruction(instruction const& instruction) {
+	// clang-format off
+	switch (instruction.type) {
+		case instruction_type::push_input:              return fmt::format("push_input inputs[{}]", instruction.index);
+		case instruction_type::push_result:             return fmt::format("push_result results[{}]", instruction.index);
+		case instruction_type::push_delay:              return fmt::format("push_delay delays[{}]", instruction.index);
+		case instruction_type::push_constant:           return fmt::format("push_constant {}", format_number(instruction.value));
+		case instruction_type::truncate:                return fmt::format("truncate {:#018x}", instruction.bit_mask);
+		case instruction_type::addition:                return "addition";
+		case instruction_type::subtraction:             return "subtraction";
+		case instruction_type::multiplication:          return "multiplication";
+		case instruction_type::division:                return "division";
+		case instruction_type::min:                     return "min";
+		case instruction_type::max:                     return "max";
+		case instruction_type::square_root:             return "square_root";
+		case instruction_type::complex_conjugate:       return "complex_conjugate";
+		case instruction_type::absolute:                return "absolute";
+		case instruction_type::constant_multiplication: return fmt::format("constant_multiplication {}", format_number(instruction.value));
+		case instruction_type::update_delay:            return fmt::format("update_delay delays[{}]", instruction.index);
+		case instruction_type::custom:                  return fmt::format("custom custom_sources[{}]", instruction.index);
+		case instruction_type::forward_value:           return "forward_value";
+	}
+	// clang-format on
+	return std::string{};
+}
+
+[[nodiscard]] inline std::string format_compiled_simulation_code_instructions(simulation_code const& code) {
+	auto result = std::string{};
+	for (auto const& [i, instruction] : enumerate(code.instructions)) {
+		auto instruction_string = format_compiled_simulation_code_instruction(instruction);
+		if (instruction.result_index < code.result_keys.size()) {
+			instruction_string = fmt::format(
+				"{:<26} -> {}: \"{}\"", instruction_string, instruction.result_index, code.result_keys[instruction.result_index]);
+		}
+		result += fmt::format("{:>2}: {}\n", i, instruction_string);
+	}
+	return result;
+}
+
+[[nodiscard]] inline std::string format_compiled_simulation_code(simulation_code const& code) {
+	return fmt::format(
+		"==============================================\n"
+		"> Code stats\n"
+		"==============================================\n"
+		"Input count: {}\n"
+		"Output count: {}\n"
+		"Instruction count: {}\n"
+		"Required stack size: {}\n"
+		"Delay count: {}\n"
+		"Result count: {}\n"
+		"Custom operation count: {}\n"
+		"Custom source count: {}\n"
+		"==============================================\n"
+		"> Delays\n"
+		"==============================================\n"
+		"{}"
+		"==============================================\n"
+		"> Result keys\n"
+		"==============================================\n"
+		"{}"
+		"==============================================\n"
+		"> Instructions\n"
+		"==============================================\n"
+		"{}"
+		"==============================================",
+		code.input_count,
+		code.output_count,
+		code.instructions.size(),
+		code.required_stack_size,
+		code.delays.size(),
+		code.result_keys.size(),
+		code.custom_operations.size(),
+		code.custom_sources.size(),
+		format_compiled_simulation_code_delays(code),
+		format_compiled_simulation_code_result_keys(code),
+		format_compiled_simulation_code_instructions(code));
+}
+
+} // namespace asic
+
+#endif // ASIC_SIMULATION_FORMAT_CODE_HPP
diff --git a/src/simulation/instruction.hpp b/src/simulation/instruction.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..22d0464d1e7cc19bf47cdd2a01fdd807ce53de6e
--- /dev/null
+++ b/src/simulation/instruction.hpp
@@ -0,0 +1,55 @@
+#ifndef ASIC_SIMULATION_INSTRUCTION_HPP
+#define ASIC_SIMULATION_INSTRUCTION_HPP
+
+#include "../number.hpp"
+
+#include <cstddef>
+#include <cstdint>
+#include <optional>
+
+namespace asic {
+
+enum class instruction_type : std::uint8_t {
+	push_input,              // push(inputs[index])
+	push_result,             // push(results[index])
+	push_delay,              // push(delays[index])
+	push_constant,           // push(value)
+	truncate,                // push(trunc(pop(), bit_mask))
+	addition,                // rhs=pop(), lhs=pop(), push(lhs + rhs)
+	subtraction,             // rhs=pop(), lhs=pop(), push(lhs - rhs)
+	multiplication,          // rhs=pop(), lhs=pop(), push(lhs * rhs)
+	division,                // rhs=pop(), lhs=pop(), push(lhs / rhs)
+	min,                     // rhs=pop(), lhs=pop(), push(min(lhs, rhs))
+	max,                     // rhs=pop(), lhs=pop(), push(max(lhs, rhs))
+	square_root,             // push(sqrt(pop()))
+	complex_conjugate,       // push(conj(pop()))
+	absolute,                // push(abs(pop()))
+	constant_multiplication, // push(pop() * value)
+	update_delay,            // delays[index] = pop()
+	custom,                  // Custom operation. Uses custom_source[index].
+	forward_value            // Forward the current value on the stack (push(pop()), i.e. do nothing).
+};
+
+using result_index_type = std::uint16_t;
+
+struct instruction final {
+	constexpr instruction() noexcept // NOLINT(cppcoreguidelines-pro-type-member-init)
+		: index(0) {}
+
+	union {
+		// Index used by push_input, push_result, delay and custom.
+		std::size_t index;
+		// Bit mask used by truncate.
+		std::int64_t bit_mask;
+		// Constant value used by push_constant and constant_multiplication.
+		number value;
+	};
+	// Index into where the result of the instruction will be stored. If the result should be ignored, this index will be one past the last valid result index.
+	result_index_type result_index = 0;
+	// Specifies what kind of operation the instruction should execute.
+	instruction_type type = instruction_type::forward_value;
+};
+
+} // namespace asic
+
+#endif // ASIC_SIMULATION_INSTRUCTION_HPP
diff --git a/src/simulation/run.cpp b/src/simulation/run.cpp
index c14fa192bf0fd52a6b661aa4fb392b78859c382f..da307533b3791e0d2e89cae99a61d107ea5a2301 100644
--- a/src/simulation/run.cpp
+++ b/src/simulation/run.cpp
@@ -1,10 +1,10 @@
-#define NOMINMAX
-#include "run.h"
+#include "run.hpp"
 
-#include "../algorithm.h"
-#include "../debug.h"
-#include "format_code.h"
+#include "../algorithm.hpp"
+#include "../debug.hpp"
+#include "format_code.hpp"
 
+#define NOMINMAX
 #include <algorithm>
 #include <complex>
 #include <cstddef>
@@ -54,7 +54,7 @@ simulation_state run_simulation(simulation_code const& code, span<number const>
 
 	// Setup stack.
 	state.stack.resize(code.required_stack_size);
-	auto stack_pointer = state.stack.data();
+	auto* stack_pointer = state.stack.data();
 
 	// Utility functions to make the stack manipulation code below more readable.
 	// Should hopefully be inlined by the compiler.
@@ -97,21 +97,33 @@ simulation_state run_simulation(simulation_code const& code, span<number const>
 					push(truncate_value(pop(), instruction.bit_mask));
 				}
 				break;
-			case instruction_type::addition:
-				push(pop() + pop());
+			case instruction_type::addition: {
+				auto const rhs = pop();
+				auto const lhs = pop();
+				push(lhs + rhs);
 				break;
-			case instruction_type::subtraction:
-				push(pop() - pop());
+			}
+			case instruction_type::subtraction: {
+				auto const rhs = pop();
+				auto const lhs = pop();
+				push(lhs - rhs);
 				break;
-			case instruction_type::multiplication:
-				push(pop() * pop());
+			}
+			case instruction_type::multiplication: {
+				auto const rhs = pop();
+				auto const lhs = pop();
+				push(lhs * rhs);
 				break;
-			case instruction_type::division:
-				push(pop() / pop());
+			}
+			case instruction_type::division: {
+				auto const rhs = pop();
+				auto const lhs = pop();
+				push(lhs / rhs);
 				break;
+			}
 			case instruction_type::min: {
-				auto const lhs = pop();
 				auto const rhs = pop();
+				auto const lhs = pop();
 				if (lhs.imag() != 0 || rhs.imag() != 0) {
 					throw std::runtime_error{"Min does not support complex numbers."};
 				}
@@ -119,8 +131,8 @@ simulation_state run_simulation(simulation_code const& code, span<number const>
 				break;
 			}
 			case instruction_type::max: {
-				auto const lhs = pop();
 				auto const rhs = pop();
+				auto const lhs = pop();
 				if (lhs.imag() != 0 || rhs.imag() != 0) {
 					throw std::runtime_error{"Max does not support complex numbers."};
 				}
@@ -173,4 +185,4 @@ simulation_state run_simulation(simulation_code const& code, span<number const>
 	return state;
 }
 
-} // namespace asic
\ No newline at end of file
+} // namespace asic
diff --git a/src/simulation/run.hpp b/src/simulation/run.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..32c56b3e764c01080bbf28b93ebb4e1fec4d9556
--- /dev/null
+++ b/src/simulation/run.hpp
@@ -0,0 +1,23 @@
+#ifndef ASIC_SIMULATION_RUN_HPP
+#define ASIC_SIMULATION_RUN_HPP
+
+#include "../number.hpp"
+#include "../span.hpp"
+#include "compile.hpp"
+
+#include <cstdint>
+#include <vector>
+
+namespace asic {
+
+struct simulation_state final {
+	std::vector<number> stack{};
+	std::vector<number> results{};
+};
+
+simulation_state run_simulation(simulation_code const& code, span<number const> inputs, span<number> delays,
+								std::optional<std::uint8_t> bits_override, bool truncate);
+
+} // namespace asic
+
+#endif // ASIC_SIMULATION_RUN_HPP
diff --git a/src/simulation/simulation.cpp b/src/simulation/simulation.cpp
index 3af24c10e62bfe09590a05153abd25468d6bee4c..540505b40d2b09fa301c030ca374102b0d63ba02 100644
--- a/src/simulation/simulation.cpp
+++ b/src/simulation/simulation.cpp
@@ -1,11 +1,11 @@
-#define NOMINMAX
-#include "simulation.h"
+#include "simulation.hpp"
 
-#include "../algorithm.h"
-#include "../debug.h"
-#include "compile.h"
-#include "run.h"
+#include "../algorithm.hpp"
+#include "../debug.hpp"
+#include "compile.hpp"
+#include "run.hpp"
 
+#define NOMINMAX
 #include <fmt/format.h>
 #include <limits>
 #include <pybind11/numpy.h>
@@ -15,9 +15,9 @@ namespace py = pybind11;
 
 namespace asic {
 
-simulation::simulation(pybind11::handle sfg, std::optional<std::vector<std::optional<input_provider_t>>> input_providers)
+simulation::simulation(pybind11::handle sfg, std::optional<std::vector<std::optional<input_provider_type>>> input_providers)
 	: m_code(compile_simulation(sfg))
-	, m_input_functions(sfg.attr("input_count").cast<std::size_t>(), [](iteration_t) -> number { return number{}; }) {
+	, m_input_functions(sfg.attr("input_count").cast<std::size_t>(), [](iteration_type) -> number { return number{}; }) {
 	m_delays.reserve(m_code.delays.size());
 	for (auto const& delay : m_code.delays) {
 		m_delays.push_back(delay.initial_value);
@@ -27,29 +27,30 @@ simulation::simulation(pybind11::handle sfg, std::optional<std::vector<std::opti
 	}
 }
 
-void simulation::set_input(std::size_t index, input_provider_t input_provider) {
+void simulation::set_input(std::size_t index, input_provider_type input_provider) {
 	if (index >= m_input_functions.size()) {
 		throw py::index_error{fmt::format("Input index out of range (expected 0-{}, got {})", m_input_functions.size() - 1, index)};
 	}
-	if (auto* const callable = std::get_if<input_function_t>(&input_provider)) {
+	if (auto* const callable = std::get_if<input_function_type>(&input_provider)) {
 		m_input_functions[index] = std::move(*callable);
 	} else if (auto* const numeric = std::get_if<number>(&input_provider)) {
-		m_input_functions[index] = [value = *numeric](iteration_t) -> number {
+		m_input_functions[index] = [value = *numeric](iteration_type) -> number {
 			return value;
 		};
 	} else if (auto* const list = std::get_if<std::vector<number>>(&input_provider)) {
 		if (!m_input_length) {
-			m_input_length = static_cast<iteration_t>(list->size());
-		} else if (*m_input_length != static_cast<iteration_t>(list->size())) {
+			m_input_length = static_cast<iteration_type>(list->size());
+		} else if (*m_input_length != static_cast<iteration_type>(list->size())) {
 			throw py::value_error{fmt::format("Inconsistent input length for simulation (was {}, got {})", *m_input_length, list->size())};
 		}
-		m_input_functions[index] = [values = std::move(*list)](iteration_t n) -> number {
+		m_input_functions[index] = [values = std::move(*list)](iteration_type n) -> number {
 			return values.at(n);
 		};
 	}
 }
 
-void simulation::set_inputs(std::vector<std::optional<input_provider_t>> input_providers) {
+void simulation::set_inputs(
+	std::vector<std::optional<input_provider_type>> input_providers) { // NOLINT(performance-unnecessary-value-param)
 	if (input_providers.size() != m_input_functions.size()) {
 		throw py::value_error{fmt::format(
 			"Wrong number of inputs supplied to simulation (expected {}, got {})", m_input_functions.size(), input_providers.size())};
@@ -65,7 +66,7 @@ std::vector<number> simulation::step(bool save_results, std::optional<std::uint8
 	return this->run_for(1, save_results, bits_override, truncate);
 }
 
-std::vector<number> simulation::run_until(iteration_t iteration, bool save_results, std::optional<std::uint8_t> bits_override,
+std::vector<number> simulation::run_until(iteration_type iteration, bool save_results, std::optional<std::uint8_t> bits_override,
 										  bool truncate) {
 	auto result = std::vector<number>{};
 	while (m_iteration < iteration) {
@@ -84,9 +85,9 @@ std::vector<number> simulation::run_until(iteration_t iteration, bool save_resul
 	return result;
 }
 
-std::vector<number> simulation::run_for(iteration_t iterations, bool save_results, std::optional<std::uint8_t> bits_override,
+std::vector<number> simulation::run_for(iteration_type iterations, bool save_results, std::optional<std::uint8_t> bits_override,
 										bool truncate) {
-	if (iterations > std::numeric_limits<iteration_t>::max() - m_iteration) {
+	if (iterations > std::numeric_limits<iteration_type>::max() - m_iteration) {
 		throw py::value_error("Simulation iteration type overflow!");
 	}
 	return this->run_until(m_iteration + iterations, save_results, bits_override, truncate);
@@ -99,7 +100,7 @@ std::vector<number> simulation::run(bool save_results, std::optional<std::uint8_
 	throw py::index_error{"Tried to run unlimited simulation"};
 }
 
-iteration_t simulation::iteration() const noexcept {
+iteration_type simulation::iteration() const noexcept {
 	return m_iteration;
 }
 
diff --git a/src/simulation/simulation.hpp b/src/simulation/simulation.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..e5624cd27a801383a07cd64b54ec5c04028ea5eb
--- /dev/null
+++ b/src/simulation/simulation.hpp
@@ -0,0 +1,55 @@
+#ifndef ASIC_SIMULATION_DOD_HPP
+#define ASIC_SIMULATION_DOD_HPP
+
+#include "../number.hpp"
+#include "compile.hpp"
+
+#define NOMINMAX
+#include <cstddef>
+#include <cstdint>
+#include <functional>
+#include <optional>
+#include <pybind11/functional.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include <variant>
+#include <vector>
+
+namespace asic {
+
+using iteration_type = std::uint32_t;
+using input_function_type = std::function<number(iteration_type)>;
+using input_provider_type = std::variant<number, std::vector<number>, input_function_type>;
+
+class simulation final {
+public:
+	simulation(pybind11::handle sfg, std::optional<std::vector<std::optional<input_provider_type>>> input_providers = std::nullopt);
+
+	void set_input(std::size_t index, input_provider_type input_provider);
+	void set_inputs(std::vector<std::optional<input_provider_type>> input_providers);
+
+	[[nodiscard]] std::vector<number> step(bool save_results, std::optional<std::uint8_t> bits_override, bool truncate);
+	[[nodiscard]] std::vector<number> run_until(iteration_type iteration, bool save_results, std::optional<std::uint8_t> bits_override,
+												bool truncate);
+	[[nodiscard]] std::vector<number> run_for(iteration_type iterations, bool save_results, std::optional<std::uint8_t> bits_override,
+											  bool truncate);
+	[[nodiscard]] std::vector<number> run(bool save_results, std::optional<std::uint8_t> bits_override, bool truncate);
+
+	[[nodiscard]] iteration_type iteration() const noexcept;
+	[[nodiscard]] pybind11::dict results() const noexcept;
+
+	void clear_results() noexcept;
+	void clear_state() noexcept;
+
+private:
+	simulation_code m_code;
+	std::vector<input_function_type> m_input_functions;
+	std::vector<number> m_delays{};
+	std::optional<iteration_type> m_input_length{};
+	iteration_type m_iteration = 0;
+	std::vector<std::vector<number>> m_results{};
+};
+
+} // namespace asic
+
+#endif // ASIC_SIMULATION_DOD_HPP
diff --git a/src/span.hpp b/src/span.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..3ccc34eaa9bd485d84808ff812d52e74a98c53f4
--- /dev/null
+++ b/src/span.hpp
@@ -0,0 +1,362 @@
+#ifndef ASIC_SPAN_HPP
+#define ASIC_SPAN_HPP
+
+#include <algorithm>
+#include <array>
+#include <cassert>
+#include <cstddef>
+#include <iterator>
+#include <limits>
+#include <type_traits>
+#include <utility>
+
+namespace asic {
+
+constexpr auto dynamic_size = static_cast<std::size_t>(-1);
+
+// C++17-compatible std::span substitute.
+template <typename T, std::size_t Size = dynamic_size>
+class span;
+
+namespace detail {
+
+template <typename T>
+struct is_span_impl : std::false_type {};
+
+template <typename T, std::size_t Size>
+struct is_span_impl<span<T, Size>> : std::true_type {};
+
+template <typename T>
+struct is_span : is_span_impl<std::remove_cv_t<T>> {};
+
+template <typename T>
+constexpr auto is_span_v = is_span<T>::value;
+
+template <typename T>
+struct is_std_array_impl : std::false_type {};
+
+template <typename T, std::size_t Size>
+struct is_std_array_impl<std::array<T, Size>> : std::true_type {};
+
+template <typename T>
+struct is_std_array : is_std_array_impl<std::remove_cv_t<T>> {};
+
+template <typename T>
+constexpr auto is_std_array_v = is_std_array<T>::value;
+
+template <std::size_t From, std::size_t To>
+struct is_size_convertible : std::bool_constant<From == To || From == dynamic_size || To == dynamic_size> {};
+
+template <std::size_t From, std::size_t To>
+constexpr auto is_size_convertible_v = is_size_convertible<From, To>::value;
+
+template <typename From, typename To>
+struct is_element_type_convertible : std::bool_constant<std::is_convertible_v<From (*)[], To (*)[]>> {};
+
+template <typename From, typename To>
+constexpr auto is_element_type_convertible_v = is_element_type_convertible<From, To>::value;
+
+template <typename T, std::size_t Size>
+struct span_base {
+	using element_type = T;
+	using pointer = element_type*;
+	using size_type = std::size_t;
+
+	constexpr span_base() noexcept = default;
+
+	constexpr span_base(pointer data, [[maybe_unused]] size_type size)
+		: m_data(data) {
+		assert(size == Size);
+	}
+
+	template <size_type N>
+	constexpr span_base(span_base<T, N> other)
+		: m_data(other.data()) {
+		static_assert(N == Size || N == dynamic_size);
+		assert(other.size() == Size);
+	}
+
+	[[nodiscard]] constexpr pointer data() const noexcept {
+		return m_data;
+	}
+
+	[[nodiscard]] constexpr size_type size() const noexcept {
+		return Size;
+	}
+
+private:
+	pointer m_data = nullptr;
+};
+
+template <typename T>
+struct span_base<T, dynamic_size> {
+	using element_type = T;
+	using pointer = element_type*;
+	using size_type = std::size_t;
+
+	constexpr span_base() noexcept = default;
+
+	constexpr span_base(pointer data, size_type size)
+		: m_data(data)
+		, m_size(size) {}
+
+	template <size_type N>
+	explicit constexpr span_base(span_base<T, N> other)
+		: m_data(other.data())
+		, m_size(other.size()) {}
+
+	[[nodiscard]] constexpr pointer data() const noexcept {
+		return m_data;
+	}
+
+	[[nodiscard]] constexpr size_type size() const noexcept {
+		return m_size;
+	}
+
+private:
+	pointer m_data = nullptr;
+	size_type m_size = 0;
+};
+
+template <typename T, std::size_t Size, std::size_t Offset, std::size_t N>
+struct subspan_type {
+	using type = span<T, (N != dynamic_size) ? N : (Size != dynamic_size) ? Size - Offset : Size>;
+};
+
+template <typename T, std::size_t Size, std::size_t Offset, std::size_t Count>
+using subspan_type_t = typename subspan_type<T, Size, Offset, Count>::type;
+
+} // namespace detail
+
+template <typename T, std::size_t Size>
+class span final : public detail::span_base<T, Size> { // NOLINT(cppcoreguidelines-special-member-functions)
+public:
+	using element_type = typename detail::span_base<T, Size>::element_type;
+	using pointer = typename detail::span_base<T, Size>::pointer;
+	using size_type = typename detail::span_base<T, Size>::size_type;
+	using value_type = std::remove_cv_t<element_type>;
+	using reference = element_type&;
+	using iterator = element_type*;
+	using const_iterator = const element_type*;
+	using reverse_iterator = std::reverse_iterator<iterator>;
+	using const_reverse_iterator = std::reverse_iterator<const_iterator>;
+
+	// Default constructor.
+	constexpr span() noexcept = default;
+
+	// Construct from pointer, size.
+	constexpr span(pointer data, size_type size)
+		: detail::span_base<T, Size>(data, size) {}
+
+	// Copy constructor.
+	template <typename U, std::size_t N, typename = std::enable_if_t<detail::is_size_convertible_v<N, Size>>,
+			  typename = std::enable_if_t<detail::is_element_type_convertible_v<U, T>>>
+	constexpr span(span<U, N> const& other)
+		: span(other.data(), other.size()) {}
+
+	// Copy assignment.
+	constexpr span& operator=(span const&) noexcept = default;
+
+	// Destructor.
+	~span() = default;
+
+	// Construct from begin, end.
+	constexpr span(pointer begin, pointer end)
+		: span(begin, end - begin) {}
+
+	// Construct from C array.
+	template <std::size_t N>
+	constexpr span(element_type (&arr)[N]) noexcept
+		: span(std::data(arr), N) {}
+
+	// Construct from std::array.
+	template <std::size_t N, typename = std::enable_if_t<N != 0>>
+	constexpr span(std::array<value_type, N>& arr) noexcept
+		: span(std::data(arr), N) {}
+
+	// Construct from empty std::array.
+	constexpr span(std::array<value_type, 0>&) noexcept
+		: span() {}
+
+	// Construct from const std::array.
+	template <std::size_t N, typename = std::enable_if_t<N != 0>>
+	constexpr span(std::array<value_type, N> const& arr) noexcept
+		: span(std::data(arr), N) {}
+
+	// Construct from empty const std::array.
+	constexpr span(std::array<value_type, 0> const&) noexcept
+		: span() {}
+
+	// Construct from other container.
+	template <
+		typename Container, typename = std::enable_if_t<!detail::is_span_v<Container>>,
+		typename = std::enable_if_t<!detail::is_std_array_v<Container>>, typename = decltype(std::data(std::declval<Container>())),
+		typename = decltype(std::size(std::declval<Container>())),
+		typename = std::enable_if_t<std::is_convertible_v<typename Container::pointer, pointer>>,
+		typename = std::enable_if_t<std::is_convertible_v<typename Container::pointer, decltype(std::data(std::declval<Container>()))>>>
+	constexpr span(Container& container)
+		: span(std::data(container), std::size(container)) {}
+
+	// Construct from other const container.
+	template <
+		typename Container, typename Element = element_type, typename = std::enable_if_t<std::is_const_v<Element>>,
+		typename = std::enable_if_t<!detail::is_span_v<Container>>, typename = std::enable_if_t<!detail::is_std_array_v<Container>>,
+		typename = decltype(std::data(std::declval<Container>())), typename = decltype(std::size(std::declval<Container>())),
+		typename = std::enable_if_t<std::is_convertible_v<typename Container::pointer, pointer>>,
+		typename = std::enable_if_t<std::is_convertible_v<typename Container::pointer, decltype(std::data(std::declval<Container>()))>>>
+	constexpr span(Container const& container)
+		: span(std::data(container), std::size(container)) {}
+
+	[[nodiscard]] constexpr iterator begin() const noexcept {
+		return this->data();
+	}
+
+	[[nodiscard]] constexpr const_iterator cbegin() const noexcept {
+		return this->data();
+	}
+
+	[[nodiscard]] constexpr iterator end() const noexcept {
+		return this->data() + this->size();
+	}
+
+	[[nodiscard]] constexpr const_iterator cend() const noexcept {
+		return this->data() + this->size();
+	}
+
+	[[nodiscard]] constexpr reverse_iterator rbegin() const noexcept {
+		return std::make_reverse_iterator(this->end());
+	}
+
+	[[nodiscard]] constexpr const_reverse_iterator crbegin() const noexcept {
+		return std::make_reverse_iterator(this->cend());
+	}
+
+	[[nodiscard]] constexpr reverse_iterator rend() const noexcept {
+		return std::make_reverse_iterator(this->begin());
+	}
+
+	[[nodiscard]] constexpr const_reverse_iterator crend() const noexcept {
+		return std::make_reverse_iterator(this->cbegin());
+	}
+
+	[[nodiscard]] constexpr reference operator[](size_type i) const noexcept {
+		assert(i < this->size());
+		return this->data()[i];
+	}
+
+	[[nodiscard]] constexpr reference operator()(size_type i) const noexcept {
+		assert(i < this->size());
+		return this->data()[i];
+	}
+
+	[[nodiscard]] constexpr size_type size_bytes() const noexcept {
+		return this->size() * sizeof(element_type);
+	}
+
+	[[nodiscard]] constexpr bool empty() const noexcept {
+		return this->size() == 0;
+	}
+
+	[[nodiscard]] constexpr reference front() const noexcept {
+		assert(!this->empty());
+		return this->data()[0];
+	}
+
+	[[nodiscard]] constexpr reference back() const noexcept {
+		assert(!this->empty());
+		return this->data()[this->size() - 1];
+	}
+
+	template <std::size_t N>
+	[[nodiscard]] constexpr span<T, N> first() const {
+		static_assert(N != dynamic_size && N <= Size);
+		return {this->data(), N};
+	}
+
+	template <std::size_t N>
+	[[nodiscard]] constexpr span<T, N> last() const {
+		static_assert(N != dynamic_size && N <= Size);
+		return {this->data() + (Size - N), N};
+	}
+
+	template <std::size_t Offset, std::size_t N = dynamic_size>
+	[[nodiscard]] constexpr auto subspan() const -> detail::subspan_type_t<T, Size, Offset, N> {
+		static_assert(Offset <= Size);
+		return {this->data() + Offset, (N == dynamic_size) ? this->size() - Offset : N};
+	}
+
+	[[nodiscard]] constexpr span<T, dynamic_size> first(size_type n) const {
+		assert(n <= this->size());
+		return {this->data(), n};
+	}
+
+	[[nodiscard]] constexpr span<T, dynamic_size> last(size_type n) const {
+		return this->subspan(this->size() - n);
+	}
+
+	[[nodiscard]] constexpr span<T, dynamic_size> subspan(size_type offset, size_type n = dynamic_size) const {
+		if constexpr (Size == dynamic_size) {
+			assert(offset <= this->size());
+			if (n == dynamic_size) {
+				return {this->data() + offset, this->size() - offset};
+			}
+			assert(n <= this->size());
+			assert(offset + n <= this->size());
+			return {this->data() + offset, n};
+		} else {
+			return span<T, dynamic_size>{*this}.subspan(offset, n);
+		}
+	}
+};
+
+template <typename T, std::size_t LhsSize, std::size_t RhsSize>
+[[nodiscard]] constexpr bool operator==(span<T, LhsSize> lhs, span<T, RhsSize> rhs) {
+	return std::equal(lhs.begin(), lhs.end(), rhs.begin(), rhs.end());
+}
+
+template <typename T, std::size_t LhsSize, std::size_t RhsSize>
+[[nodiscard]] constexpr bool operator!=(span<T, LhsSize> lhs, span<T, RhsSize> rhs) {
+	return !(lhs == rhs);
+}
+
+template <typename T, std::size_t LhsSize, std::size_t RhsSize>
+[[nodiscard]] constexpr bool operator<(span<T, LhsSize> lhs, span<T, RhsSize> rhs) {
+	return std::lexicographical_compare(lhs.begin(), lhs.end(), rhs.begin(), rhs.end());
+}
+
+template <typename T, std::size_t LhsSize, std::size_t RhsSize>
+[[nodiscard]] constexpr bool operator<=(span<T, LhsSize> lhs, span<T, RhsSize> rhs) {
+	return !(lhs > rhs);
+}
+
+template <typename T, std::size_t LhsSize, std::size_t RhsSize>
+[[nodiscard]] constexpr bool operator>(span<T, LhsSize> lhs, span<T, RhsSize> rhs) {
+	return rhs < lhs;
+}
+
+template <typename T, std::size_t LhsSize, std::size_t RhsSize>
+[[nodiscard]] constexpr bool operator>=(span<T, LhsSize> lhs, span<T, RhsSize> rhs) {
+	return !(lhs < rhs);
+}
+
+template <typename Container>
+span(Container&) -> span<typename Container::value_type>;
+
+template <typename Container>
+span(Container const&) -> span<typename Container::value_type const>;
+
+template <typename T, std::size_t N>
+span(T (&)[N]) -> span<T, N>;
+
+template <typename T, std::size_t N>
+span(std::array<T, N>&) -> span<T, N>;
+
+template <typename T, std::size_t N>
+span(std::array<T, N> const&) -> span<T const, N>;
+
+template <typename T, typename Dummy>
+span(T, Dummy &&) -> span<std::remove_reference_t<decltype(std::declval<T>()[0])>>;
+
+} // namespace asic
+
+#endif // ASIC_SPAN_HPP