Initial (with old code)HEADmaster

This includes old code too, perhaps from 2021-05-05, judging from the
birth timestamp on the directory. The old code included option parsing
and equation solving for fancy mutually-dependent options, but no actual
rendering (imagine that).

12 files changed, 1060 insertions, 0 deletions

diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..f211247
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,6 @@
+/.obj/
+/buddhabrot
+/compile_commands.json
+/.ccls-cache/
+
+*.png
diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..ed2da60
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,35 @@
+CXX = g++
+CXXFLAGS = -Wall -Wextra -std=c++20 -march=native
+CXXFLAGS += -D_DEFAULT_SOURCE
+LDFLAGS = -lpng
+BIN = buddhabrot
+
+OBJDIR = .obj
+
+ifneq ($(DEBUG),)
+	CXXFLAGS += -Og -g -fsanitize=address
+	LDFLAGS += -fsanitize=address
+else
+	CXXFLAGS += -O3
+endif
+
+SOURCES := $(wildcard *.cpp)
+OBJECTS := $(patsubst %.cpp,$(OBJDIR)/%.o,$(SOURCES))
+HEADERS := $(wildcard *.h)
+
+.PHONY: all clean
+
+all: $(BIN)
+
+clean:
+	rm -rf $(BIN) $(OBJDIR)
+
+
+$(BIN): $(OBJECTS)
+	$(CXX) -o $@ $^ $(LDFLAGS)
+
+$(OBJDIR)/%.o: %.cpp $(HEADERS) | $(OBJDIR)
+	$(CXX) $(CXXFLAGS) -c -o $@ $<
+
+$(OBJDIR):
+	mkdir $@
diff --git a/bench.sh b/bench.sh
new file mode 100755
index 0000000..90a4cc5
--- /dev/null
+++ b/bench.sh
@@ -0,0 +1,2 @@
+#!/usr/bin/env bash
+time ./buddhabrot -L -2.3 -R 1.3 -Y 0 -w 2560 -h 1440 -n 6e8 -m 40 -t 5 "$@" -o /dev/null
diff --git a/compute.cpp b/compute.cpp
new file mode 100644
index 0000000..d3f5e42
--- /dev/null
+++ b/compute.cpp
@@ -0,0 +1,191 @@
+#include "compute.h"
+#include <complex>
+#include <random>
+#include <thread>
+#include <memory>
+#include <cstdlib>
+#include <cassert>
+#include <immintrin.h>
+
+struct Pointgen {
+	Pointgen(std::mt19937::result_type seed) : randgen{seed} {}
+	double operator()() {
+		return dist(randgen);
+	}
+
+private:
+	std::mt19937 randgen;
+	std::uniform_real_distribution<double> dist{-1.5, 1.5};
+};
+
+std::ostream& operator<<(std::ostream &os, __m256d x) {
+	union {
+		__m256d vec;
+		double arr[4];
+	} u;
+	u.vec = x;
+	return os << '<' << u.arr[0] << ',' << u.arr[1] << ',' << u.arr[2] << ',' << u.arr[3] << '>';
+}
+
+std::ostream& operator<<(std::ostream &os, __m256i x) {
+	union {
+		__m256i vec;
+		int64_t arr[4];
+	} u;
+	u.vec = x;
+	return os << '<' << u.arr[0] << ',' << u.arr[1] << ',' << u.arr[2] << ',' << u.arr[3] << '>';
+}
+
+union int64_4vec {
+	__m256i i;
+	int64_t a[4];
+};
+
+static HistImage accumulate_brot_simple(const Options &options, int64_t numorbits_override, std::mt19937::result_type seed) {
+	HistImage img{(size_t)options.imgwidth.value, (size_t)options.imgheight.value};
+
+	Pointgen gen{seed};
+
+	// trace[0] is never set; considered equal to c
+	std::vector<std::pair<double, double>> trace(options.maxiter.value + 1);
+
+	for (int64_t orbiti = 0; orbiti < numorbits_override; orbiti++) {
+		const double cx = gen(), cy = gen();
+		// std::cout << "cx=" << cx << " cy=" << cy << std::endl;
+
+		int iteri = 0;
+		{
+			double a = cx, b = cy;
+			if (a * a + b * b >= 4) continue;
+			while (iteri < options.maxiter.value) {
+				const double newb = 2 * a * b + cy;
+				a = a * a - b * b + cx;
+				b = newb;
+				iteri++;
+				trace[iteri] = std::make_pair(a, b);
+
+				// std::cout << "a=" << a << " b=" << b << std::endl;
+				if (a * a + b * b >= 4) goto escape;
+			}
+			continue;  // did not escape, so in set, so ignore
+		}
+
+	escape:
+		if (options.orbit_burnin.value == 0) trace[0] = std::make_pair(cx, cy);
+		// std::cout << "iteri=" << iteri << std::endl;
+		for (int i = options.orbit_burnin.value; i <= iteri; i++) {
+			const double zx = trace[i].first, zy = trace[i].second;
+			const int imgx = (zx - options.xmin.value) / options.cplxwidth.value * (options.imgwidth.value - 1);
+			const int imgy = (zy - options.ymin.value) / options.cplxheight.value * (options.imgheight.value - 1);
+			// std::cerr << " zx=" << zx << " zy=" << zy << " img=(" << imgx << ',' << imgy << ")" << std::endl;
+			if (0 <= imgx && imgx < options.imgwidth.value &&
+					0 <= imgy && imgy < options.imgheight.value) {
+				img.data[options.imgwidth.value * imgy + imgx]++;
+			}
+		}
+	}
+
+	return img;
+}
+
+static HistImage accumulate_brot_vectorised(const Options &options, int64_t numorbits_override, std::mt19937::result_type seed) {
+	HistImage img{(size_t)options.imgwidth.value, (size_t)options.imgheight.value};
+
+	Pointgen gen{seed};
+
+	using F4 = __m256d;
+	using I4 = __m256i;
+
+	// trace[0] is never set; considered equal to c
+	std::vector<std::pair<F4, F4>> trace(options.maxiter.value + 1);
+
+	for (int64_t orbiti = 0; orbiti < numorbits_override; orbiti += 4) {
+		const F4 cx = _mm256_set_pd(gen(), gen(), gen(), gen());
+		const F4 cy = _mm256_set_pd(gen(), gen(), gen(), gen());
+		// std::cout << "cx=" << cx << " cy=" << cy << std::endl;
+
+		I4 escaped_mask = _mm256_setzero_si256();
+		I4 escaped_at = _mm256_setzero_si256();
+
+		int iteri = 0;
+		{
+			F4 a = cx, b = cy;
+			if (_mm256_movemask_pd(a * a + b * b < 4) == 0) continue;
+			while (iteri < options.maxiter.value) {
+				const F4 newb = 2 * a * b + cy;
+				a = a * a - b * b + cx;
+				b = newb;
+				iteri++;
+				trace[iteri] = std::make_pair(a, b);
+
+				const I4 remaining = a * a + b * b < 4;
+				const I4 escaped = ~remaining;
+				escaped_at |= _mm256_set1_epi64x(iteri) & ~escaped_mask & escaped;
+				escaped_mask |= escaped;
+				// std::cout << "a=" << a << " b=" << b << " rem=" << remaining << " esc=" << escaped << " at=" << escaped_at << " mask=" << escaped_mask << std::endl;
+				if (_mm256_movemask_epi8(remaining) == 0) goto escape;
+				// a = _mm256_castsi256_pd(_mm256_castpd_si256(a) & remaining);
+				// b = _mm256_castsi256_pd(_mm256_castpd_si256(b) & remaining);
+			}
+		}
+		if (_mm256_movemask_epi8(escaped_mask) == 0) continue;  // if none escaped, nothing to do
+
+	escape:
+		if (options.orbit_burnin.value == 0) trace[0] = std::make_pair(cx, cy);
+		union int64_4vec escaped_at_un;
+		_mm256_store_si256(&escaped_at_un.i, escaped_at);
+		for (int tracei = 0; tracei < 4; tracei++) {
+			const int limit = escaped_at_un.a[tracei];
+			// std::cout << "tracei=" << tracei << " limit=" << limit << std::endl;
+			for (int i = options.orbit_burnin.value; i <= limit; i++) {
+				const double zx = ((const double*)&trace[i].first)[tracei];
+				const double zy = ((const double*)&trace[i].second)[tracei];
+				const int imgx = (zx - options.xmin.value) / options.cplxwidth.value * (options.imgwidth.value - 1);
+				const int imgy = (zy - options.ymin.value) / options.cplxheight.value * (options.imgheight.value - 1);
+				// std::cerr << " zx=" << zx << " zy=" << zy << " img=(" << imgx << ',' << imgy << ")" << std::endl;
+				if (0 <= imgx && imgx < options.imgwidth.value &&
+						0 <= imgy && imgy < options.imgheight.value) {
+					img.data[options.imgwidth.value * imgy + imgx]++;
+				}
+			}
+		}
+	}
+
+	return img;
+}
+
+static HistImage accumulate_brot(const Options &options, int64_t numorbits_override, std::mt19937::result_type seed) {
+	if (options.algorithm.value == "simple") return accumulate_brot_simple(options, numorbits_override, seed);
+	if (options.algorithm.value == "vectorised") return accumulate_brot_vectorised(options, numorbits_override, seed);
+	std::cerr << "Invalid algorithm '" << options.algorithm.value << "'" << std::endl;
+	exit(1);
+}
+
+HistImage compute_brot(const Options &options) {
+	std::vector<std::unique_ptr<HistImage>> images(options.num_threads.value);
+
+	std::mt19937 randgen(std::random_device{}());
+
+	std::vector<std::thread> ths(options.num_threads.value);
+	for (size_t i = 0; i < (size_t)options.num_threads.value; i++) {
+		std::mt19937::result_type seed = randgen();
+		ths[i] = std::thread{[&options, &images, i, seed]() {
+			images[i] = std::make_unique<HistImage>(
+				accumulate_brot(options, options.numorbits.value / options.num_threads.value, seed));
+		}};
+	}
+
+	for (std::thread &th : ths) th.join();
+
+	HistImage first = std::move(*images[0]);
+	images[0].reset();
+
+	for (int i = 1; i < options.num_threads.value; i++) {
+		first.add(*images[i]);
+		images[i].reset();
+	}
+
+	return first;
+}
+
+// vim: set sw=4 ts=4 noet:
diff --git a/compute.h b/compute.h
new file mode 100644
index 0000000..5cd8915
--- /dev/null
+++ b/compute.h
@@ -0,0 +1,9 @@
+#pragma once
+
+#include "image.h"
+#include "options.h"
+
+
+HistImage compute_brot(const Options &options);
+
+// vim: set sw=4 ts=4 noet:
diff --git a/eqsystem_solve.cpp b/eqsystem_solve.cpp
new file mode 100644
index 0000000..878ddb6
--- /dev/null
+++ b/eqsystem_solve.cpp
@@ -0,0 +1,330 @@
+#include <iostream>
+#include <sstream>
+#include <unordered_set>
+#include <stdexcept>
+#include <algorithm>
+#include <cassert>
+#include "eqsystem_solve.h"
+
+
+const Product Product::times(double x) const {
+	return Product{cnst * x, vars};
+}
+
+void Product::times_inplace(double x) {
+	cnst *= x;
+}
+
+const Product Product::times(const Product &p) const {
+	std::vector<std::string> vars2{vars};
+	vars2.insert(vars2.end(), p.vars.begin(), p.vars.end());
+	return Product{cnst * p.cnst, vars2};
+}
+
+bool Product::contains(const std::string &var) const {
+	for (const std::string &v : vars) {
+		if (v == var) return true;
+	}
+	return false;
+}
+
+const Product Product::without(const std::string &var) const {
+	Product result{cnst, {}};
+	result.vars.reserve(vars.size());
+	for (const std::string &v : vars) {
+		if (v != var) result.vars.push_back(v);
+	}
+	return result;
+}
+
+bool Product::operator==(const Product &other) const {
+	if (cnst != other.cnst) return false;
+	if (vars.size() != other.vars.size()) return false;
+	std::unordered_set<std::string> seen(vars.begin(), vars.end());
+	for (const std::string &v : other.vars) {
+		auto it = seen.find(v);
+		if (it == seen.end()) return false;
+		seen.erase(it);
+	}
+	return true;
+}
+
+const Sum Sum::times(double x) const {
+	Sum result;
+	result.terms.reserve(terms.size());
+	for (const Product &p : terms) result.terms.push_back(p.times(x));
+	return result;
+}
+
+const Sum Sum::substitute(const std::string &var, Sum sum) const {
+	Sum result;
+
+	for (const Product &p : terms) {
+		if (p.contains(var)) {
+			Product rest = p.without(var);
+			for (const Product &term : sum.terms) {
+				result.terms.push_back(rest.times(term));
+			}
+		} else {
+			result.terms.push_back(p);
+		}
+	}
+
+	return result;
+}
+
+void Sum::simplify() {
+	double tot = 0.0;
+	std::vector<Product> terms2;
+
+	for (Product &p : terms) {
+		if (p.vars.empty()) {
+			tot += p.cnst;
+		} else {
+			sort(p.vars.begin(), p.vars.end());
+			bool added = false;
+			for (Product &q : terms2) {
+				if (q.vars == p.vars) {
+					q.cnst += p.cnst;
+					added = true;
+					break;
+				}
+			}
+			if (!added) terms2.push_back(std::move(p));
+		}
+	}
+
+	if (tot != 0.0) terms2.push_back(Product{tot, {}});
+	terms = terms2;
+}
+
+std::optional<double> Sum::evaluate() const {
+	double tot = 0.0;
+	for (const Product &p : terms) {
+		if (!p.vars.empty()) return std::nullopt;
+		tot += p.cnst;
+	}
+	return tot;
+}
+
+std::optional<std::string> Equation::validate() const {
+	std::unordered_set<std::string> seen;
+
+	for (int i = 0; i < 2; i++) {
+		const Sum &sum = i == 0 ? lhs : rhs;
+
+		for (const Product &p : sum.terms) {
+			for (const std::string &var : p.vars) {
+				auto it = seen.find(var);
+				if (it != seen.end()) {
+					std::stringstream ss;
+					ss << "In equation " << *this << ": variable " << var << " occurs multiple times";
+					return ss.str();
+				}
+				seen.insert(var);
+			}
+		}
+	}
+
+	return std::nullopt;
+}
+
+const Equation Equation::substitute(const std::string &var, Sum sum) const {
+	return Equation{lhs.substitute(var, sum), rhs.substitute(var, sum)};
+}
+
+void Equation::substitute_inplace(const std::string &var, Sum sum) {
+	*this = substitute(var, sum);
+}
+
+int Equation::num_vars() const {
+	int count = 0;
+	for (const Product &p : lhs.terms) count += p.vars.size();
+	for (const Product &p : rhs.terms) count += p.vars.size();
+	return count;
+}
+
+std::vector<std::string> Equation::all_vars() const {
+	std::vector<std::string> result;
+	for (const Product &p : lhs.terms) result.insert(result.end(), p.vars.begin(), p.vars.end());
+	for (const Product &p : rhs.terms) result.insert(result.end(), p.vars.begin(), p.vars.end());
+	return result;
+}
+
+bool Equation::isolatable(const std::string &target) const {
+	int num_terms = 0;
+
+	for (int i = 0; i < 2; i++) {
+		const Sum &side = i == 0 ? lhs : rhs;
+		for (const Product &p : side.terms) {
+			if (p.contains(target)) {
+				num_terms++;
+				if (num_terms > 1) return false;
+				if (p.vars.size() > 1) return false;
+			}
+		}
+	}
+
+	return true;
+}
+
+void Equation::isolate_left(const std::string &target) {
+	Sum newlhs, newrhs;
+	for (int i = 0; i < 2; i++) {
+		const Sum &side = i == 0 ? lhs : rhs;
+		for (const Product &p : side.terms) {
+			if (p.contains(target)) newlhs.terms.push_back(p.times(i == 0 ? 1 : -1));
+			else newrhs.terms.push_back(p.times(i == 1 ? 1 : -1));
+		}
+	}
+
+	lhs = newlhs;
+	rhs = newrhs;
+
+	if (lhs.terms.size() == 1) {
+		const double factor = lhs.terms.at(0).cnst;
+		if (factor != 0) {
+			for (Product &p : rhs.terms) p.times_inplace(1.0 / factor);
+		}
+		lhs.terms.at(0).cnst = 1.0;
+	}
+}
+
+void Equation::simplify() {
+	lhs.simplify();
+	rhs.simplify();
+}
+
+// xmax = xmin + cplxwidth
+// ymax = ymin + cplxheight
+// xmin + xmax = 2 * cx
+// ymin + ymax = 2 * cy
+// imgwidth * cplxheight = cplxwidth * imgheight
+
+std::variant<
+	std::string,  // description of the error
+	std::vector<std::pair<std::string, double>>  // result assignment
+> System::solve_inplace(bool debug) {
+	for (const Equation &eq : eqs) {
+		if (auto err = eq.validate()) {
+			throw std::invalid_argument("Invalid equation passed to System::solve(): " + *err);
+		}
+	}
+
+	std::vector<std::pair<std::string, double>> assignment;
+
+	std::vector<bool> assigned(eqs.size());
+	std::vector<bool> substituted(eqs.size());
+
+success_try_again:
+	if (debug) {
+		std::cerr << std::endl << "\x1B[1mSolving:\x1B[0m" << std::endl;
+		std::cerr << *this << std::endl << std::endl;
+	}
+
+	for (int stage = 1; stage <= 2; stage++) {
+		for (size_t i = 0; i < eqs.size(); i++) {
+			if (assigned[i] || (stage == 2 && substituted[i])) continue;
+
+			Equation &eq = eqs[i];
+
+			if (stage >= 2 || eq.num_vars() == 1) {
+				const std::string var = eq.all_vars().at(0);
+				if (eq.isolatable(var)) {
+					eq.isolate_left(var);
+					eq.rhs.simplify();
+					if (eq.lhs.terms.size() == 1 &&
+							eq.lhs.terms.at(0) == Product{1.0, std::vector<std::string>{var}}) {
+						if (debug) std::cerr << "\x1B[1m[" << stage << "] Substituting " << var;
+						if (auto value = eq.rhs.evaluate()) {
+							if (debug) std::cerr << " = " << *value;
+							assignment.emplace_back(
+								eq.lhs.terms.at(0).vars.at(0),
+								*value
+							);
+							assigned[i] = true;
+						}
+						if (debug) std::cerr << "\x1B[0m" << std::endl;
+
+						for (size_t j = 0; j < eqs.size(); j++) {
+							if (j == i || assigned[j]) continue;
+							eqs[j].substitute_inplace(var, eq.rhs);
+							eqs[j].simplify();
+						}
+						substituted[i] = true;
+						goto success_try_again;
+					}
+				}
+			}
+		}
+	}
+
+	std::stringstream errmsg;
+	errmsg << "Equations could not be solved:\n";
+	bool have_error = false;
+
+	// If we arrive here, no more equations could be solved
+	for (size_t i = 0; i < eqs.size(); i++) {
+		if (assigned[i]) continue;
+		have_error = true;
+		errmsg << "  " << eqs[i] << '\n';
+	}
+
+	if (have_error) return errmsg.str();
+	else return assignment;
+}
+
+void System::substitute_inplace(const std::string &var, Sum sum) {
+	for (Equation &eq : eqs) {
+		eq.substitute_inplace(var, sum);
+	}
+}
+
+
+std::ostream& operator<<(std::ostream &os, const Product &p) {
+	bool first = true;
+	if (p.cnst != 1 || p.vars.empty()) {
+		os << p.cnst;
+		first = false;
+	}
+
+	for (const std::string &v : p.vars) {
+		if (first) first = false;
+		else os << "*";
+		os << v;
+	}
+
+	return os;
+}
+
+std::ostream& operator<<(std::ostream &os, const Sum &s) {
+	if (s.terms.empty()) {
+		os << "0";
+		return os;
+	}
+
+	bool first = true;
+	for (const Product &p : s.terms) {
+		if (first) first = false;
+		else os << " + ";
+		os << p;
+	}
+	return os;
+}
+
+std::ostream& operator<<(std::ostream &os, const Equation &eq) {
+	os << eq.lhs << " = " << eq.rhs;
+	return os;
+}
+
+std::ostream& operator<<(std::ostream &os, const System &sys) {
+	bool first = true;
+	for (const Equation &eq : sys.eqs) {
+		if (first) first = false;
+		else os << '\n';
+		os << eq;
+	}
+	return os;
+}
+
+// vim: set sw=4 ts=4 noet:
diff --git a/eqsystem_solve.h b/eqsystem_solve.h
new file mode 100644
index 0000000..34cd4cc
--- /dev/null
+++ b/eqsystem_solve.h
@@ -0,0 +1,62 @@
+#pragma once
+
+#include <optional>
+#include <string>
+#include <variant>
+#include <vector>
+
+
+struct Product {
+	double cnst;
+	std::vector<std::string> vars;
+
+	const Product times(double x) const;
+	void times_inplace(double x);
+	const Product times(const Product &p) const;
+	bool contains(const std::string &var) const;
+	const Product without(const std::string &var) const;
+
+	bool operator==(const Product &other) const;
+};
+
+struct Sum {
+	std::vector<Product> terms;
+
+	const Sum times(double x) const;
+	const Sum substitute(const std::string &var, Sum sum) const;
+	void simplify();
+	// Returns std::nullopt if there are variables in the terms
+	std::optional<double> evaluate() const;
+};
+
+// An equation is not allowed to contain a variable name more than once.
+struct Equation {
+	Sum lhs, rhs;
+
+	std::optional<std::string> validate() const;
+	const Equation substitute(const std::string &var, Sum sum) const;
+	void substitute_inplace(const std::string &var, Sum sum);
+	int num_vars() const;
+	std::vector<std::string> all_vars() const;
+	bool isolatable(const std::string &target) const;
+	void isolate_left(const std::string &target);
+	void simplify();
+};
+
+struct System {
+	std::vector<Equation> eqs;
+
+	std::variant<
+		std::string,  // description of the error
+		std::vector<std::pair<std::string, double>>  // result assignment
+	> solve_inplace(bool debug);
+
+	void substitute_inplace(const std::string &var, Sum sum);
+};
+
+std::ostream& operator<<(std::ostream &os, const Product &p);
+std::ostream& operator<<(std::ostream &os, const Sum &s);
+std::ostream& operator<<(std::ostream &os, const Equation &eq);
+std::ostream& operator<<(std::ostream &os, const System &sys);
+
+// vim: set sw=4 ts=4 noet:
diff --git a/image.cpp b/image.cpp
new file mode 100644
index 0000000..675e22f
--- /dev/null
+++ b/image.cpp
@@ -0,0 +1,56 @@
+#include "image.h"
+#include <iostream>
+#include <array>
+#include <string>
+#include <cstring>
+#include <cassert>
+#include <cmath>
+#include <png.h>
+
+
+static std::array<uint8_t, 3> render_pixel(double x) {
+	const double red = 1 - std::pow(1 - x, 3);
+	const double green = 1 - std::pow(1 - x, 3);
+	const double blue = 1 - std::pow(1 - x, 6);
+	return {(uint8_t)(red * 255), (uint8_t)(green * 255), (uint8_t)(blue * 255)};
+}
+
+HistImage::HistImage(size_t width, size_t height)
+	: width{width}, height{height}, data(width * height, 0) {}
+
+void HistImage::add(const HistImage &other) {
+	assert(width == other.width && height == other.height);
+	for (size_t i = 0; i < width * height; i++) data[i] += other.data[i];
+}
+
+void HistImage::render(const char *fname) const {
+	uint64_t maxval = 0;
+	for (size_t i = 0; i < width * height; i++) maxval = std::max(maxval, data[i]);
+	const double maxvald = maxval;
+
+	std::vector<uint8_t> buffer(3 * width * height);
+	for (size_t i = 0; i < width * height; i++) {
+		const std::array<uint8_t, 3> pix = render_pixel((double)data[i] / maxvald);
+		memcpy(&buffer[3 * i], &pix[0], 3);
+	}
+
+	png_image png;
+	memset(&png, 0, sizeof png);
+	png.version = PNG_IMAGE_VERSION;
+	png.width = width;
+	png.height = height;
+	png.format = PNG_FORMAT_RGB;
+
+	int ret = png_image_write_to_file(&png, fname, 0, buffer.data(), 3 * width, NULL);
+	if (ret != 0 && (png.warning_or_error & 2)) {
+		std::string errmsg;
+		char *p = (char*)memchr(png.message, ' ', sizeof png.message);
+		if (p != NULL) { errmsg = std::string(png.message, p - png.message); }
+		else errmsg = std::string(png.message, sizeof png.message);
+		std::cerr << "PNG write failed: " << errmsg << std::endl;
+	}
+
+	png_image_free(&png);
+}
+
+// vim: set sw=4 ts=4 noet:
diff --git a/image.h b/image.h
new file mode 100644
index 0000000..2d02fae
--- /dev/null
+++ b/image.h
@@ -0,0 +1,18 @@
+#pragma once
+
+#include <vector>
+#include <cstddef>
+#include <cstdint>
+
+
+struct HistImage {
+	size_t width, height;
+	std::vector<uint64_t> data;
+
+	HistImage(size_t width, size_t height);
+
+	void add(const HistImage &other);
+	void render(const char *fname) const;
+};
+
+// vim: set sw=4 ts=4 noet:
diff --git a/main.cpp b/main.cpp
new file mode 100644
index 0000000..765735c
--- /dev/null
+++ b/main.cpp
@@ -0,0 +1,14 @@
+#include "options.h"
+#include "compute.h"
+#include <iostream>
+
+
+int main(int argc, char **argv) {
+	Options options = parse_options(argc, argv);
+	std::cout << "Options:" << std::endl << options << std::endl;
+	HistImage img{compute_brot(options)};
+	if (options.outfname.value != "/dev/null")
+		img.render(options.outfname.value.data());
+}
+
+// vim: set sw=4 ts=4 noet:
diff --git a/options.cpp b/options.cpp
new file mode 100644
index 0000000..3d7786b
--- /dev/null
+++ b/options.cpp
@@ -0,0 +1,200 @@
+#include "options.h"
+#include <thread>
+#include <cctype>
+#include <cassert>
+
+
+static int64_t ipow(int64_t b, int64_t e) {
+	if (e == 0) return 1;
+	assert(e >= 0);
+	int64_t x = 1;
+	while (e > 0) {
+		if (e & 1) x *= b;
+		b *= b;
+		e >>= 1;
+	}
+	return x;
+}
+
+std::optional<const char*> parse_int_param(int64_t *dest, const char *str) {
+	if (!str || !*str) return "empty string";
+	char *endp;
+	int64_t whole = strtoll(str, &endp, 10);
+
+	int64_t fractional = 0, fract_offset = 0;
+	if (*endp == '.' && isdigit(endp[1])) {
+		char *endp2;
+		fractional = strtoll(endp + 1, &endp2, 10);
+		fract_offset = endp2 - (endp + 1);
+		endp = endp2;
+	}
+
+	int64_t exponent = 0;
+	if ((*endp == 'e' || *endp == 'E') && isdigit(endp[1])) {
+		exponent = strtoll(endp + 1, &endp, 10);
+		// log2 (10^20) ~= 66 > 64
+		if (exponent < 0 || exponent >= 20) return "exponent too large";
+	}
+
+	if (fract_offset > exponent) return "not integer (more fractional digits than exponent)";
+
+	if (*endp != '\0') return "unexpected character";
+	*dest = (whole * ipow(10, fract_offset) + fractional) * ipow(10, exponent - fract_offset);
+	return std::nullopt;
+}
+
+void Options::fill_dimensions() {
+	System sys;
+
+	sys.eqs.push_back(
+		Equation{
+			Sum{{xmax.toProduct()}},
+			Sum{{xmin.toProduct(), cplxwidth.toProduct()}}
+		}
+	);
+
+	sys.eqs.push_back(
+		Equation{
+			Sum{{ymax.toProduct()}},
+			Sum{{ymin.toProduct(), cplxheight.toProduct()}}
+		}
+	);
+
+	sys.eqs.push_back(
+		Equation{
+			Sum{{xmin.toProduct(), xmax.toProduct()}},
+			Sum{{cx.toProduct().times(2.0)}}
+		}
+	);
+
+	sys.eqs.push_back(
+		Equation{
+			Sum{{ymin.toProduct(), ymax.toProduct()}},
+			Sum{{cy.toProduct().times(2.0)}}
+		}
+	);
+
+	sys.eqs.push_back(
+		Equation{
+			Sum{{imgwidth.toProduct().times(cplxheight.toProduct())}},
+			Sum{{cplxwidth.toProduct().times(imgheight.toProduct())}}
+		}
+	);
+
+	auto res = sys.solve_inplace(debug_paramsolver);
+	if (auto errmsg = std::get_if<std::string>(&res)) {
+		std::cerr << *errmsg;
+		exit(1);
+	} else {
+		auto assignment = std::get<1>(res);
+		xmin.from_computed(assignment);
+		xmax.from_computed(assignment);
+		ymin.from_computed(assignment);
+		ymax.from_computed(assignment);
+		cplxwidth.from_computed(assignment);
+		cplxheight.from_computed(assignment);
+		cx.from_computed(assignment);
+		cy.from_computed(assignment);
+		imgwidth.from_computed(assignment);
+		imgheight.from_computed(assignment);
+	}
+}
+
+void Options::validate() {
+	if (maxiter.value < 0) {
+		std::cerr << "Negative maxiter" << std::endl;
+		exit(1);
+	}
+
+	if (numorbits.value < 0) {
+		std::cerr << "Negative numorbits" << std::endl;
+		exit(1);
+	}
+
+	if (num_threads.value <= 0) {
+		std::cerr << "Non-positive number of threads" << std::endl;
+		exit(1);
+	}
+
+	if (orbit_burnin.value < 0) {
+		std::cerr << "Negative orbit_burnin" << std::endl;
+		exit(1);
+	}
+
+	if (!outfname.valid()) {
+		std::cerr << "Output filename (-o) not given" << std::endl;
+		exit(1);
+	}
+
+	if (imgwidth.value <= 0 || imgheight.value <= 0) {
+		std::cerr << "Zero or negative image size; check window options" << std::endl;
+		exit(1);
+	}
+}
+
+std::ostream& operator<<(std::ostream &os, const Options &opts) {
+	os << "window: [" << opts.xmin.value << "," << opts.xmax.value << "] (" << opts.cplxwidth.value << ") x [" << opts.ymin.value << "," << opts.ymax.value << "] (" << opts.cplxheight.value << ") around (" << opts.cx.value << "," << opts.cy.value << ")\n";
+	os << "image: " << opts.imgwidth.value << " x " << opts.imgheight.value << "; ";
+	os << "maxiter=" << opts.maxiter.value << "; ";
+	os << "numorbits=" << opts.numorbits.value << "; ";
+	os << "num_threads=" << opts.num_threads.value << "; ";
+	os << "algorithm=" << opts.algorithm.value;
+	return os;
+}
+
+void usage(const char *argv0) {
+	std::cout << "Usage: " << argv0 << " [OPTIONS]\n";
+	std::cout << "  -L [float]    Left boundary of window in C (xmin)\n";
+	std::cout << "  -R [float]    Right boundary of window in C (xmax)\n";
+	std::cout << "  -T [float]    Top boundary of window in C (ymax)\n";
+	std::cout << "  -B [float]    Bottom boundary of window in C (ymin)\n";
+	std::cout << "  -W [float]    Width of window in C (cplxwidth)\n";
+	std::cout << "  -H [float]    Height of window in C (cplxheight)\n";
+	std::cout << "  -X [float]    Horizontal middle of window in C (cx)\n";
+	std::cout << "  -Y [float]    Vertical middle of window in C (cy)\n";
+	std::cout << "  -w [int]      Image width (imgwidth)\n";
+	std::cout << "  -h [int]      Image height (imgheight)\n";
+	std::cout << "  -m [int]      Number of iterations before escape (maxiter)\n";
+	std::cout << "  -n [int]      Number of orbits to trace (numorbits)\n";
+	std::cout << "  -t [int]      Number of threads (num_threads)\n";
+	std::cout << "  -o [string]   Output file name (PNG file) (outfname)\n";
+	std::cout << "  -al [string]  Algorithm to use (see compute.cpp) (algorithm)\n";
+	std::cout << "  --help        This text\n";
+	std::cout << "\n";
+	std::cout << "  -ds           Debug window equation solving\n";
+	std::cout << std::flush;
+}
+
+Options parse_options(int argc, char **argv) {
+	Options options;
+	for (int i = 1; i < argc; i++) {
+		if (strcmp(argv[i], "-L") == 0) options.xmin.parse(argv[i+1], "-L"), i++;
+		else if (strcmp(argv[i], "-R") == 0) options.xmax.parse(argv[i+1], "-R"), i++;
+		else if (strcmp(argv[i], "-T") == 0) options.ymax.parse(argv[i+1], "-T"), i++;
+		else if (strcmp(argv[i], "-B") == 0) options.ymin.parse(argv[i+1], "-B"), i++;
+		else if (strcmp(argv[i], "-W") == 0) options.cplxwidth.parse(argv[i+1], "-W"), i++;
+		else if (strcmp(argv[i], "-H") == 0) options.cplxheight.parse(argv[i+1], "-H"), i++;
+		else if (strcmp(argv[i], "-X") == 0) options.cx.parse(argv[i+1], "-X"), i++;
+		else if (strcmp(argv[i], "-Y") == 0) options.cy.parse(argv[i+1], "-Y"), i++;
+		else if (strcmp(argv[i], "-w") == 0) options.imgwidth.parse(argv[i+1], "-w"), i++;
+		else if (strcmp(argv[i], "-h") == 0) options.imgheight.parse(argv[i+1], "-h"), i++;
+		else if (strcmp(argv[i], "-m") == 0) options.maxiter.parse(argv[i+1], "-m"), i++;
+		else if (strcmp(argv[i], "-n") == 0) options.numorbits.parse(argv[i+1], "-n"), i++;
+		else if (strcmp(argv[i], "-ds") == 0) options.debug_paramsolver = true;
+		else if (strcmp(argv[i], "-t") == 0) options.num_threads.parse(argv[i+1], "-t"), i++;
+		else if (strcmp(argv[i], "-o") == 0) options.outfname.parse(argv[i+1], "-o"), i++;
+		else if (strcmp(argv[i], "-al") == 0) options.algorithm.parse(argv[i+1], "-o"), i++;
+		else if (strcmp(argv[i], "--help") == 0) { usage(argv[0]); exit(0); }
+		else {
+			std::cerr << "Unexpected argument '" << argv[i] << "'" << std::endl;
+			exit(1);
+		}
+	}
+
+	options.fill_dimensions();
+	options.validate();
+
+	return options;
+}
+
+// vim: set sw=4 ts=4 noet:
diff --git a/options.h b/options.h
new file mode 100644
index 0000000..831af3e
--- /dev/null
+++ b/options.h
@@ -0,0 +1,137 @@
+#pragma once
+
+#include "eqsystem_solve.h"
+#include <iostream>
+#include <optional>
+#include <thread>
+#include <cstring>
+#include <cstdlib>
+#include <cstdint>
+#include <cassert>
+#include <unistd.h>
+
+
+template <typename T>
+struct Param {
+	bool given = false, computed = false;
+	T value;  // valid if given || computed
+	std::string name;
+
+	Param(std::string name) : name{name} {}
+	Param(std::string name, T value) : computed{true}, value{std::move(value)}, name{name} {}
+
+	Product toProduct() const;
+	void from_computed(const std::vector<std::pair<std::string, double>> &assign);
+	void parse(const char *str, const char *option);
+	void set(T value);
+	bool valid() const;
+};
+
+template <typename T>
+Product Param<T>::toProduct() const {
+	static_assert(std::is_integral_v<T> || std::is_floating_point_v<T>);
+	if (given || computed) return Product{static_cast<double>(value), {}};
+	else return Product{1.0, {name}};
+}
+
+template <typename T>
+void Param<T>::from_computed(const std::vector<std::pair<std::string, double>> &assign) {
+	if (given || computed) return;
+	for (const auto &p : assign) {
+		if (p.first == name) {
+			computed = true;
+			if constexpr (std::is_integral_v<T>) {
+				value = p.second;
+				if (p.second != value) {
+					std::cerr << "Non-integral value computed for integral option " << name << ": " << p.second << std::endl;
+					exit(1);
+				}
+			} else {
+				static_assert(std::is_same_v<T, double>);
+				value = p.second;
+			}
+			value = p.second;
+			return;
+		}
+	}
+
+	if (!computed) {
+		std::cerr << "Value for option " << name << " ambiguous! Was a dimension zero?" << std::endl;
+		exit(1);
+	}
+}
+
+std::optional<const char*> parse_int_param(int64_t *dest, const char *str);
+
+template <typename T>
+void Param<T>::parse(const char *str, const char *option) {
+	if (str == nullptr) {
+		std::cerr << "Missing argument for option '" << option << "'" << std::endl;
+		exit(1);
+	}
+
+	if constexpr (std::is_same_v<T, double>) {
+		char *endp;
+		value = strtod(str, &endp);
+		if (!*str || *endp) {
+			std::cerr << "Could not parse floating-point value from argument '" << str << "'" << std::endl;
+			exit(1);
+		}
+		given = true;
+	} else if constexpr (std::is_same_v<T, std::string>) {
+		value = str;
+		given = true;
+	} else {
+		static_assert(std::is_same_v<T, int> || std::is_same_v<T, int64_t>);
+		int64_t res;
+		if (std::optional<const char*> err = parse_int_param(&res, str)) {
+			std::cerr << "Could not parse integer from argument '" << str << "': " << *err << std::endl;
+			exit(1);
+		} else {
+			value = (T)res;
+			given = true;
+		}
+	}
+}
+
+template <typename T>
+void Param<T>::set(T value) {
+	computed = true;
+	this->value = std::move(value);
+}
+
+template <typename T>
+bool Param<T>::valid() const {
+	return given || computed;
+}
+
+struct Options {
+	// Mutually-inferrable parameters
+	Param<double> xmin{"xmin"}, xmax{"xmax"}, ymin{"ymin"}, ymax{"ymax"};
+	Param<double> cplxwidth{"cplxwidth"}, cplxheight{"cplxheight"};
+	Param<double> cx{"cx"}, cy{"cy"};
+	Param<int> imgwidth{"imgwidth"}, imgheight{"imgheight"};
+
+	// Normal parameters
+	Param<int> maxiter{"maxiter", 1024};
+	Param<int64_t> numorbits{"numorbits", 1000000};
+	Param<int> num_threads{"num_threads", (int)std::thread::hardware_concurrency()};
+	Param<int> orbit_burnin{"orbit_burnin", 4};
+	bool debug_paramsolver = false;
+	Param<std::string> outfname{"outfname"};
+	Param<std::string> algorithm{"algorithm", "vectorised"};
+
+	// Infer the mutually-inferrable parameters
+	void fill_dimensions();
+
+	// Check presence of normal parameters and consistency of all parameters
+	void validate();
+};
+
+std::ostream& operator<<(std::ostream &os, const Options &opts);
+
+void usage(const char *argv0);
+
+Options parse_options(int argc, char **argv);
+
+// vim: set sw=4 ts=4 noet: