Update documentation and perform benchmarking

Author: Sekqies

README.md

@@ -113,3 +113,6 @@ An exhaustive list of features supported by this tool
 
 ### [Advanced](docs/advanced.md)
 An architectural explanation of how this tool works, internally.
+
+### [Optimizations](docs/optimization.md)
+A list of optimizations used in this project's developments

docs/advanced.md

@@ -95,9 +95,18 @@ $$
 
 Where $|l(z)| = \frac{2}{\pi}\arctan(|z|)$
 
+Currently, we also use a gamma factor, $k$, to make the light ascent from black to white slower, as per request of one of the users. The new definition is simply
+
+$|l(z)| = \frac{2}{\pi}\arctan(|z|^k)$
+
 ## 3D mode
 3D mode is rendered by making an NxN grid mesh, handled by the `Mesh` struct and `create_grid_mesh` function, and shaping it around the function's "real" shape (picking points in set spaces and placing vertices there). 
 Height is directly mapped to |f(z)|, so information is redundant. Movement is done through a custom `Camera` struct.
 
 ## Picker
-Whenever you hover over a value and get a number back, this is done through a "picker" shader. It's a 1x1 grid rendered off-screen that calculates the value of f(z) precisely where you are hovering. The result is then shared as an `out vec4`, the first two floats being `z`, and second two, `f(z)`.
+Whenever you hover over a value and get a number back, this is done through a "picker" shader. It's a 1x1 grid rendered off-screen that calculates the value of f(z) precisely where you are hovering. The result is then shared as an `out vec4`, the first two floats being `z`, and second two, `f(z)`.
+
+## Ultra-High Precision Mode
+Ultra-High Precision Mode, or Arbitrary Precision Mode, works by running this program's shaders in the CPU. They are first transpiled from GLSL to C++ in the build step, done in `math_transpiler.py`, and a map from a GLSL function definition to a C++ function definition is set. 
+At runtime, the function stack is read and evaluated. This is done through multiple threads, that concurrently update their rows. This allows for the user to see the function as it's being updated.
+There exists a custom made math library, but the final product uses Boost::Multiprecision for performance. It is set for 50 digits of precision.

docs/assets/high_precision_mandelbrot.png

@@ -65,3 +65,7 @@ Time is a supported variable, `t`. Modulo and trigonometric operators are suppor
 - `(1-sin(t))*z + sin(t)*z^2`
 - `tan(z)^(sin(t)) % sin(z)`
 ![Animation of tan(z)^(sin(t)) % sin(z)](assets/gifs/Animation.gif)
+
+## Ultra high precision mode
+Plots can be rendered in arbitrary precision through the "Arbitrary Precision" subheader. This means that functions will be rendered up to a user-defined value of decimal digits, which defaults to 50. This is particularly useful for seeing fractal-like functions in high zoom values
+![Render of mandelbrot set with UHPM versus standard mode](assets/high_precision_mandelbrot.png)

docs/features.md

