diff options
author | Tom Smeding <tom@tomsmeding.com> | 2025-09-03 23:31:40 +0200 |
---|---|---|
committer | Tom Smeding <tom@tomsmeding.com> | 2025-09-03 23:37:53 +0200 |
commit | 7e60437d3c064bca402d486be967c43bf4326067 (patch) | |
tree | 8b4307b90c942861176d9c705cd760b9fe7b6c6b |
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).
-rw-r--r-- | .gitignore | 6 | ||||
-rw-r--r-- | Makefile | 35 | ||||
-rwxr-xr-x | bench.sh | 2 | ||||
-rw-r--r-- | compute.cpp | 191 | ||||
-rw-r--r-- | compute.h | 9 | ||||
-rw-r--r-- | eqsystem_solve.cpp | 330 | ||||
-rw-r--r-- | eqsystem_solve.h | 62 | ||||
-rw-r--r-- | image.cpp | 56 | ||||
-rw-r--r-- | image.h | 18 | ||||
-rw-r--r-- | main.cpp | 14 | ||||
-rw-r--r-- | options.cpp | 200 | ||||
-rw-r--r-- | options.h | 137 |
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: @@ -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: |