From 417b46e1ac7c1795cd12a0d57114fc83569be4da Mon Sep 17 00:00:00 2001 From: shikhin Date: Tue, 5 Dec 2023 15:00:17 -0500 Subject: [PATCH] init --- public/main.js | 544 ++++++++++++++++++++++++++++++++++++++++++++++ public/tjsm.html | 60 ++++++ src/main.cpp | 548 +++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 1152 insertions(+) create mode 100644 public/main.js create mode 100644 public/tjsm.html create mode 100644 src/main.cpp 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); +}