Interactive visualization of polynomial roots in the complex plane, inspired by 2swap's video on the quintic formula.
This is a C application I wrote after viewing 2swap's video explaining the non-existence of the quintic formula. The visuals in the video are nothing short of incredible, being mesmerising and beautiful. Compelled to create my own version to play around with the polynomial visualisations myself, this project came into existence.
This project recreates two of the different visualisations shown in the video. The first plots the behavior of the polynomial given inputs across the complex plane, and allows the manipulation of roots and coefficients to observe the effect one has on another. The second mode renders the roots of all the combinations of
roots mode: shows a single polynomial's output (argument → hue, magnitude → brightness) given inputs on the complex plane. drag roots or coefficients around to see how things change.

cloud mode: generates all 
Standard CMake project. Tested to build on macOS and Windows (using the MSVC compiler). (It is because I tried to support MSVC that the complex number operations look so ugly—MSVC does not support the typical complex keyword and intrinsic complex arithmetic).
This project relies on SDL3 and cJSON (fetched automatically in CMake).
- OpenCL for gpu rendering
- LAPACK (or Intel MKL on Windows) for faster root finding
- OpenMP for multithreading (highly recommended).
mkdir build && cd build
cmake .. -DCMAKE_BUILD_TYPE=Release
cmake --build . --config Release| Option | Default | Description |
|---|---|---|
STATIC_LINK |
ON | Static link runtime libraries |
PORTABLE_BUILD |
OFF | Disable -march=native for portable binaries |
BUILD_BENCHMARK |
ON | Build benchmark executables |
BUILD_TESTS |
ON | Build test suite |
Pretty self explanatory since it is shown in the UI.
- drag points to move roots/coefficients
- scroll to zoom, right-click drag to pan
y/h- increase/decrease degree+/-- add/remove base coefficients (cloud mode)t- toggle between roots and cloud moden/m- adjust gamma (roots) or point radius (cloud)r- reset viewF1- toggle info overlay
and a few others that I do not remember
(yes, this supports animations)
pass a json file as argument to load an animation script:
./polynomial animation.jsonspace to play/pause, x to reset.
Animation scripts are JSON files with keyframe-based interpolation. Example structure:
{
"mode": "roots",
"degree": 3,
"loop": true,
"loop_delay": 0.5,
"keyframes": [
{
"time": 0.0,
"points": [[1.0, 0.0], [-1.0, 0.0], [0.0, 1.0]],
"ease": "smooth"
},
{
"time": 2.0,
"points": [[-1.0, 0.0], [1.0, 0.0], [0.0, -1.0]],
"ease": "inout"
}
]
}Each point is [real, imaginary]. The mode field determines how points are interpreted:
| Mode | Description |
|---|---|
roots |
Points animate polynomial roots directly. Coefficients are computed from roots. |
coeffs |
Points animate polynomial coefficients. Roots are computed via solver. |
cloud |
Points animate base coefficients used to generate all polynomial combinations. |
The degree field explicitly sets the polynomial degree. If omitted, degree is inferred from the number of points in keyframes (for roots mode: degree = num_points; for coeffs mode: degree = num_points - 1). Explicit degree is useful when animating fewer points than the polynomial degree.
| Field | Type | Default | Description |
|---|---|---|---|
mode |
string | "cloud" |
"roots", "coeffs", or "cloud" |
degree |
int | inferred | Polynomial degree (explicit override) |
num_coeffs |
int | inferred | Number of base coefficients (cloud mode) |
loop |
bool | false |
Restart animation after completion |
loop_delay |
float | 0.0 |
Seconds to pause before looping |
keyframes |
array | required | Array of keyframe objects |
Polynomials are stored with both coefficient and root representations. The polynomial_t struct maintains:
- Coefficients ordered as
$k_n x^n + k_{n-1} x^{n-1} + \cdots + k_1 x + k_0$ - Distinct roots with their multiplicities, in the case of repeated roots
When manipulating a root, the coefficients are recomputed by expanding from the factored form. When manipulating a coefficient, the roots are found via the solver. The colours on the screen are rendered using OpenCL if available, or CPU fallback if not available.
The solver implements a scaled variant of the Jenkins-Traub algorithm, but when available the LAPACK eigenvalue finding function zgeev is preferred. This is because any iterative method for root finding where each root is 'deflated' (divided) out of the polynomial to find the next root invariably will accumulate numerical errors, especially due to numerical precision issues with very large or small roots or coefficients. The companion matrix eigenvalue method avoids this by solving for all the roots together.
For a monic polynomial
has eigenvalues equal to the polynomial's roots.
Given
It would be prohibitively computationally expensive to perform root-finding on tens of thousands of polynomials every frame during drag operations. Several optimizations make this tractable:
The observation is that the roots have to morph continuously. When a base coefficient moves slightly, most roots also move slightly. Instead of full re-solving every frame, we can:
- Estimates new root positions using first-order approximation (implicit differentiation of
$P(r) = 0$ ) - Uses only Stage 3 of Jenkins-Traub to refine from these estimates
- Falls back to full solve if refinement fails to converge
This exploits temporal coherence—during a drag, we have excellent initial guesses from the previous frame.
Even with incremental solving, numerical drift accumulates. The solver maintains a since_last_update counter per combination, and the initial full-solves are staggered. By doing so we spread out the full root-finding such that only a small subset of all the polynomials get a full-solve in any one frame. It is known that when the coefficients near zero, the incremental solve will have a higher rate of failure and will fallback to a full solve, which is a potential source of lag.
When OpenCL is available, rendering runs entirely on the GPU via four kernels:
-
evaluate_polynomial— Evaluates$P(z)$ at each pixel using Horner's method. Each work item handles one pixel. -
compute_max_magnitude/reduce_max_final— Two-stage parallel reduction to find the global maximum magnitude for normalization. -
render_frame— Color mapping kernel. Reads the evaluation results, computes HSL from magnitude/argument, applies gamma correction, and outputs RGBA pixels.
For cloud mode, these kernels handle the point rendering:
-
render_point_cloud— Each work item handles one polynomial's roots. Uses atomic operations to accumulate hue vectors (as cos/sin components) into a per-pixel buffer to average hues when multiple points overlap. -
render_point_cloud_hue— Converts accumulated hue vectors to final pixel colors and outputs to a RGBA buffer.
When OpenCL is not available, rendering uses OpenMP-parallelized CPU evaluation.
The code jumps through some hoops for MSVC compatibility (seriously why...):
- No
_Complexkeyword—uses explicit struct-based complex numbers with function wrappers (cxadd,cxmul, etc.) - Long double precision is optional (
COMPAT_COMPLEX_DOUBLE_ONLY)—MSVC'slong doubleis justdoubleanyway