init

public/main.js

@@ -0,0 +1,544 @@
+// CC BY-NC-SA 4.0, gtrrebel, shikhin, zgrep
+
+var dr_canvas, p_canvas
+var dr_ctx, p_ctx
+
+const high_resolution = 1000, low_resolution = 100
+const enhance_resolutions = [low_resolution, 400, 800, high_resolution, 1.5 * high_resolution]
+const enhance_text = ["enhance", "enhance more", "enhance further", "ENHANCE", "no more enhance"]
+
+var enhance_level = 0
+var dr_resolution = high_resolution
+
+var enhance_button, enhance_enabled = false
+
+var dim = 2
+var matrix_entries = [];//[-3, 0, 0, 0, 0, 0, 0, -1,
+						//0, 0, -2, 0, 0, 0, 0, -1,
+						//0, 0, 0, 0, 1, -2, 0, -2,
+						//0, -1, 0, -1, 0, -2, 2, 2]
+var tangents = []
+
+const tangent_style = "rgb(30,144,255)"
+const tangent_thickness = 3
+const tangent_thinness = 3
+
+const p_radius_const =  0.06
+const p_scale = 3.0
+const p_font = "40px serif"
+var p_radius, p_radius_sq
+
+var zoom_scale = 1
+var dr_scale = zoom_scale * dim * p_scale
+
+var mouse_moving = false, touch_moving = false, moving_id = 0
+var touch_identifier
+var moving_origin = []
+var moving_epsilon = 0.04
+
+function p_to_canv_coords(x, y) {
+	return [p_canvas.width * (x + p_scale) / (2 * p_scale), p_canvas.height * (-y + p_scale) / (2 * p_scale)]
+}
+
+function p_coords(i, j) {
+	return [p_scale * ((i * 2.0 / p_canvas.width) - 1.0), -p_scale * ((j * 2.0 / p_canvas.height) - 1.0)]
+}
+
+function dr_to_canv_coords(x, y) {
+	return [dr_canvas.width * (x + dr_scale) / (2 * dr_scale), dr_canvas.height * (-y + dr_scale) / (2 * dr_scale)]
+}
+
+function init() {
+	enhance_button = document.getElementById("enhance")
+
+	dr_canvas = document.getElementById("drawing_canvas")
+	dr_ctx = dr_canvas.getContext("2d")
+
+	p_canvas = document.getElementById("picker_canvas")
+	p_ctx = p_canvas.getContext("2d")
+
+	p_canvas.style.width = p_canvas.style.height = "400px"
+	p_canvas.width = p_canvas.height = 800
+
+	dr_canvas.style.width = dr_canvas.style.height = "1000px"
+	dr_canvas.width = dr_canvas.height = dr_resolution
+
+	Module.set_resolution(dr_resolution)
+
+	p_radius = p_canvas.width * p_radius_const
+	p_radius_sq = p_radius * p_radius
+
+	p_canvas.addEventListener("mousedown", mouse_down)
+	p_canvas.addEventListener("mousemove", mouse_move)
+	p_canvas.addEventListener("mouseup", mouse_up)
+	p_canvas.addEventListener("mouseout", mouse_up)
+
+	p_canvas.addEventListener("touchstart", touch_start)
+	p_canvas.addEventListener("touchmove", touch_move)
+	p_canvas.addEventListener("touchend", touch_end)
+	p_canvas.addEventListener("touchcancel", touch_end)
+
+	dim_update()
+}
+
+// const known_letters = 'abcdefghijklmnopqrstuvwxyzαβγδεζηθικλμνξπρστυφχψωABCDEFGHIJKLMNOPQRSTUVWXYZΓΔΘΛΞΠΣΦΨΩςאבגדהזחטכלמנעפצקשת';
+function matrix_letter(i) {
+	return 'C_{' + Math.floor(i / dim) + ',' + (i % dim) + '}'
+}
+
+function p_draw() {
+	p_ctx.clearRect(0, 0, p_canvas.width, p_canvas.height)
+
+	// axes
+	p_ctx.beginPath()
+	p_ctx.moveTo(p_canvas.width / 2, 0)
+	p_ctx.lineTo(p_canvas.width / 2, p_canvas.height)
+	p_ctx.stroke()
+
+	p_ctx.beginPath()
+	p_ctx.moveTo(0, p_canvas.height / 2)
+	p_ctx.lineTo(p_canvas.width, p_canvas.height / 2)
+	p_ctx.stroke()
+
+	for (let i = 0; i < 2 * dim * dim; i += 2) {
+		const [x, y] = p_to_canv_coords(matrix_entries[i], matrix_entries[i + 1])
+		p_ctx.beginPath()
+		p_ctx.arc(x, y, p_radius, 0, 2 * Math.PI, false)
+		p_ctx.fillStyle = "rgba(149, 165, 166,0.9)"
+		p_ctx.fill()
+
+		p_ctx.fillStyle = "black"
+		p_ctx.font = p_font
+		p_ctx.fillText(matrix_letter(i / 2), x + 1.5 * p_radius, y)
+	}
+
+	// text input grid
+	update_grid()
+}
+
+function set_enhance_level(level) {
+	enhance_level = level
+	enhance_button.textContent = enhance_text[level]
+
+	dr_resolution = enhance_resolutions[level]
+	dr_canvas.width = dr_canvas.height = dr_resolution
+	Module.set_resolution(dr_resolution)
+}
+
+function enhance_disable(dim_switch = false) {
+	if (dim_switch) {
+		enhance_button.textContent = enhance_text[0]
+		dr_resolution = high_resolution
+		dr_canvas.width = dr_canvas.height = dr_resolution
+		Module.set_resolution(dr_resolution)
+	}
+
+	if (!enhance_enabled) { return }
+
+	enhance_button.setAttribute('disabled', true)
+	enhance_enabled = false
+
+	if (enhance_level > 0 && !dim_switch) {
+		set_enhance_level(0)
+	}
+
+	dr_draw()
+}
+
+function enhance_enable() {
+	if (dim < 5) { return }
+	enhance_button.removeAttribute('disabled')
+	enhance_enabled = true
+}
+
+function enhance() {
+	if (enhance_level == enhance_resolutions.length - 1) { return }
+
+	enhance_button.textContent = "enhancing..."
+	setTimeout(function() {
+		set_enhance_level(enhance_level + 1)
+		dr_draw()
+		enhance_button.textContent = enhance_text[enhance_level]
+	}, 10);
+}
+
+function to_float_arr(arr) {
+	const res = new Float32Array(arr.length)
+	for (let i = 0; i < arr.length; i++) {
+		res[i] = arr[i]
+	}
+	return res
+}
+
+function transfer_matrix_to_heap() {
+	const typed_array = to_float_arr(matrix_entries)
+	const heap_pointer = Module._malloc(typed_array.length * typed_array.BYTES_PER_ELEMENT)
+
+	Module.HEAPF32.set(typed_array, heap_pointer >> 2)
+	return heap_pointer
+}
+
+function dr_draw() {
+	tangents = []
+
+	let heap_array
+	try {
+		heap_array = transfer_matrix_to_heap()
+		Module.set_matrix(dim, heap_array)
+	} finally {
+		Module._free(heap_array)
+	}
+
+	let arr = new Uint8ClampedArray(Module.HEAPU8.buffer, Module.get_output_buffer(),
+									dr_canvas.width * dr_canvas.height * 4)
+	var imageData = dr_ctx.createImageData(dr_canvas.width, dr_canvas.height)
+	imageData.data.set(arr)
+
+	dr_ctx.clearRect(0, 0, dr_canvas.width, dr_canvas.height);
+	dr_ctx.putImageData(imageData, 0, 0)
+
+	const [or_x, or_y] = dr_to_canv_coords(0, 0)
+	dr_ctx.strokeStyle = tangent_style
+
+	dr_ctx.beginPath()
+	
+	if (dim >= 5 && enhance_level == 0) {
+		dr_ctx.lineWidth = tangent_thinness
+	} else {
+		dr_ctx.lineWidth = tangent_thickness
+	}
+	if (dim < 5 || (dim >= 5 && enhance_level > 0)) {
+		dr_ctx.fillStyle = tangent_style
+		dr_ctx.arc(or_x, or_y, tangent_thickness, 0, 2 * Math.PI, false)
+		dr_ctx.fill()
+	}
+
+	for (let i = 0; i < tangents.length; i++) {
+		dr_ctx.moveTo(or_x, or_y)
+		dr_ctx.lineTo(tangents[i][0], tangents[i][1])
+	}
+	dr_ctx.stroke()
+}
+
+function add_tangent(a, b) {
+	tangents.push(dr_to_canv_coords(a, b))
+}
+
+function update_dr_scale() {
+	dr_scale = zoom_scale * dim * p_scale
+}
+
+function update_zoom() {
+	update_dr_scale()
+	Module.set_zoom(zoom_scale)
+	dr_draw()
+}
+
+function zoom_in_dr() {
+	zoom_scale *= 0.8
+	update_zoom()
+}
+
+function zoom_out_dr() {
+	zoom_scale /= 0.8
+	update_zoom()
+}
+
+function randomize_matrix() {
+	matrix_entries = []
+	for (let i = 0; i < dim * dim; i++) {
+		matrix_entries.push((6.0 * Math.random()) - 3.0)
+		matrix_entries.push((6.0 * Math.random()) - 3.0)
+	}
+
+	p_draw()
+	dr_draw()
+}
+
+async function copy_matrix() {
+	matrix = "["
+	for (let i = 0; i < dim; i++) {
+		matrix += "["
+		for (let j = 0; j < dim; j++) {
+			matrix += matrix_entries[2 * (dim * i + j)].toFixed(3)
+			matrix += "+"
+			matrix += matrix_entries[2 * (dim * i + j) + 1].toFixed(3)
+			matrix += "i"
+			if (j < dim - 1) matrix += ","
+		}
+		matrix += "]"
+		if (i < dim - 1) matrix += "\n"
+	}
+	matrix += "]\n"
+	matrix += "base64 = "
+	matrix += btoa(JSON.stringify(matrix_entries))
+
+	let copy_timeout = undefined;
+	dom_copy = document.getElementById("copy_button")
+	try {
+		await navigator.clipboard.writeText(matrix);
+		dom_copy.textContent = 'copied';
+	} catch (e) {
+		dom_copy.textContent = 'failed to copy';
+	} finally {
+		copy_timeout && clearTimeout(copy_timeout);
+		copy_timeout = setTimeout(function() {
+			dom_copy.textContent = 'copy';
+			copy_timeout = undefined;
+		}, 1500);
+	}
+}
+
+function enter_matrix() {
+	var matrix = prompt('base64?')
+
+	try {
+		entries = JSON.parse(atob(matrix))
+		assert(Array.isArray(entries))
+		assert(entries.length == 2 * dim * dim)
+		for (let i = 0; i < 2 * dim * dim; i++) {
+			assert(!isNaN(entries[i]) && -3 < entries[i] && entries[i] < 3)
+		}
+
+		matrix_entries = entries
+		dr_draw()
+	} catch(error) { }
+}
+
+function dim_update() {
+	var old_dim = dim
+	dim = Number(document.getElementById("dim").value)
+	update_dr_scale()
+
+	matrix_entries = matrix_entries.slice(0, 2 * dim * dim)
+
+	let diff = 2 * dim * dim - matrix_entries.length
+	for (let i = 0; i < diff; i += 2) {
+		matrix_entries.push((6.0 * Math.random()) - 3.0)
+		matrix_entries.push((6.0 * Math.random()) - 3.0)
+	}
+
+	let latex_matrix = [];
+	for (let row = 0; row < dim; row++) {
+		let current = [];
+		latex_matrix.push(current);
+		for (let col = 0; col < dim; col++) {
+			current.push('{' + matrix_letter(row*dim + col) + '}');
+		}
+	}
+	latex_matrix = ',\\hspace{2em} A + iB = \\begin{bmatrix}' + latex_matrix.map(row => row.join('&')).join('\\\\') + '\\end{bmatrix} =\\ '
+	katex.render(latex_matrix, document.getElementById('matrix'))
+	make_grid()
+
+	if (old_dim < 5 && dim >= 5) {
+		dr_resolution = low_resolution
+		dr_canvas.width = dr_canvas.height = dr_resolution
+		Module.set_resolution(dr_resolution)
+	}
+
+	if (old_dim >= 5 && dim < 5) {
+		dr_resolution = high_resolution
+		dr_canvas.width = dr_canvas.height = dr_resolution
+		Module.set_resolution(dr_resolution)
+	}
+
+	enhance_disable(old_dim >= 5 && dim < 5)
+	enhance_enable()
+
+	p_draw()
+	dr_draw()
+}
+
+var input_timers = {};
+
+function input_update(e) {
+	const input = e.srcElement;
+	const row = input.dataset.row;
+	const col = input.dataset.col;
+	const i = row * dim * 2 + col * 2;
+	try {
+		let [real, imag] = input.value.split('+');
+		real = Number(real);
+		imag = Number(imag.replace(/i$/, ''));
+		assert(!isNaN(real) && -3 < real && real < 3);
+		assert(!isNaN(imag) && -3 < imag && imag < 3);
+		matrix_entries[i] = real;
+		matrix_entries[i+1] = imag;
+		p_draw()
+		dr_draw()
+	} catch (error) {
+		update_grid();
+		const old = input.style.borderColor;
+		input.style.borderColor = 'red';
+		input.blur();
+		input_timers[i] = setTimeout(function() {
+			input.style.borderColor = old;
+			delete input_timers[i];
+		}, 1500);
+	}
+}
+
+function make_grid() {
+	for (const [i, timer] of Object.entries(input_timers)) {
+		clearTimeout(timer);
+		delete input_timers[i];
+	}
+	let grid = document.createElement('table');
+	grid.id = 'matrix_table';
+	for (let row = 0; row < dim; row++) {
+		let section = document.createElement('tr');
+		grid.appendChild(section);
+		for (let col = 0; col < dim; col++) {
+			let input = document.createElement('input');
+			input.type = 'text';
+			input.onchange = input_update;
+			input.dataset.row = row;
+			input.dataset.col = col;
+			let td = document.createElement('td');
+			td.appendChild(input);
+			section.appendChild(td);
+		}
+	}
+	document.getElementById('matrix_table').replaceWith(grid);
+}
+
+function update_grid() {
+	let grid = document.getElementById('matrix_table');
+	for (let row = 0; row < dim; row++) {
+		let section = grid.children[row];
+		for (let col = 0; col < dim; col++) {
+			let input = section.children[col].children[0];
+			let i = row*dim*2 + col*2;
+			input.value = matrix_entries[i] + '+' + matrix_entries[i+1] + 'i';
+		}
+	}
+}
+
+// mouse
+function mouse_event_get_coords(e) {
+	rect = p_canvas.getBoundingClientRect(); 
+	x = e.clientX - rect.left; 
+	y = e.clientY - rect.top;
+	
+	return [2 * x, 2 * y]
+}
+
+function touch_event_get_coords(e) {
+	rect = p_canvas.getBoundingClientRect(); 
+	for (let i = 0; i < e.changedTouches.length; i++) {
+		if (e.changedTouches[i].identifier == touch_identifier) {
+			x = e.changedTouches[i].clientX - rect.left; 
+			y = e.changedTouches[i].clientY - rect.top;
+	
+			return [2 * x, 2 * y]
+		}
+	}
+	return false
+}
+
+function input_start(x, y) {
+	for (let i = 0; i < 2 * dim * dim; i += 2) {
+		const [point_x, point_y] = p_to_canv_coords(matrix_entries[i], matrix_entries[i + 1])
+		
+		if (((x - point_x)*(x - point_x) + (y - point_y)*(y - point_y)) <= p_radius_sq) {
+			moving_id = i
+			moving_origin[0] = matrix_entries[i]
+			moving_origin[1] = matrix_entries[i + 1]
+
+			enhance_disable()
+			return true
+		}
+	}
+	return false
+}
+
+function touch_start(e) {
+	if (touch_moving || mouse_moving) { return }
+
+	rect = p_canvas.getBoundingClientRect()
+	for (let i = 0; i < e.touches.length; i++) {
+		x = e.touches[i].clientX - rect.left
+		y = e.touches[i].clientY - rect.top
+		touch_moving = input_start(2 * x, 2 * y)
+		if (touch_moving) {
+			touch_identifier = e.touches[i].identifier
+			return
+		}
+	}
+}
+
+function mouse_down(e) {
+	if (touch_moving) { return }
+	const [x, y] = mouse_event_get_coords(e)
+	mouse_moving = input_start(x, y)
+} 
+
+function input_moving(x, y) {
+	if (x < p_radius || x > (p_canvas.width - p_radius) ||
+		y < p_radius || y > (p_canvas.height - p_radius)) {
+		return
+	}
+
+	[x, y] = p_coords(x, y)
+	matrix_entries[moving_id] = x
+	matrix_entries[moving_id + 1] = y
+	p_draw()
+
+	if (dim <= 2 || (Math.abs(matrix_entries[moving_id] - moving_origin[0]) +
+			Math.abs(matrix_entries[moving_id + 1] - moving_origin[1])) > moving_epsilon) {
+		moving_origin[0] = matrix_entries[moving_id]
+		moving_origin[1] = matrix_entries[moving_id + 1]
+		dr_draw()
+	}
+}
+
+function touch_move(e) {
+	if (touch_moving == false) { return }
+	const ret = touch_event_get_coords(e)
+	if (ret == false) { return }
+	input_moving(ret[0], ret[1])
+}
+
+function mouse_move(e) {
+	if (mouse_moving == false) { return }
+	const [x, y] = mouse_event_get_coords(e)
+	input_moving(x, y)
+}
+
+function input_end(x, y) {
+	if (x < p_radius) { x = p_radius }
+	if (x > (p_canvas.width - p_radius)) { x = p_canvas.width - p_radius }
+	if (y < p_radius) { y = p_radius }
+	if (y > (p_canvas.height - p_radius)) { y = p_canvas.width - p_radius }
+
+	[x, y] = p_coords(x, y)
+	matrix_entries[moving_id] = x
+	matrix_entries[moving_id + 1] = y
+
+	enhance_enable()
+
+	p_draw()
+	dr_draw()
+}
+
+function touch_end(e) {
+	if (touch_moving == true) {
+		const ret = touch_event_get_coords(e)
+		if (ret == false) { return }
+		input_end(ret[0], ret[1])
+		touch_moving = false
+	}
+}
+
+function mouse_up(e) {
+	//for (let i = 0; i < 4; i += 1) {
+	//	document.getElementById("number_"+i).value = matrix_entries[i][0].toFixed(3) + 
+	//												 (matrix_entries[i][1] >= 0 ? "+" : "") + matrix_entries[i][1].toFixed(3) + "i"
+	//}
+
+	if (mouse_moving == true) {
+		const [x, y] = mouse_event_get_coords(e)
+		input_end(x, y)
+		mouse_moving = false
+	}
+}

