InauguralSystems · InauguralPhysicist · Aug 27, 2026 · Aug 27, 2026
@@ -4,6 +4,44 @@ All notable changes to EigenScript are documented here.
 
 ## [Unreleased]
 
+### Added
+
+- **`lib/complex.eigs` — complex arithmetic, polar form and polynomial roots
+  (#1043).** Mode analysis, signal processing and root finding are all
+  intrinsically complex-valued, and every consumer was re-rolling the same
+  arithmetic on `[re, im]` pairs — including the stdlib itself, where
+  `engineering.magnitude_spectrum` and `power_spectrum` were hand-rolled
+  moduli. Adds `add`/`sub`/`mul`/`div`/`neg`/`conj`/`scale`/`mag`/`mag2`/
+  `arg`, `to_polar`/`from_polar`, `eq_near`, and `poly_eval`/`poly_roots`
+  (Durand-Kerner). `div` returns `null` on a zero denominator rather than
+  raising, which is what lets `poly_roots` hold an iterate at a repeated root
+  instead of dying on it.
+
+  Its tests are DIRECT unit checks against hand-computed values, deliberately:
+  a differential oracle is structurally blind here — phugoid measured a halved
+  `div` surviving a full root-finding differential, because Durand-Kerner
+  self-corrects under a scaled delta. That exact fault is planted in
+  `tests/test_complex.eigs`.
+
+- **`linalg.charpoly` and `linalg.eigenvalues` — general eigenvalues (#1042).**
+  `linalg` stopped at `eigenvalues_2x2`, which returns `null` for any complex
+  spectrum (a rotation matrix's ±i read as "no answer"), and `numerics` offered
+  only the dominant real eigenvalue. `charpoly` is Faddeev-LeVerrier for any
+  square n; `eigenvalues` returns all n as complex pairs. Graded on a 4x4 with
+  a known complex spectrum by the identities that tie the answer back to the
+  matrix rather than the solver — sum of eigenvalues = `mat_trace`, product =
+  `mat_det` — plus a residual check per eigenvalue.
+
+- **`engineering.phase_spectrum` (#1043).** The missing sibling of
+  `magnitude_spectrum`/`power_spectrum`, which now delegate to `complex.mag`
+  and `complex.mag2` (verified byte-identical to the hand-rolled versions
+  across four signals). A spectrum has a phase as well as a magnitude; it was
+  absent because the arithmetic to express it did not exist.
+
+  `linalg` and `engineering` are the first stdlib modules to import another
+  (`complex`). Every `lib/*.eigs` is installed together, so the dependency
+  travels with them.
+
 ## [0.41.0] - 2026-08-23
 
 ### Added

@@ -255,7 +255,7 @@ The minimal build (`make build`) sets all flags to 0. The full build
 
 ## Standard Library
 
-The 76 modules in `lib/` are pure EigenScript — no C code. They are loaded at
+The 77 modules in `lib/` are pure EigenScript — no C code. They are loaded at
 runtime via `load_file of "lib/module.eigs"`. Path resolution searches in
 order: the current working directory, the script file's directory, the script's
 parent directory, directories relative to the executable (`exe_dir/..` and the

@@ -22,7 +22,8 @@ mode is re-inventing what already ships (see "Before you hand-roll" below).
 | average / median / stddev / quantiles / correlation | `stats.mean`, `stats.median`, `stats.std_dev`, `stats.quantile`, `stats.correlation`, `stats.describe` | `lib/stats.eigs` |
 | tables, CSV, group-by, join, select/where | `data.df_from_csv`, `df_where`, `df_group_by`, `df_join`, `df_select` | `lib/data.eigs` |
 | probability distributions, combinatorics, Bayes | `probability.normal_pdf`, `binomial_pmf`, `combinations`, `bayes` | `lib/probability.eigs` |
-| matrices / vectors / solve Ax=b / eigenvalues | `linalg.mat_mul`, `mat_inverse`, `solve_linear`, `eigenvalues_2x2` | `lib/linalg.eigs` |
+| matrices / vectors / solve Ax=b / eigenvalues | `linalg.mat_mul`, `mat_inverse`, `solve_linear`, `charpoly`, `eigenvalues` | `lib/linalg.eigs` |
+| complex arithmetic, polar form, polynomial roots | `complex.mul`, `complex.div`, `complex.to_polar`, `complex.poly_roots` | `lib/complex.eigs` |
 | element-wise / GPU-style tensor math (matmul, softmax) | `matmul`, `softmax`, `add`, `mean` (**builtins**) | BUILTINS.md → Tensor Math |
 | derivatives / integrals / root-finding / ODEs | `calculus.derivative`, `integrate_simpson`, `newton_raphson`, `rk4_method` | `lib/calculus.eigs` |
 | minimize a function (gradient descent, annealing, GA) | `optimize.gradient_descent`, `simulated_annealing`, `genetic_optimize` | `lib/optimize.eigs` |
@@ -1001,7 +1002,7 @@ functions with **no `# name of args -> type` comment** (or only a section
 banner) above them — signatures in the tables above were reconstructed from
 the `define name(params)` line and should be back-filled in the lib files:
 `bcd`, `checksum`, `format` (`hexdump`), `datetime` (civil-math half),
-`eigen`, `harness`, `observer_slots`, `store`, `queue`, `lab`, `linalg`
+`eigen`, `harness`, `observer_slots`, `store`, `queue`, `lab`, `linalg`, `complex`
 (`mat_inverse`), `engineering` (unit converters), `geometry` (vector/solid
 helpers), and the `ui`/`ui_w_*`/`ui_theme`/`ui_anim`/`ui_draw` widget modules
 (which document widgets in the module header instead). Adding the per-function
@@ -1148,6 +1149,31 @@ Pure-EigenScript matrices (lists of lists) and vectors: transpose/multiply/deter
 | `least_squares` | `least_squares of [A, b]` | Solve overdetermined Ax ~ b via normal equations |
 | `eigenvalues_2x2` | `eigenvalues_2x2 of A` | eigenvalues of 2x2 matrix via characteristic polynomial |
 | `eigenvectors_2x2` | `eigenvectors_2x2 of A` | eigenvectors for each eigenvalue of 2x2 matrix |
+| `charpoly` | `charpoly of A` | characteristic polynomial coefficients of any square A (Faddeev-LeVerrier); monic, leading 1 implied |
+| `eigenvalues` | `eigenvalues of A` | all n eigenvalues of any square A as complex `[re, im]` pairs — handles conjugate pairs, which `eigenvalues_2x2` refuses |
+
+### lib/complex.eigs
+
+Complex numbers as two-element `[re, im]` lists — the shape
+`engineering.dft` already returns.
+
+| Function | Signature | Notes |
+|---|---|---|
+| `add` / `sub` / `mul` | `mul of [a, b]` | complex arithmetic |
+| `div` | `div of [a, b]` | returns `null` when b is zero (division by zero raises in EigenScript) |
+| `neg` / `conj` | `conj of a` | negation flips both components, conjugation only the imaginary one |
+| `scale` | `scale of [a, s]` | multiply by a REAL scalar |
+| `mag` / `mag2` | `mag of a` | modulus, and its square without the `sqrt` |
+| `arg` | `arg of a` | argument in radians, `(-pi, pi]`, via `atan2` |
+| `to_polar` / `from_polar` | `to_polar of a` | `[re, im]` <-> `[modulus, argument]` |
+| `eq_near` | `eq_near of [a, b, tol]` | compares BOTH components — not a modulus test |
+| `poly_eval` | `poly_eval of [coeffs, z]` | Horner on the monic polynomial |
+| `poly_roots` | `poly_roots of coeffs` | all n complex roots (Durand-Kerner) |
+| `poly_roots_tuned` | `poly_roots_tuned of [coeffs, iters, tol]` | same, with the iteration cap and tolerance exposed |
+
+Coefficient convention, shared with `linalg.charpoly`: `[c1 ... cn]` means
+the monic polynomial `z^n + c1 z^(n-1) + ... + cn`, so `[3, 2]` is
+`z^2 + 3z + 2`. The leading 1 is implied and never appears in the list.
 
 ### lib/calculus.eigs
 
@@ -1463,6 +1489,7 @@ Unit conversions, signal processing (DFT/IDFT, convolution, spectrum), control (
 | `dft` | `dft of signal` | Returns list of [real, imaginary] pairs |
 | `idft` | `idft of spectrum` | Inverse DFT |
 | `magnitude_spectrum` | `magnitude_spectrum of dft_result` | \|X[k]\| = sqrt(re^2 + im^2) |
+| `phase_spectrum` | `phase_spectrum of dft_result` | arg X[k] per bin, radians, via atan2 |
 | `power_spectrum` | `power_spectrum of dft_result` | \|X[k]\|^2 |
 | `convolve` | `convolve of [signal, kernel]` | Linear convolution |
 | `moving_average` | `moving_average of [signal, window]` | Moving average filter |

@@ -0,0 +1,202 @@
+# ============================================================
+# Complex Numbers — arithmetic, polar form, and polynomial roots
+# ============================================================
+#
+# How to use:
+#   import complex
+#
+#   complex.mul of [[3, 4], [1, -2]]      -> [11, -2]
+#   complex.div of [[11, -2], [1, -2]]    -> [3, 4]
+#   complex.mag of [3, 4]                 -> 5
+#   complex.to_polar of [3, 4]            -> [5, 0.9272...]
+#   complex.poly_roots of [3, 2]          -> roots of z^2 + 3z + 2
+#
+# A complex number is a two-element list [re, im]. That is the shape
+# lib/engineering.eigs's `dft` already returns, so its output feeds these
+# functions directly.
+#
+# WHY THIS MODULE EXISTS. Mode analysis, signal processing and root
+# finding are all intrinsically complex-valued, and before this every
+# consumer re-rolled the same five functions on [re, im] pairs
+# (EigenScript#1043, found by the phugoid flight-simulator consumer).
+#
+# AND WHY ITS TESTS ARE DIRECT UNIT CHECKS. A differential oracle is
+# structurally BLIND to this arithmetic: phugoid measured a HALVED `div`
+# surviving a full root-finding differential, because Durand-Kerner
+# self-corrects under a scaled delta — a wrong step size still converges
+# to the right root, just more slowly. Every function here is therefore
+# checked against a hand-computed value in tests/test_complex.eigs, and
+# the halved-div fault is planted there explicitly. Do not replace those
+# with "the consumer still works" checks; that is the exact test which
+# was measured not to work.
+
+# ============================================================
+# ARITHMETIC
+# ============================================================
+
+# add: (a + b)
+define add(a, b) as:
+    return [a[0] + b[0], a[1] + b[1]]
+
+# sub: (a - b)
+define sub(a, b) as:
+    return [a[0] - b[0], a[1] - b[1]]
+
+# mul: (a * b) = (ar*br - ai*bi) + (ar*bi + ai*br)i
+define mul(a, b) as:
+    return [a[0] * b[0] - a[1] * b[1], a[0] * b[1] + a[1] * b[0]]
+
+# div: (a / b), by multiplying through by the conjugate of b.
+#
+# Returns null when b is zero. Division by zero RAISES in EigenScript
+# (measured: "Error line N: division by zero"), so without this guard a
+# zero denominator would kill the caller's program from inside library
+# code. Returning null makes it the caller's decision — and `poly_roots`
+# below DEPENDS on that: at a repeated root every other iterate collides
+# with the current one, the Durand-Kerner denominator is exactly zero, and
+# the step must be skipped rather than raised on. An unguarded `div` makes
+# every polynomial with a repeated root a crash.
+define div(a, b) as:
+    local d is b[0] * b[0] + b[1] * b[1]
+    if d == 0:
+        return null
+    return [(a[0] * b[0] + a[1] * b[1]) / d, (a[1] * b[0] - a[0] * b[1]) / d]
+
+# neg: -a
+define neg(a) as:
+    return [0 - a[0], 0 - a[1]]
+
+# conj: the complex conjugate, re - im*i
+define conj(a) as:
+    return [a[0], 0 - a[1]]
+
+# scale: a * s for a REAL scalar s
+define scale(a, s) as:
+    return [a[0] * s, a[1] * s]
+
+# mag: |a|, the modulus
+define mag(a) as:
+    return sqrt of (a[0] * a[0] + a[1] * a[1])
+
+# mag2: |a|^2, without the square root — cheaper, and exact for integer
+# inputs where `mag` would round.
+define mag2(a) as:
+    return a[0] * a[0] + a[1] * a[1]
+
+# arg: the argument (phase angle) in radians, via atan2, so the range is
+# [-pi, pi] and all four quadrants are correct. Note the closed lower end:
+# a negative real with a NEGATIVE ZERO imaginary part gives -pi, not +pi.
+# The two denote the same angle and the sign is decided by the sign of a
+# zero, so do not assert one branch of it.
+define arg(a) as:
+    return atan2 of [a[1], a[0]]
+
+# ============================================================
+# POLAR FORM
+# ============================================================
+
+# to_polar: [re, im] -> [modulus, argument]
+define to_polar(a) as:
+    return [mag of a, arg of a]
+
+# from_polar: [modulus, argument] -> [re, im]
+define from_polar(r, theta) as:
+    return [r * (cos of theta), r * (sin of theta)]
+
+# ============================================================
+# COMPARISON
+# ============================================================
+
+# eq_near: within tol on BOTH components. Complex equality has two axes
+# and comparing only the modulus would call 3+4i and 5 equal.
+define eq_near(a, b, tol) as:
+    local dr is a[0] - b[0]
+    local di is a[1] - b[1]
+    if dr < 0:
+        dr is 0 - dr
+    if di < 0:
+        di is 0 - di
+    if dr > tol:
+        return 0
+    if di > tol:
+        return 0
+    return 1
+
+# ============================================================
+# POLYNOMIALS
+# ============================================================
+#
+# COEFFICIENT CONVENTION, used by every function below and by
+# linalg.charpoly: a list of n coefficients describes the MONIC degree-n
+# polynomial
+#
+#     z^n + c[0] z^(n-1) + c[1] z^(n-2) + ... + c[n-1]
+#
+# so [3, 2] is z^2 + 3z + 2, whose roots are -1 and -2. The leading 1 is
+# implied and never appears in the list.
+
+# poly_eval: evaluate the monic polynomial at complex z, by Horner.
+define poly_eval(coeffs, z) as:
+    local n is len of coeffs
+    local p is [1.0, 0.0]
+    local k is 0
+    loop while k < n:
+        p is add of [mul of [p, z], [coeffs[k], 0.0]]
+        k is k + 1
+    return p
+
+# poly_roots: all n complex roots of the monic polynomial, by
+# Durand-Kerner. Returns a list of [re, im] pairs.
+define poly_roots(coeffs) as:
+    return poly_roots_tuned of [coeffs, 200, 1e-14]
+
+# poly_roots_tuned: poly_roots with the iteration cap and convergence
+# tolerance exposed. `iters = 0` returns the untouched starting guesses,
+# which is what lets a test assert that the iteration is doing the work
+# rather than the initialisation.
+define poly_roots_tuned(coeffs, iters, tol) as:
+    local n is len of coeffs
+    if n == 0:
+        return []
+    # The standard (0.4 + 0.9i)^k spiral: off the real axis, so a
+    # polynomial with real coefficients cannot trap every iterate in the
+    # reals and miss a conjugate pair.
+    local base is [0.4, 0.9]
+    local roots is []
+    local cur is [1.0, 0.0]
+    local k is 0
+    loop while k < n:
+        append of [roots, cur]
+        cur is mul of [cur, base]
+        k is k + 1
+    local it is 0
+    local done is 0
+    loop while it < iters and done == 0:
+        local maxd is 0.0
+        local i is 0
+        loop while i < n:
+            local zi is roots[i]
+            local pnum is poly_eval of [coeffs, zi]
+            local den is [1.0, 0.0]
+            local j is 0
+            loop while j < n:
+                if j != i:
+                    den is mul of [den, sub of [zi, roots[j]]]
+                j is j + 1
+            # A REPEATED ROOT drives the denominator to zero: every other
+            # iterate has collided with this one. `div` returns null there
+            # instead of raising, and the step is simply skipped — the
+            # iterate is already AT the multiple root, so holding it is
+            # the correct move, and without the skip the whole call would
+            # die with "division by zero" on any repeated root.
+            local delta is div of [pnum, den]
+            if delta != null:
+                roots[i] is sub of [zi, delta]
+                local dm is mag of delta
+                if dm > maxd:
+                    maxd is dm
+            i is i + 1
+        if maxd < tol:
+            done is 1
+        it is it + 1
+    return roots
@@ -12,6 +12,8 @@
 #
 # All SI units unless otherwise noted.
 
+import complex
+
 load_file of "lib/math.eigs"
 
 # ============================================================
@@ -117,21 +119,34 @@ define idft(spectrum) as:
     return result
 
 # |X[k]| = sqrt(re^2 + im^2)
+#
+# `dft` returns [re, im] pairs, which is exactly lib/complex.eigs's
+# representation — so these three are the complex modulus, its square, and
+# the argument, applied down the spectrum. They were hand-rolled here
+# before that module existed, which is the duplication EigenScript#1043
+# was raised about: the stdlib was re-rolling the arithmetic too.
 define magnitude_spectrum(dft_result) as:
     result is []
     for i in range of (len of dft_result):
-        re is dft_result[i][0]
-        im is dft_result[i][1]
-        append of [result, sqrt of (re * re + im * im)]
+        append of [result, complex.mag of (dft_result[i])]
     return result
 
-# |X[k]|^2
+# |X[k]|^2 — no sqrt, so it is exact for integer-valued spectra
 define power_spectrum(dft_result) as:
     result is []
     for i in range of (len of dft_result):
-        re is dft_result[i][0]
-        im is dft_result[i][1]
-        append of [result, re * re + im * im]
+        append of [result, complex.mag2 of (dft_result[i])]
+    return result
+
+# arg X[k], the phase of each bin in radians, (-pi, pi]
+#
+# The missing sibling of the two above: a spectrum has a phase as well as
+# a magnitude, and this was absent because the complex arithmetic to
+# express it did not exist. Uses atan2, so all four quadrants are right.
+define phase_spectrum(dft_result) as:
+    result is []
+    for i in range of (len of dft_result):
+        append of [result, complex.arg of (dft_result[i])]
     return result
 
 # Linear convolution