diff --git a/public/main.js b/public/main.js
new file mode 100644
index 0000000..3ccce37
--- /dev/null
+++ b/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
+ }
+}
diff --git a/public/tjsm.html b/public/tjsm.html
new file mode 100644
index 0000000..1cded98
--- /dev/null
+++ b/public/tjsm.html
@@ -0,0 +1,60 @@
+
+
+
+
+
+
+tjsm
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/src/main.cpp b/src/main.cpp
new file mode 100644
index 0000000..fe2378e
--- /dev/null
+++ b/src/main.cpp
@@ -0,0 +1,548 @@
+// CC BY-NC-SA 4.0, gtrrebel, shikhin, zgrep
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+const std::complex If(0.0f, 1.0f);
+
+const std::complex Omega(-1.0f/2.0f, std::sqrt(3.0f) / 2);
+const std::complex 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, 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 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, 2> quadratic_solver(float a, float b, float c) {
+ if (std::abs(a) < eps) {
+ // todo
+ return std::array, 2> {0, 0};
+ }
+
+ auto D = b*b - 4*a*c;
+ if (D > eps) {
+ return std::array, 2> {(-b + std::sqrt(D)) / (2 * a),
+ (-b - std::sqrt(D)) / (2 * a)};
+ } else if (D < -eps) {
+ return std::array, 2> {(-b + If * std::sqrt(-D)) / (2 * a),
+ (-b - If * std::sqrt(-D)) / (2 * a)};
+ } else {
+ return std::array, 2> {(-b) / (2 * a), (-b) / (2 * a)};
+ }
+}
+
+// first root returned is real
+std::array, 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, 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, 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, 3> { (A + B) - a,
+ real_part, real_part };
+ } else {
+ real_roots = 1;
+ return std::array, 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, 3> {x - rp2 / x + p, Omega*x - rp2 / (Omega*x) + p,
+ Omega_sq*x - rp2 / (Omega_sq*x) + p};*/
+}
+
+std::array, 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, 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 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("getElementById", val("drawing_canvas"));
+ val ctx = canvas.call("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);
+}