public/tjsm.html

@@ -0,0 +1,60 @@
+<!DOCTYPE html>
+<html>
+
+<!-- CC BY-NC-SA 4.0, gtrrebel, shikhin, zgrep -->
+
+<head>
+<title>tjsm</title>
+<link rel="stylesheet" href="https://cdn.jsdelivr.net/npm/katex@0.16.9/dist/katex.min.css" integrity="sha384-n8MVd4RsNIU0tAv4ct0nTaAbDJwPJzDEaqSD1odI+WdtXRGWt2kTvGFasHpSy3SV" crossorigin="anonymous">
+<script src="https://cdn.jsdelivr.net/npm/katex@0.16.9/dist/katex.min.js" integrity="sha384-XjKyOOlGwcjNTAIQHIpgOno0Hl1YQqzUOEleOLALmuqehneUG+vnGctmUb0ZY0l8" crossorigin="anonymous"></script>
+<script src="wasm-main.js" type="text/javascript"></script>
+<script src="main.js" type="text/javascript"></script>
+
+<script>
+	Module.postRun = init
+</script>
+
+<style>
+table {
+	display: inline-table;
+	vertical-align: middle;
+}
+canvas {
+	border: 1px solid grey;
+	vertical-align: top;
+	max-width: 100%;
+	max-height: 100%;
+	aspect-ratio: 1 / 1;
+}
+#matrix_table input {
+	width: 280px;
+}
+</style>
+</head>
+
+<body>
+
+<div id="grid_input">
+<label id="d_eq"></label><input type="number" id="dim" onChange="dim_update()" value="2" min="2" style="width: 100px;">
+<span id="matrix"></span><table id="matrix_table"></table>
+</div>
+<script type="text/javascript">katex.render("d =\\ ", document.getElementById('d_eq'))</script>
+
+<br>
+
+<div>
+<canvas id="picker_canvas" width="400" height="400"></canvas>
+<canvas id="drawing_canvas" width="1000" height="1000"></canvas>
+
+
+<button id="enhance" type="button" onclick="enhance()" disabled>enhance</button>
+<button id="randomize_button" type="button" onclick="randomize_matrix()">randomize</button>
+<button id="zoom_in_button" type="button" onclick="zoom_in_dr()">+</button>
+<button id="zoom_out_button" type="button" onclick="zoom_out_dr()">-</button>
+<button id="copy_button" type="button" onclick="copy_matrix()">copy</button>
+<button id="enter_button" type="button" onclick="enter_matrix()">enter matrix</button>
+
+</div>
+
+</body>
+</html>

