Mandelbrot Orbit Differences

README.md

Obtain the JSON and npy libraries.

https://github.com/nlohmann/json

https://github.com/rogersce/cnpy

Compile the binary.

g++ -O3 -mavx2 -o orbits orbits.cc -lcnpy

Allow this to run and accumulate points. It runs indefinitely and produces the output file every time it prints a line to the terminal.

./orbits params.json images.npy

Regularly process the output into a sequence of images. With the given params.json file, this will produce 16 images named diff_stack_%08d.png

python3 process_stack.py images.npy diff_stack

orbits.cc

#include <xmmintrin.h>
#include <immintrin.h>
#include <smmintrin.h>
#include <complex>
#include <fstream>
#include <iostream>
#include <cnpy.h>
#include <vector>
#include "nlohmann/json.hpp"

using namespace std::complex_literals;
using json = nlohmann::json;
int main(int argc, char* argv[]) {
	typedef long long S;
	typedef double T;
	typedef long long A;
	typedef int I;
	const int components = 4;

	typedef T vec __attribute__((vector_size(sizeof(T) * components)));
	typedef I veci __attribute__((vector_size(sizeof(I) * components)));

	if(argc <= 2) {
		return -1;
	}

	const std::string& out_name = argv[2];
	const std::string& in_name = argv[1];



	json args;
	std::ifstream infile;
	infile.open(in_name);
	infile >> args;
	infile.close();


	const long long iterations = args["iterations"];

	const size_t width = args["width"];
	const size_t height = args["height"];
	T x_0, x_1, y_0, y_1;
	x_0 = args["canvas range"][0][0];
	y_0 = args["canvas range"][0][1];
	x_1 = args["canvas range"][1][0];
	y_1 = args["canvas range"][1][1];

	const size_t size = width * height;
	const size_t channels = args["channels"];
	A* img;
	img = new A[width * height * channels]();
	S p_0, q_0, p_1, q_1;
	p_0 = args["circle range"][0][0];
	q_0 = args["circle range"][0][1];
	p_1 = args["circle range"][1][0];
	q_1 = args["circle range"][1][1];

	T radius = args["radius"];
	T threshold = args["threshold"];

	S p_2, q_2;
	p_2 = args["starting rational"][0];
	q_2 = args["starting rational"][1];

	S s = 0;

	vec x_scale_vec;
	vec y_scale_vec;
	vec x_0_vec;
	vec y_0_vec;
	veci zero_vec;
	veci width_vec;
	veci height_vec;

	for(int i = 0; i < components; i++) {
		x_0_vec[i] = x_0;
		y_0_vec[i] = y_0;
		x_scale_vec[i] = width / (x_1 - x_0);
		y_scale_vec[i] = height / (y_1 - y_0);
		zero_vec[i] = 0;
		width_vec[i] = width;
		height_vec[i] = height;
	}

	S total_orbits;

	for(S total_orbits = 0; !args.contains("total orbits") || total_orbits < args["total orbits"]; total_orbits += components) {
		vec c_real;
		vec c_imag;
		veci i_vec;
		veci qs2_vec;

		S iterations2 = 0;

		for(int i = 0; i < components; i++) {
			i_vec[i] = 0;
			S p = p_0 * q_2 + p_1 * p_2;
			S q = q_0 * q_2 + q_1 * p_2;

			s += q;

			if(q > iterations2) {
				iterations2 = q;
			}

			{
				T t = 2 * M_PI * T(p) / T(q);
				std::complex<T> mu = std::polar(radius, t);
				const auto P = mu / T(2);
				const auto c = P - P * P;
				qs2_vec[i] = q * iterations;
				c_real[i] = c.real();
				c_imag[i] = c.imag();
			}
		}

		iterations2 *= iterations;

		vec z_real = c_real;
		vec z_imag = c_imag;

		vec z0_real[channels];
		vec z0_imag[channels];

		for(int i = 0; i < channels; i++) {
			vec z_real_sq = z_real * z_real;
			vec z_imag_sq = z_imag * z_imag;
			vec z_real_imag = z_real * z_imag;
			vec z_abs_sq = z_real_sq + z_imag_sq;
			z_real = z_real_sq - z_imag_sq + c_real;
			z_imag = z_real_imag + z_real_imag + c_imag;
			z0_real[i] = z_real;
			z0_imag[i] = z_imag;
		}

		size_t size = width * height;
		{
			for(int i = 0; i < iterations2; i++, i_vec += 1) {
				vec z_real_sq = z_real * z_real;
				vec z_imag_sq = z_imag * z_imag;
				vec z_real_imag = z_real * z_imag;
				vec z_abs_sq = z_real_sq + z_imag_sq;
				z_real = z_real_sq - z_imag_sq + c_real;
				z_imag = z_real_imag + z_real_imag + c_imag;

				for(int j = 0; j < channels; j++) {
					A* img2 = &img[j * size];
					vec x = z_real - z0_real[(i - j - 1) % channels];
					vec y = z_imag - z0_imag[(i - j - 1) % channels];

					veci x_s = __builtin_convertvector((x - x_0_vec) * y_scale_vec, veci);
					veci y_s = __builtin_convertvector((y - y_0_vec) * y_scale_vec, veci);
					veci index = y_s * width_vec + x_s;
					veci mask = (zero_vec <= x_s) && (x_s < width_vec) && (zero_vec <= y_s) && (y_s < height_vec) && (i_vec < qs2_vec);

					for(int k = 0; k < components; k++) {
						if(mask[k]) {
							size_t index_k = index[k];
							img2[index_k] ++;
						}
					}
				}

				z0_real[i % channels] = z_real;
				z0_imag[i % channels] = z_imag;
			}
		}

		if(args.contains("refresh iterations") && s >= args["refresh iterations"]) {
			json status;
			status["current rational"] = {p_2, q_2};
			std::ofstream o("status.json");
			o << status << std::endl;
			o.close();
			s = 0;
			std::cerr << p_2 << "/" << q_2 << std::endl;
			cnpy::npy_save(out_name, img, {channels, height, width}, "w");
		}

		{
			S t = ((p_2 / q_2) * 2 + 1) * q_2 - p_2;
			p_2 = q_2;
			q_2 = t;
		}

	}

	cnpy::npy_save(out_name, img, {channels, height, width}, "w");
	delete [] img;

	return 0;
}

params.json

{"iterations": 1, "width": 2048, "height": 2048, "channels": 16, "threshold": 2, 
    "radius":1,
    "starting rational": [0,1],
    "refresh iterations": 268435456,
     "endpoints": [false, false], "canvas range": [[-1,-1], [1, 1]], "circle range": [[0,1],[1,2]]}

process_stack.py

import datashader.transfer_functions as tf
import datashader.utils
import colorcet
import numpy as np
import xarray
import sys


def main():
    # Usage: process_bin.py images.npy /path/to/out
    image_stack = np.load(sys.argv[1])
    out_prefix = sys.argv[2]
    cmap = colorcet.fire
    for i in range(image_stack.shape[0]):
        agg = image_stack[i]
        agg += agg[::-1]
        agg = agg.T
        if np.max(agg) < 2**28:
            how = 'eq_hist'
        else:
            how = 'log'
        agg = xarray.DataArray(agg)
        img = tf.shade(agg, cmap=cmap, how=how)
        img = tf.set_background(img, 'black')
        with open(f'{out_prefix}_{i:08}.png', 'wb') as outfile:
            img.to_pil().save(outfile)


if __name__ == '__main__':
    main()