summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--.gitignore6
-rw-r--r--Makefile35
-rwxr-xr-xbench.sh2
-rw-r--r--compute.cpp191
-rw-r--r--compute.h9
-rw-r--r--eqsystem_solve.cpp330
-rw-r--r--eqsystem_solve.h62
-rw-r--r--image.cpp56
-rw-r--r--image.h18
-rw-r--r--main.cpp14
-rw-r--r--options.cpp200
-rw-r--r--options.h137
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:
generated by cgit v1.2.3-70-g09d2 (git 2.51.0) at 2025-10-31 16:45:48 +0100