src/main.cpp

@@ -0,0 +1,548 @@
+// CC BY-NC-SA 4.0, gtrrebel, shikhin, zgrep
+
+#include <cmath>
+#include <complex>
+#include <cstdlib>
+#include <Eigen/Dense>
+#include <emscripten/bind.h>
+#include <emscripten/emscripten.h>
+#include <iostream>
+#include <thread>
+
+const std::complex<float> If(0.0f, 1.0f); 
+
+const std::complex<float> Omega(-1.0f/2.0f, std::sqrt(3.0f) / 2);
+const std::complex<float> Omega_sq(-1.0f/2.0f, -std::sqrt(3.0f) / 2);
+
+const float eps = 1e-12;
+const float PI = 3.141592653589793238463L;
+const float M_2PI = 2*PI;
+
+float p_scale = 3;
+float zoom_scale = 1;
+float dr_scale;
+
+using emscripten::val;
+
+int d;
+
+size_t width, height;
+size_t num_threads;
+
+float *output = nullptr;
+uint32_t *output_buffer = nullptr;
+
+thread_local const val document = val::global("document");
+
+std::array<std::complex<float>, 2> quadratic_solver(float a, float b, float c);
+
+size_t EMSCRIPTEN_KEEPALIVE get_output_buffer() {
+	return (size_t) output_buffer;
+}
+
+float density_default(Eigen::MatrixXcf &A, Eigen::MatrixXcf &B, float a, float b) {
+	Eigen::MatrixXcf first = Eigen::MatrixXcf::Identity(A.rows(), A.cols()) - (((a / (a*a + b*b)) * A) + (b / (a*a + b*b) * B));
+	Eigen::MatrixXcf second = (b * A - a * B).inverse();
+
+	Eigen::ComplexEigenSolver<Eigen::MatrixXcf> eigensolver(first * second, /* computeEigenvectors = */ false);
+	if (eigensolver.info() != Eigen::Success) { return 0; }
+
+	auto eigenvalues = eigensolver.eigenvalues();
+	float value = 0;
+	for(size_t i = 0; i < eigenvalues.rows(); i++) {
+		if (std::abs(eigenvalues[i].imag()) > eps) {
+			value += std::abs(eigenvalues[i].imag());
+		}
+	}
+
+	return value;
+}
+
+float density_2(Eigen::MatrixXcf &A, Eigen::MatrixXcf &B, float a, float b, float precomp_coeffs[7]) {
+	//float tt0 = precomp_coeffs[0], tt1 = precomp_coeffs[1], tt00 = precomp_coeffs[2], tt01 = precomp_coeffs[3],
+	//		tt11 = precomp_coeffs[4];
+	//auto roots = quadratic_solver(tt1*tt1 - tt11 - 2*tt1*b + 2*b*b,
+	//								-2*tt01 + 2*tt0*tt1 - 2*tt1*a - 2*tt0*b + 4*a*b,
+	//								tt0*tt0 - tt00 - 2*tt0*a + 2*a*a);
+	//float value = 0;
+	//for (const auto &root : roots) {
+	//	if (root.imag() > eps) {
+	//		value += root.imag() / ((b * root + a) * std::conj(b * root + a)).real();
+	//	}
+	//}
+	//return value;
+
+	float value = precomp_coeffs[0] + precomp_coeffs[1]*b + precomp_coeffs[2]*b*b + precomp_coeffs[3]*a 
+			+ precomp_coeffs[4]*a*b + precomp_coeffs[5]*a*a;
+	if (value < eps) { return 0.0f; }
+
+	float x = (precomp_coeffs[6]*b*b + precomp_coeffs[7]*a*b + precomp_coeffs[8]*a*a);
+	value /= (x * x);
+	return std::sqrt(value);
+}
+
+float cubic_solver_cheat(float values[4]) {
+	if (std::abs(values[3]) < eps) { 
+		if (std::abs(values[2]) < eps) {
+			return 0;
+		} 
+		
+		auto D = values[1] * values[1] - 4 * values[0] * values[2];
+		if (D > eps) { return 0; }
+
+		return std::sqrt(-D) / (2.0f * values[2]);
+	}
+
+	auto p = -values[2] / (3.0f * values[3]), 
+		q = p*p*p + (values[1]*values[2] - 3.0f*values[0]*values[3]) / (6.0f*(values[3]*values[3])),
+		r = values[1] / (3.0f*values[3]);
+
+	auto rp2 = (r - p*p), h = q*q + rp2*rp2*rp2;
+
+	if (h < eps) { return 0; }
+
+	h = std::sqrt(h);
+	return std::sqrt(3.0f) * std::abs(std::cbrt(q + h) - std::cbrt(q - h)) / (2.0f);	
+}
+
+float density_3(Eigen::MatrixXcf &A, Eigen::MatrixXcf &B, float a, float b, float precomp_coeffs_1[7][7],
+				float precomp_coeffs_2[7][7], float precomp_coeffs_3[7][7], float precomp_coeffs_4[7][7]) {
+	float values[4] = { 0, 0, 0, 0 };
+	float a_powers[7] = { 1, a, 0, 0, 0, 0, 0 };
+	float b_powers[7] = { 1, b, 0, 0, 0, 0, 0 };
+
+	for (size_t i = 2; i < 7; i += 1) {
+		a_powers[i] = a_powers[i - 1] * a;
+		b_powers[i] = b_powers[i - 1] * b;
+	}
+
+	for (size_t i = 0; i < 7; i += 1) {
+		for (size_t j = 0; j < 7; j += 1) {
+			float prod = a_powers[j] * b_powers[i];
+			values[0] += precomp_coeffs_1[i][j] * prod;
+			values[1] += precomp_coeffs_2[i][j] * prod;
+			values[2] += precomp_coeffs_3[i][j] * prod;
+			values[3] += precomp_coeffs_4[i][j] * prod;
+		}
+	}
+
+	return cubic_solver_cheat(values) / (a*a + b*b);
+}
+
+std::array<std::complex<float>, 2> quadratic_solver(float a, float b, float c) {
+	if (std::abs(a) < eps) {
+		// todo
+		return std::array<std::complex<float>, 2> {0, 0};
+	}
+
+	auto D = b*b - 4*a*c;
+	if (D > eps) {
+		return std::array<std::complex<float>, 2> {(-b + std::sqrt(D)) / (2 * a),
+										(-b - std::sqrt(D)) / (2 * a)};
+	} else if (D < -eps) {
+		return std::array<std::complex<float>, 2> {(-b + If * std::sqrt(-D)) / (2 * a),
+										(-b - If * std::sqrt(-D)) / (2 * a)};
+	} else {
+		return std::array<std::complex<float>, 2> {(-b) / (2 * a), (-b) / (2 * a)};
+	}
+}
+
+// first root returned is real
+std::array<std::complex<float>, 3> cubic_solver(float d, float a, float b, float c, size_t &real_roots) {
+	// stolen from https://github.com/sasamil/Quartic/blob/master/quartic.cpp
+	if (std::abs(d) < eps) { 
+		auto roots = quadratic_solver(a, b, c);
+		return std::array<std::complex<float>, 3> { roots[0], roots[1], 0 };
+	}
+
+	a /= d; b /= d; c /= d;
+
+	auto a2 = a*a, q = (a2 - 3*b) / 9, r = (a*(2*a2 - 9*b) + 27*c) / 54,
+			r2 = r*r, q3 = q*q*q;
+
+	if (r2 < q3) {
+		auto t = r/std::sqrt(q3);
+		if (t < -1) { t = -1; }
+		else if (t > 1) { t = 1; }
+
+		t = std::acos(t); a /= 3; q = -2*std::sqrt(q);
+
+		real_roots = 3;
+		return std::array<std::complex<float>, 3> { q*std::cos(t/3) - a, 
+							q*std::cos((t + M_2PI) / 3) - a,
+							q*std::cos((t - M_2PI) / 3) - a };
+	} else {
+		float A = -std::pow(std::abs(r) + std::sqrt(r2 - q3), 1.0f/3);
+		float B;
+		if (r < 0) { A = -A; }
+		if (std::abs(A) < eps) { B = 0; } else { B = q/A; }
+
+		a /= 3;
+
+		float imaginary_part = 0.5 * std::sqrt(3.0) * (A - B);
+		float real_part = -0.5 * (A + B) - a;
+		if (std::abs(imaginary_part) < eps) {
+			real_roots = 3;
+			return std::array<std::complex<float>, 3> { (A + B) - a,
+									real_part, real_part };
+		} else {
+			real_roots = 1;
+			return std::array<std::complex<float>, 3> { (A + B) - a,
+									real_part + If * imaginary_part, real_part - If * imaginary_part };
+		}
+	}
+/*	auto p = -b / (3.0f * a), 
+		q = p*p*p + (b*c - 3.0f*a*d) / (6.0f*(a*a)),
+		r = c / (3.0f*a);
+
+	auto rp2 = (r - p*p), h = q*q + rp2*rp2*rp2;
+
+	auto h_root = std::sqrt(h + If * 0.0f);
+	auto x = std::pow(q + h_root, 1.0f/3.0f);
+	return std::array<std::complex<float>, 3> {x - rp2 / x + p, Omega*x - rp2 / (Omega*x) + p,
+												Omega_sq*x - rp2 / (Omega_sq*x) + p};*/
+}
+
+std::array<std::complex<float>, 4> quartic_solver(float e, float a, float b, float c, float d) {
+	// stolen from https://github.com/sasamil/Quartic/blob/master/quartic.cpp
+	if (std::abs(e) < eps) {
+		size_t trash;
+		auto roots = cubic_solver(a, b, c, d, trash);
+		return  std::array<std::complex<float>, 4> { roots[0], roots[1], roots[2], 0 };
+	}
+
+	a /= e; b /= e; c /= e; d /= e;
+
+	float a3 = -b;
+	float b3 =  a*c -4.*d;
+	float c3 = -a*a*d - c*c + 4.*b*d;
+
+	size_t iZeroes = 0;
+	auto x3 = cubic_solver(1, a3, b3, c3, iZeroes);
+
+	float q1, q2, p1, p2, D, sqD, y = x3[0].real();
+
+	if (iZeroes > 1) {
+		if (std::abs(x3[1].real()) > std::abs(y)) y = x3[1].real();
+		if (std::abs(x3[2].real()) > std::abs(y)) y = x3[2].real();
+	}
+
+	D = y*y - 4*d;
+	if (std::abs(D) < eps) {
+		q1 = q2 = y * 0.5;
+		D = a*a - 4*(b-y);
+		if(std::abs(D) < eps) {
+			p1 = p2 = a * 0.5;
+		} else {
+			sqD = std::sqrt(D);
+			p1 = (a + sqD) * 0.5;
+			p2 = (a - sqD) * 0.5;
+		}
+	} else {
+		sqD = std::sqrt(D);
+		q1 = (y + sqD) * 0.5;
+		q2 = (y - sqD) * 0.5;
+		p1 = (a*q1-c)/(q1-q2);
+		p2 = (c-a*q2)/(q1-q2);
+	}
+
+	auto first = quadratic_solver(1, p1, q1);
+	auto second = quadratic_solver(1, p2, q2);
+
+    return { first[0], first[1], second[0], second[1] };
+}
+
+float density_4(Eigen::MatrixXcf &A, Eigen::MatrixXcf &B, float a, float b, float precomp_coeffs_1[5][5],
+				float precomp_coeffs_2[5][5], float precomp_coeffs_3[5][5], float precomp_coeffs_4[5][5],
+				float precomp_coeffs_5[5][5]) {
+	float values[5] = { 0, 0, 0, 0, 0 };
+	float a_powers[5] = { 1, a, a*a, 0, 0 };
+	float b_powers[5] = { 1, b, b*b, 0, 0 };
+	for (size_t i = 3; i < 5; i++) {
+		a_powers[i] = a_powers[i - 1] * a;
+		b_powers[i] = b_powers[i - 1] * b;
+	}
+
+	for (size_t i = 0; i < 5; i += 1) {
+		for (size_t j = 0; j < 5; j += 1) {
+			float prod = a_powers[i] * b_powers[j];
+			values[0] += precomp_coeffs_1[i][j] * prod;
+			values[1] += precomp_coeffs_2[i][j] * prod;
+			values[2] += precomp_coeffs_3[i][j] * prod;
+			values[3] += precomp_coeffs_4[i][j] * prod;
+			values[4] += precomp_coeffs_5[i][j] * prod;
+		}
+	}
+
+	auto roots = quartic_solver(values[4], values[3], values[2], values[1], values[0]);
+	float value = 0;
+	for (const auto &root : roots) {
+		if (root.imag() > eps) {
+			value += root.imag() / ((a * root + b) * std::conj(a * root + b)).real();
+		}
+	}
+	return value;
+}
+
+void compute_output_thread(Eigen::MatrixXcf &A, Eigen::MatrixXcf &B, size_t idx) {
+	for (size_t i = (height / num_threads) * idx; i < (height / num_threads) * (idx + 1); i += 1) {
+		for (size_t j = 0; j < width; j += 1) {
+			output[i * width + j] = density_default(A, B, dr_scale * ((j * 2.0f / width) - 1.0f), 
+													-dr_scale * ((i * 2.0f / height) - 1.0f));
+		}
+	}
+}
+
+void compute_output_thread_2(Eigen::MatrixXcf &A, Eigen::MatrixXcf &B, float precomp_coeffs[7], size_t idx) {
+	for (size_t i = (height / num_threads) * idx; i < (height / num_threads) * (idx + 1); i += 1) {
+		for (size_t j = 0; j < width; j += 1) {
+			output[i * width + j] = density_2(A, B, dr_scale * ((j * 2.0f / width) - 1.0f), 
+													-dr_scale * ((i * 2.0f / height) - 1.0f), precomp_coeffs);
+		}
+	}
+}
+
+void compute_output_thread_3(Eigen::MatrixXcf &A, Eigen::MatrixXcf &B, float precomp_coeffs_1[7][7],
+								float precomp_coeffs_2[7][7], float precomp_coeffs_3[7][7], float precomp_coeffs_4[7][7],
+								size_t idx) {
+	for (size_t i = (height / num_threads) * idx; i < (height / num_threads) * (idx + 1); i += 1) {
+		for (size_t j = 0; j < width; j += 1) {
+			output[i * width + j] = density_3(A, B, dr_scale * ((j * 2.0f / width) - 1.0f), 
+													-dr_scale * ((i * 2.0f / height) - 1.0f), precomp_coeffs_1,
+													precomp_coeffs_2, precomp_coeffs_3, precomp_coeffs_4);
+		}
+	}
+}
+
+void compute_output_thread_4(Eigen::MatrixXcf &A, Eigen::MatrixXcf &B, float precomp_coeffs_1[5][5],
+								float precomp_coeffs_2[5][5], float precomp_coeffs_3[5][5], float precomp_coeffs_4[5][5],
+								float precomp_coeffs_5[5][5],
+								size_t idx) {
+	for (size_t i = (height / num_threads) * idx; i < (height / num_threads) * (idx + 1); i += 1) {
+		for (size_t j = 0; j < width; j += 1) {
+			output[i * width + j] = density_4(A, B, dr_scale * ((j * 2.0f / width) - 1.0f), 
+													-dr_scale * ((i * 2.0f / height) - 1.0f), precomp_coeffs_1,
+													precomp_coeffs_2, precomp_coeffs_3, precomp_coeffs_4, precomp_coeffs_5);
+		}
+	}
+}
+
+void compute_output(Eigen::MatrixXcf &C) {
+	Eigen::MatrixXcf A = (C + C.adjoint()) / 2.0f;
+	Eigen::MatrixXcf B = -If * (C - C.adjoint()) / 2.0f;
+
+	std::thread *threads = new std::thread[num_threads];
+	if (d == 2) {
+		float A_tr = A.trace().real();
+		float AA_tr = (A*A).trace().real();
+		float B_tr = B.trace().real();
+		float BB_tr = (B*B).trace().real();
+		float AB_tr = (A*B).trace().real();
+
+		float precomp_coeffs[9] = { -AB_tr*AB_tr + 2*A_tr*B_tr*AB_tr - A_tr*A_tr*BB_tr - B_tr*B_tr*AA_tr + AA_tr*BB_tr,
+			2*AA_tr*B_tr - 2*A_tr*AB_tr, A_tr*A_tr - 2*AA_tr, 
+			2*BB_tr*A_tr - 2*B_tr*AB_tr, 4*AB_tr - 2*A_tr*B_tr,
+			B_tr*B_tr - 2*BB_tr,
+			A_tr*A_tr - AA_tr, 2*AB_tr - 2*A_tr*B_tr, B_tr*B_tr - BB_tr };
+
+		for (size_t idx = 0; idx < num_threads; idx++) {
+			threads[idx] = std::thread(compute_output_thread_2, std::ref(A), std::ref(B), precomp_coeffs, idx);
+		}
+	} else if (d == 3) {
+		float tr[4][4];
+		Eigen::MatrixXcf A_powers[4] = { Eigen::MatrixXcf::Identity(A.rows(), A.cols()), A, A*A, Eigen::MatrixXcf::Identity(A.rows(), A.cols()) };
+		Eigen::MatrixXcf B_powers[4] = { Eigen::MatrixXcf::Identity(B.rows(), B.cols()), B, B*B, Eigen::MatrixXcf::Identity(A.rows(), A.cols()) };
+		A_powers[3] = A_powers[2] * A;
+		B_powers[3] = B_powers[2] * B;
+		for (size_t i = 0; i < 4; i += 1) {
+			for (size_t j = 0; j < 4; j += 1) {
+				tr[i][j] = (A_powers[i] * B_powers[j]).trace().real();
+			}
+		}
+
+		float precomp_coeffs_1[7][7] = {{0, 0, 0, -(tr[1][0]*tr[1][0]*tr[1][0]) + 3*tr[1][0]*tr[2][0] -
+							2*tr[3][0], 3*(tr[1][0]*tr[1][0]) - 3*tr[2][0], -6*tr[1][0], 6},
+							{0, 0, -3*tr[0][1]*(tr[1][0]*tr[1][0]) + 6*tr[1][0]*tr[1][1] +
+							3*tr[0][1]*tr[2][0] - 6*tr[2][1], 6*tr[0][1]*tr[1][0] - 6*tr[1][1],
+							-6*tr[0][1], 0, 0},
+							{0, -3*(tr[0][1]*tr[0][1])*tr[1][0] + 3*tr[0][2]*tr[1][0] +
+							6*tr[0][1]*tr[1][1] - 6*tr[1][2], 3*(tr[0][1]*tr[0][1]) - 3*tr[0][2] +
+							3*(tr[1][0]*tr[1][0]) - 3*tr[2][0], -12*tr[1][0], 18, 0, 0},
+							{-(tr[0][1]*tr[0][1]*tr[0][1]) + 3*tr[0][1]*tr[0][2] - 2*tr[0][3],
+							6*tr[0][1]*tr[1][0] - 6*tr[1][1], -12*tr[0][1], 0, 0, 0, 0},
+							{3*(tr[0][1]*tr[0][1]) - 3*tr[0][2], -6*tr[1][0], 18, 0, 0, 0, 0},
+							{-6*tr[0][1], 0, 0, 0, 0, 0, 0},
+							{6, 0, 0, 0, 0, 0, 0}};
+
+		float precomp_coeffs_2[7][7] = {{0, 0, 0, -3*tr[0][1]*tr[1][0]*tr[1][0] +
+							6*tr[1][0]*tr[1][1] + 3*tr[0][1]*tr[2][0] - 6*tr[2][1],
+							6*tr[0][1]*tr[1][0] - 6*tr[1][1], -6*tr[0][1], 0},
+							{0, 0, -6*tr[0][1]*tr[0][1]*tr[1][0] + 6*tr[0][2]*tr[1][0] + 3*tr[1][0]*tr[1][0]*tr[1][0] +
+							12*tr[0][1]*tr[1][1] - 12*tr[1][2] - 9*tr[1][0]*tr[2][0] + 6*tr[3][0],
+							6*tr[0][1]*tr[0][1] - 6*tr[0][2] - 6*tr[1][0]*tr[1][0] + 6*tr[2][0], 6*tr[1][0],
+							0, 0},
+							{0, -3*tr[0][1]*tr[0][1]*tr[0][1] + 9*tr[0][1]*tr[0][2] - 6*tr[0][3] +
+							6*tr[0][1]*tr[1][0]*tr[1][0] - 12*tr[1][0]*tr[1][1] - 6*tr[0][1]*tr[2][0] +
+							12*tr[2][1], 0, -12*tr[0][1], 0, 0, 0},
+							{3*tr[0][1]*tr[0][1]*tr[1][0] - 3*tr[0][2]*tr[1][0] - 6*tr[0][1]*tr[1][1] +
+							6*tr[1][2], 6*tr[0][1]*tr[0][1] - 6*tr[0][2] - 6*tr[1][0]*tr[1][0] + 6*tr[2][0],
+							12*tr[1][0], 0, 0, 0, 0},
+							{-6*tr[0][1]*tr[1][0] + 6*tr[1][1], -6*tr[0][1], 0, 0, 0, 0, 0},
+							{6*tr[1][0], 0, 0, 0, 0, 0, 0},
+							{0, 0, 0, 0, 0, 0, 0}};
+
+		float precomp_coeffs_3[7][7] = {{0, 0, 0, -3*tr[0][1]*tr[0][1]*tr[1][0] +
+							3*tr[0][2]*tr[1][0] + 6*tr[0][1]*tr[1][1] - 6*tr[1][2], 3*tr[0][1]*tr[0][1]
+							- 3*tr[0][2], 0, 0},
+							{0, 0, -3*tr[0][1]*tr[0][1]*tr[0][1] + 9*tr[0][1]*tr[0][2] - 6*tr[0][3] +
+							6*tr[0][1]*tr[1][0]*tr[1][0] - 12*tr[1][0]*tr[1][1] - 6*tr[0][1]*tr[2][0] +
+							12*tr[2][1], -6*tr[0][1]*tr[1][0] + 6*tr[1][1], 0, 0, 0},
+							{0, 6*tr[0][1]*tr[0][1]*tr[1][0] - 6*tr[0][2]*tr[1][0] - 3*tr[1][0]*tr[1][0]*tr[1][0] -
+							12*tr[0][1]*tr[1][1] + 12*tr[1][2] + 9*tr[1][0]*tr[2][0] - 6*tr[3][0],
+							3*tr[0][1]*tr[0][1] - 3*tr[0][2] + 3*tr[1][0]*tr[1][0] - 3*tr[2][0], 0, 0, 0, 0},
+							{-3*tr[0][1]*tr[1][0]*tr[1][0] + 6*tr[1][0]*tr[1][1] + 3*tr[0][1]*tr[2][0] -
+							6*tr[2][1], -6*tr[0][1]*tr[1][0] + 6*tr[1][1], 0, 0, 0, 0, 0},
+							{3*tr[1][0]*tr[1][0] - 3*tr[2][0], 0, 0, 0, 0, 0, 0},
+							{0, 0, 0, 0, 0, 0, 0},
+							{0, 0, 0, 0, 0, 0, 0}};
+
+		float precomp_coeffs_4[7][7] = {{0, 0, 0, -(tr[0][1]*tr[0][1]*tr[0][1]) + 3*tr[0][1]*tr[0][2] -
+							2*tr[0][3], 0, 0, 0},
+							{0, 0, 3*tr[0][1]*tr[0][1]*tr[1][0] - 3*tr[0][2]*tr[1][0] -
+							6*tr[0][1]*tr[1][1] + 6*tr[1][2], 0, 0, 0, 0},
+							{0, -3*tr[0][1]*tr[1][0]*tr[1][0] + 6*tr[1][0]*tr[1][1] +
+							3*tr[0][1]*tr[2][0] - 6*tr[2][1], 0, 0, 0, 0, 0},
+							{tr[1][0]*tr[1][0]*tr[1][0] - 3*tr[1][0]*tr[2][0] + 2*tr[3][0], 0, 0, 0, 0, 0, 0},
+							{0, 0, 0, 0, 0, 0, 0},
+							{0, 0, 0, 0, 0, 0, 0},
+							{0, 0, 0, 0, 0, 0, 0}};
+							
+		for (size_t idx = 0; idx < num_threads; idx++) {
+			threads[idx] = std::thread(compute_output_thread_3, std::ref(A), std::ref(B), precomp_coeffs_1,
+										precomp_coeffs_2, precomp_coeffs_3, precomp_coeffs_4, idx);
+		}
+	} else if (d == 4) {
+		float tt0 = A.trace().real();
+		float tt1 = B.trace().real();
+		float tt00 = (A*A).trace().real();
+		float tt01 = (A*B).trace().real();
+		float tt11 = (B*B).trace().real();
+		float tt000 = (A*A*A).trace().real();
+		float tt001 = (A*A*B).trace().real();
+		float tt011 = (A*B*B).trace().real();
+		float tt111 = (B*B*B).trace().real();
+		float tt0000 = (A*A*A*A).trace().real();
+		float tt0001 = (A*A*A*B).trace().real();
+		float tt0011 = (A*A*B*B).trace().real();
+		float tt0101 = (A*B*A*B).trace().real();
+		float tt0111 = (A*B*B*B).trace().real();
+		float tt1111 = (B*B*B*B).trace().real();
+
+		float precomp0[5][5] = {{tt1*tt1*tt1*tt1 - 6*tt1*tt1*tt11 + 3*tt11*tt11 + 8*tt1*tt111 - 6*tt1111,-4*tt1*tt1*tt1 + 12*tt1*tt11 - 8*tt111,12*tt1*tt1 - 12*tt11,-24*tt1,24},
+		{0,0,0,0,0},
+		{0,0,0,0,0},
+		{0,0,0,0,0},
+		{0,0,0,0,0}};
+
+		float precomp1[5][5] = {{-24*tt0111 + 24*tt011*tt1 - 12*tt01*tt1*tt1 + 4*tt0*tt1*tt1*tt1 + 12*tt01*tt11 - 12*tt0*tt1*tt11 + 8*tt0*tt111,- 24*tt011 + 24*tt01*tt1 - 12*tt0*tt1*tt1 + 12*tt0*tt11,-24*tt01 + 24*tt0*tt1,-24*tt0,0},{-4*tt1*tt1*tt1 + 12*tt1*tt11 - 8*tt111,24*tt1*tt1 - 24*tt11,-72*tt1,96,0},
+		{0,0,0,0,0},
+		{0,0,0,0,0},
+		{0,0,0,0,0}};
+
+		float precomp2[5][5] = {{-24*tt0011 + 12*tt01*tt01 - 12*tt0101 + 24*tt0*tt011 + 24*tt001*tt1 - 24*tt0*tt01*tt1 + 6*tt0*tt0*tt1*tt1 - 6*tt00*tt1*tt1 - 6*tt0*tt0*tt11 + 6*tt00*tt11,-24*tt001 + 24*tt0*tt01 - 12*tt0*tt0*tt1 + 12*tt00*tt1,12*tt0*tt0 - 12*tt00,0,0},
+		{- 24*tt011 + 24*tt01*tt1 - 12*tt0*tt1*tt1 + 12*tt0*tt11,-48*tt01 + 48*tt0*tt1,-72*tt0,0,0},
+		{12*tt1*tt1 - 12*tt11,-72*tt1,144,0,0},
+		{0,0,0,0,0},
+		{0,0,0,0,0}};
+
+		float precomp3[5][5] = {{-24*tt0001 + 24*tt0*tt001 - 12*tt0*tt0*tt01 + 12*tt00*tt01 + 4*tt0*tt0*tt0*tt1 - 12*tt0*tt00*tt1 + 8*tt000*tt1,-4*tt0*tt0*tt0 + 12*tt0*tt00 - 8*tt000,0,0,0},
+		{-24*tt001 + 24*tt0*tt01 - 12*tt0*tt0*tt1 + 12*tt00*tt1,24*tt0*tt0 - 24*tt00,0,0,0},{-24*tt01 + 24*tt0*tt1,-72*tt0,0,0,0},
+		{-24*tt1,96,0,0,0},
+		{0,0,0,0,0}};
+
+		float precomp4[5][5] = {{tt0*tt0*tt0*tt0 - 6*tt0*tt0*tt00 + 3*tt00*tt00 + 8*tt0*tt000 - 6*tt0000,0,0,0,0},
+		{-4*tt0*tt0*tt0 + 12*tt0*tt00 - 8*tt000,0,0,0,0},{12*tt0*tt0 - 12*tt00,0,0,0,0},
+		{-24*tt0,0,0,0,0},
+		{24,0,0,0,0}};
+
+		for (size_t idx = 0; idx < num_threads; idx++) {
+			threads[idx] = std::thread(compute_output_thread_4, std::ref(A), std::ref(B), precomp0,
+										precomp1, precomp2, precomp3, precomp4, idx);
+		}
+	} else {
+		for (size_t idx = 0; idx < num_threads; idx++) {
+			threads[idx] = std::thread(compute_output_thread, std::ref(A), std::ref(B), idx);
+		}
+	}
+
+	for (size_t idx = 0; idx < num_threads; idx++) {
+		threads[idx].join();
+	}
+
+	for (size_t i = 0; i < height; i += 1) {
+		for (size_t j = 0; j < width; j += 1) {
+			float comp = 40.0f * output[i * width + j] / (1 + 40.0f * output[i * width + j]);
+			uint32_t comp_255 = (uint32_t)((uint8_t)(comp * 255.0f));
+			output_buffer[i * width + j] = (0x000000FF |
+				((255 - comp_255) << 8) | ((255 - comp_255) << 16) | (comp_255 << 24)) ^ (comp_255 >> 1);
+		}
+	}
+
+	Eigen::ComplexEigenSolver<Eigen::MatrixXcf> eigensolver(A.inverse() * B);
+	if (eigensolver.info() != Eigen::Success) { return; }
+
+	auto eigenvectors  = eigensolver.eigenvectors();
+	for(size_t i = 0; i < eigenvectors.cols(); i++) {
+		auto v = eigenvectors.col(i);
+		auto norm = (v.dot(v)).real();
+		EM_ASM({
+			add_tangent($0, $1);
+		}, (v.dot(A * v)).real() / norm, (v.dot(B * v)).real() / norm);
+	}
+}
+
+void EMSCRIPTEN_KEEPALIVE set_zoom(float zoom) {
+	zoom_scale = zoom;
+}
+
+void EMSCRIPTEN_KEEPALIVE set_matrix(size_t dim, size_t pointer) {
+	d = dim;
+	dr_scale = zoom_scale * dim * p_scale;
+
+	float *entries = (float*)pointer;
+	Eigen::MatrixXcf C = Eigen::MatrixXcf::Zero(d, d);
+	for (size_t i = 0; i < d; i++) {
+		for (size_t j = 0; j < d; j++) {
+			C(i, j) = entries[2 * (i * d + j)] + entries[2 * (i * d + j) + 1] * If;
+		}
+	}
+
+	compute_output(C);
+}
+
+void set_resolution(size_t resolution) {
+	delete output;
+	delete output_buffer;
+
+	width = height = resolution;
+	output = new float[width * height];
+	output_buffer = new uint32_t[width * height];
+	std::memset(output, 0, width * height * sizeof(float));
+}
+
+int main() {
+	val canvas = document.call<val>("getElementById", val("drawing_canvas"));
+	val ctx = canvas.call<val>("getContext", val("2d"));
+	ctx.set("fillStyle", "red");
+
+	width = EM_ASM_INT({ return document.getElementById("drawing_canvas").width; });
+	height = EM_ASM_INT({ return document.getElementById("drawing_canvas").height; });
+	num_threads = EM_ASM_INT({ return navigator.hardwareConcurrency; });
+}
+
+EMSCRIPTEN_BINDINGS(my_module) {
+	emscripten::function("get_output_buffer", &get_output_buffer);
+	emscripten::function("set_matrix", &set_matrix);
+	emscripten::function("set_resolution", &set_resolution);
+	emscripten::function("set_zoom", &set_zoom);
+}