@@ -0,0 +1,88 @@
+# Optimizations
+There are many little optimizations done all throughout this project, that makes it able to run in very fast speeds. 
+
+## Benchmarking
+First, it is important to show that the optimizations used in this tool actually provide a substantial performance increase. For such, 4 distinct publicly available complex plotter were used for comparison:
+- Samuel Li's [Complex Function Plotter](https://samuelj.li/complex-function-plotter/). The main inspiration for this project.
+- Peter E. Francis' [Complex Function Plot](https://peterefrancis.com/complex-function-plot/plotter.html). A CPU based plotter
+- David Bau's [Conformal Map Plotter](https://davidbau.com/conformal/). While this function does not use regular domain coloring, it _still_ is a complex function plotter
+- Fernando Theodoro & Mateus Bastazini's [Complex Functions](https://www2.fc.unesp.br/matematicaearte/plotter/). A project I contributed to previously
+
+The computer for benchmarking had the following specs:
+CPU: Intel Core Ultra 7 255HX
+Random Memory: 16GB DDR5 5600 MT/s
+Integrated Graphics Card: NVIDIA Geforce RTX 5060
+
+All benchmarking measures the _parsing and drawing_ time for different complex functions. 
+1. Polynomial: $z^{10} - z^9 - z^8 - z^7 - z^6 - z^5 - z^4 - z^3 - z^2 - z - 1$
+2. Trigonometric and exponential: $\sin(\cos(\tan(z))) \cdot e^z$
+3. Non-elementary: $\Gamma(\zeta(z))$
+
+| Creator | Uses GPU | Polynomial | Trigonometric and exponential | Non-elementary |
+| :--- | :---: | :--- | :--- | :--- |
+| **Peter E. Francis** | No | 11336 ms | 7098 ms | *Not supported* |
+| **Samuel Li** | Yes | 279 ms | 459 ms | 11369 ms |
+| **David Bau** | Yes | 369 ms | 120 ms | *Not supported* |
+| **Theodoro & Bastazini** | Yes | 226 ms | 293 ms | 194 ms |
+| **This (Web)** | Yes | 4.14 ms | **2.42 ms** | **2.27 ms** |
+| **This (Desktop)**| Yes | **2.91 ms** | 3.93 ms | 3.14 ms |
+
+It is important to notice that the web version outperforms the desktop version for small functions. This is likely due to optimizations introduced in the Emscripten to WebAssembly optimizations, besides the transpilation from OpenGL to WebGL. 
+
+It can be observed from the table that this tool outperforms all of its other competitors. This is due to our optimizations!
+
+## Architectural Optimizations
+
+### Compiled and Interpreted Modes
+One of the main bottlenecks of GPU based renderers is dynamically writing the shader code, and compiling it at runtime, which introduces a short stutter. This tool uses **per-request** compiling, meaning that it only compiles the shader if the user splicitly requests it.
+Otherwise, the parsed expression is transformed into a sequence of bytecode instructions in Reverse-Polish Notation and sent to the shader as a texture. The shader then evaluates this as a stack. This requires no recompilations. 
+For higher performance, the expression can also be turned into a GLSL string, and recompiled.
+
+### Web Assembly and C++
+This tool is mainly desktop-focused, but the web version is compiled directly from C++ into WebAssembly using Emscripten. Web Assembly is a far lower-level, higher-performance alternative to Javascript, which makes its use a big player in outperforming other javascript-based plotters.
+
+### GPU Rendering
+All simple calculations are done in the GPU, rather than the CPU. This allows the plots to be drawn far faster than their CPU counterparts (as evidenced in the Benchmarking section)
+
+### Threading
+For the high precision CPU renders, the plot is drawn in multiple concurrent threads. 
+
+## Numeric optimizations
+### Constant Folding
+Constant expressions are evaluated before the drawing and rendering logic. For instance, `(13 * 7)*z` will be sent for rendering as `91 * z`
+
+### Simplification
+Expressions that can be evaluated to constants will be, during runtime. Such examples include:
+1. Division by itself `z/z = 1` (with an added singularity at `z=0`)
+2. Composition of inverses `sin(arcsin(z)) = z`
+
+### Analyic Differentiation
+Derivatives are calculated analytically rather than numerically. This means that an expression such as `d/dz(z^2 + z)` is correctly transformed into `2z + 1` in parse-time, rather than numerically evaluated through approximations in render-time.
+
+## Memory Optimizations
+
+### Manual Register Management and Mutable Math
+In the now deprecated arbitrary-precision math library for the GPU, registers were manually managed. Rather than using local variables, 16 registers of shared memory was used between every complex function. This was transpiled automatically from the low precision GLSL. 
+Also, this library uses non-functional math. This means that rather than functions returning a new copy to their result, they modify a given register. Such an example is the hp_add function:
+```glsl
+number R[16];
+number hp_add(number a, number b){
+    void hp_add(in number a, in number b, out number res) {
+    if (a.sign == b.sign) {
+        abs_sum(a, b, R[0]);
+        res.sign = a.sign;
+        return;
+    }
+    int cmp = 0;
+    compare_abs(a, b, cmp);
+    if (cmp >= 0) {
+        abs_hp_sub(a, b, R[0]);
+        res.sign = a.sign;
+        return;
+    }
+    abs_hp_sub(b, a, R[0]);
+    R[0].sign = b.sign;
+    res = R[0];
+}
+}
+```

docs/optimization.md

@@ -22,6 +22,7 @@
 #include <preprocessor/transpiler.h>
 #include <cpu_drawing/cpu_render.h>
 #include <glsl_generated/generated_math_mapper.h>
+#include <chrono>
 
 
 #include <preprocessor/string_builder.h>
@@ -140,7 +141,6 @@ void export_to_png(AppContext* ctx, int target_width, int target_height, const c
 	if (ctx->function_state->is_3d) glEnable(GL_DEPTH_TEST);
 	else glDisable(GL_DEPTH_TEST);
 
-	draw_scene(ctx, (float)target_width, (float)target_height);
 
 	glPixelStorei(GL_PACK_ALIGNMENT, 1);
 	unsigned char* pixels = new unsigned char[target_width * target_height * 4];
@@ -199,7 +199,6 @@ void main_loop_step(AppContext* ctx) {
 
 	ctx->function_state->is_3d = ctx->view_state->is_3d;
 	init_imgui_loop();
-
 	if(ctx->view_state->wants_high_precision){
 		const int hp_width = ctx->view_state->hp_width;
 		const int hp_height = ctx->view_state->hp_height;
@@ -245,6 +244,7 @@ void main_loop_step(AppContext* ctx) {
 		render_and_update(*(ctx->function_state), *(ctx->view_state), stack_tbo_texture, constants_tbo_texture, *(ctx->shader_program), *(ctx->compiled_shader));
 	}
 
+	
 	draw_scene(ctx, ctx->view_state->width, ctx->view_state->height);
 
 	if (is_3d != ctx->view_state->is_3d) {
@@ -338,10 +338,10 @@ int main() {
 	const string frag_source_3d = get_source("shaders/plotter3d.frag");
 	try {
 		vert_source_3d = build_shader_string(vert_source_3d, frag_source);
-		std::cout << vert_source_3d;
+		//std::cout << vert_source_3d;
 	}
 	catch (std::runtime_error& er) {
-		std::cout << er.what();
+		//std::cout << er.what();
 	}
 	shader_3d.compile(vert_source_3d, frag_source_3d);