diff --git a/CHANGELOG.md b/CHANGELOG.md index b3428cf6..f5811c56 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -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 diff --git a/docs/ARCHITECTURE.md b/docs/ARCHITECTURE.md index e90b2ca9..a15c608b 100644 --- a/docs/ARCHITECTURE.md +++ b/docs/ARCHITECTURE.md @@ -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 diff --git a/docs/STDLIB.md b/docs/STDLIB.md index 8123792a..db0a9f26 100644 --- a/docs/STDLIB.md +++ b/docs/STDLIB.md @@ -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 | diff --git a/lib/complex.eigs b/lib/complex.eigs new file mode 100644 index 00000000..aca6c69b --- /dev/null +++ b/lib/complex.eigs @@ -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 diff --git a/lib/engineering.eigs b/lib/engineering.eigs index 835a4278..b2b5f705 100644 --- a/lib/engineering.eigs +++ b/lib/engineering.eigs @@ -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 diff --git a/lib/linalg.eigs b/lib/linalg.eigs index eed1ec3c..284b4238 100644 --- a/lib/linalg.eigs +++ b/lib/linalg.eigs @@ -12,12 +12,19 @@ # linalg.vec_cross of [[1,0,0], [0,1,0]] # # Matrices are lists of lists (row-major). Vectors are flat lists. +# +# `charpoly` and `eigenvalues` return complex results, so this module +# imports lib/complex.eigs. It is the first stdlib module to depend on +# another; every lib/*.eigs is installed together by the Makefile, so the +# dependency travels with it. # ============================================================ # Matrix Construction # ============================================================ # mat_new: create rows x cols matrix filled with fill value +import complex + define mat_new(rows, cols, fill) as: result is [] for i in range of rows: @@ -416,6 +423,12 @@ define least_squares(A, b) as: # eigenvalues_2x2: eigenvalues of 2x2 matrix via characteristic polynomial # det(A - lambda*I) = 0 => lambda^2 - trace*lambda + det = 0 +# +# REAL EIGENVALUES ONLY. A negative discriminant means a complex conjugate +# pair, and this returns `null` for it — so a rotation matrix, whose +# eigenvalues are the perfectly ordinary +-i, reads here as "no answer". +# Use `eigenvalues` below for anything that can oscillate; it returns +# complex pairs and handles the 2x2 case too. define eigenvalues_2x2(A) as: tr is A[0][0] + A[1][1] det is A[0][0] * A[1][1] - A[0][1] * A[1][0] @@ -454,3 +467,57 @@ define eigenvectors_2x2(A) as: else: append of [vecs, [1, 0]] return vecs + + +# ============================================================ +# CHARACTERISTIC POLYNOMIAL AND GENERAL EIGENVALUES +# ============================================================ +# +# Added for EigenScript#1042, found by the phugoid flight-simulator +# consumer: extracting flight modes needs the eigenvalues of a 4x4 state +# matrix, and this module stopped at `eigenvalues_2x2` (real only) while +# lib/numerics.eigs offered only `power_iteration` (dominant, real only). + +# charpoly: coefficients of the characteristic polynomial of A, by the +# Faddeev-LeVerrier recurrence +# +# M_0 = I, c_k = -tr(A M_(k-1)) / k, M_k = A M_(k-1) + c_k I +# +# Returns [c1, c2, ... cn] for the MONIC polynomial +# +# lambda^n + c1 lambda^(n-1) + ... + cn +# +# i.e. the same convention lib/complex.eigs's poly_eval/poly_roots use, so +# `complex.poly_roots of (linalg.charpoly of A)` composes directly. The +# leading 1 is implied and is not in the list. +# +# Works for any square n. Returns null for a non-square or empty matrix +# rather than indexing off the end of a row. +define charpoly(A) as: + local n is len of A + if n == 0: + return null + if (len of A[0]) != n: + return null + local M is mat_identity of n + local coeffs is [] + local k is 1 + loop while k <= n: + local AM is mat_mul of [A, M] + local ck is (0 - (mat_trace of AM)) / k + append of [coeffs, ck] + M is mat_add of [AM, mat_scale of [mat_identity of n, ck]] + k is k + 1 + return coeffs + +# eigenvalues: all n eigenvalues of a square matrix A, as complex +# [re, im] pairs — the characteristic polynomial's roots. +# +# Unlike `eigenvalues_2x2` this does not refuse a complex spectrum, which +# is the whole point: oscillatory systems have conjugate pairs, and their +# imaginary parts are the frequencies. +define eigenvalues(A) as: + local coeffs is charpoly of A + if coeffs == null: + return null + return complex.poly_roots of coeffs diff --git a/tests/run_all_tests.sh b/tests/run_all_tests.sh index 5f5f9daa..776b258c 100755 --- a/tests/run_all_tests.sh +++ b/tests/run_all_tests.sh @@ -4785,6 +4785,7 @@ check_eigs_suite "numerics" test_numerics.eigs "All tests passed." 1 check_eigs_suite "optimize" test_optimize.eigs "All tests passed." 1 check_eigs_suite "simulation" test_simulation.eigs "All tests passed." 1 check_eigs_suite "linalg" test_linalg.eigs "All tests passed." 1 +check_eigs_suite "complex" test_complex.eigs "All tests passed." 1 check_eigs_suite "probability" test_probability.eigs "All tests passed." 1 check_eigs_suite "biology" test_biology.eigs "All tests passed." 1 check_eigs_suite "calculus" test_calculus.eigs "All tests passed." 1 diff --git a/tests/test_complex.eigs b/tests/test_complex.eigs new file mode 100644 index 00000000..d049b0dc --- /dev/null +++ b/tests/test_complex.eigs @@ -0,0 +1,168 @@ +# ============================================================ +# Complex Number Tests +# ============================================================ +# +# EVERY VALUE HERE IS HAND-COMPUTED. That is deliberate and it is the +# whole design of this file. +# +# A differential oracle is structurally BLIND to this arithmetic. The +# phugoid consumer 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. +# So "the root finder still gets the right answer" proves nothing about +# `div`, and the halved-div fault is planted explicitly below. +# +# If you add a function to lib/complex.eigs, add its hand-computed check +# here. Do not check it by running a consumer. +load_file of "lib/test.eigs" +import complex + +print of "=== COMPLEX NUMBER TESTS ===" + +# ---- Arithmetic, against values computed by hand ---- +print of "--- Arithmetic ---" + +# (3+4i) + (1-2i) = 4+2i +s is complex.add of [[3, 4], [1, -2]] +assert_eq of [s[0], 4, "add re"] +assert_eq of [s[1], 2, "add im"] + +# (3+4i) - (1-2i) = 2+6i +d is complex.sub of [[3, 4], [1, -2]] +assert_eq of [d[0], 2, "sub re"] +assert_eq of [d[1], 6, "sub im"] + +# (3+4i)(1-2i) = 3 - 6i + 4i - 8i^2 = 3 - 2i + 8 = 11 - 2i +m is complex.mul of [[3, 4], [1, -2]] +assert_eq of [m[0], 11, "mul re"] +assert_eq of [m[1], -2, "mul im"] + +# i * i = -1 exactly: the sign of the i^2 term, which a dropped minus +# would silently invert. +ii is complex.mul of [[0, 1], [0, 1]] +assert_eq of [ii[0], -1, "i*i re"] +assert_eq of [ii[1], 0, "i*i im"] + +# (11-2i) / (1-2i) = 3+4i — the inverse of the mul above +q is complex.div of [[11, -2], [1, -2]] +assert_near of [q[0], 3, 1e-12, "div re"] +assert_near of [q[1], 4, 1e-12, "div im"] + +# THE PLANTED-FAULT DISCRIMINATOR. (1+i)/(2+0i) = 0.5+0.5i. A `div` whose +# result is halved returns 0.25+0.25i and still points the same DIRECTION, +# which is exactly why the root-finding differential could not see it. +# This check can. +h is complex.div of [[1, 1], [2, 0]] +assert_near of [h[0], 0.5, 1e-12, "div magnitude re (halved-div discriminator)"] +assert_near of [h[1], 0.5, 1e-12, "div magnitude im (halved-div discriminator)"] + +# Division by zero returns null rather than saturating at 1e308. +assert_true of [(complex.div of [[1, 1], [0, 0]]) == null, "div by zero is null"] + +# -a and conj(a) differ: negation flips BOTH components, conjugation only +# the imaginary one. Confusing them survives any magnitude-only check. +n is complex.neg of [3, 4] +assert_eq of [n[0], -3, "neg re"] +assert_eq of [n[1], -4, "neg im"] +c is complex.conj of [3, 4] +assert_eq of [c[0], 3, "conj re (unchanged)"] +assert_eq of [c[1], -4, "conj im (flipped)"] + +sc is complex.scale of [[3, 4], 2] +assert_eq of [sc[0], 6, "scale re"] +assert_eq of [sc[1], 8, "scale im"] + +# ---- Modulus and argument ---- +print of "--- Modulus and argument ---" + +# |3+4i| = 5, the 3-4-5 triangle +assert_eq of [complex.mag of [3, 4], 5, "mag 3-4-5"] +assert_eq of [complex.mag2 of [3, 4], 25, "mag2 is the square, no sqrt"] +assert_eq of [complex.mag of [0, 0], 0, "mag of zero"] + +# arg(i) = pi/2, arg(-1) = pi, arg(1) = 0 — the quadrants atan alone +# gets wrong and atan2 gets right. +assert_near of [complex.arg of [0, 1], 1.5707963267948966, 1e-12, "arg i = pi/2"] +assert_near of [complex.arg of [-1, 0], 3.141592653589793, 1e-12, "arg -1 = pi"] +assert_eq of [complex.arg of [1, 0], 0, "arg 1 = 0"] + +# ---- Polar form ---- +print of "--- Polar form ---" +p is complex.to_polar of [3, 4] +assert_eq of [p[0], 5, "to_polar modulus"] +assert_near of [p[1], 0.9272952180016122, 1e-12, "to_polar argument"] + +# Round trip back to rectangular. +r is complex.from_polar of [5, 0.9272952180016122] +assert_near of [r[0], 3, 1e-12, "from_polar re"] +assert_near of [r[1], 4, 1e-12, "from_polar im"] + +# ---- Comparison ---- +print of "--- Comparison ---" +assert_eq of [complex.eq_near of [[3, 4], [3, 4], 1e-9], 1, "eq_near identical"] +assert_eq of [complex.eq_near of [[3, 4], [3.0000000001, 4], 1e-6], 1, "eq_near within tol"] +assert_eq of [complex.eq_near of [[3, 4], [3, 4.1], 1e-6], 0, "eq_near differs in im"] +# Two axes, not one: 3+4i and 5+0i have the SAME modulus and are not equal. +assert_eq of [complex.eq_near of [[3, 4], [5, 0], 1e-6], 0, "eq_near is not a modulus comparison"] + +# ---- Polynomials ---- +print of "--- Polynomials ---" + +# z^2 + 3z + 2 = (z+1)(z+2); evaluate at its roots and at a known point. +assert_eq of [(complex.poly_eval of [[3, 2], [-1, 0]])[0], 0, "poly_eval at root -1"] +assert_eq of [(complex.poly_eval of [[3, 2], [-2, 0]])[0], 0, "poly_eval at root -2"] +# p(1) = 1 + 3 + 2 = 6 +assert_eq of [(complex.poly_eval of [[3, 2], [1, 0]])[0], 6, "poly_eval at 1"] + +# Roots of z^2 + 3z + 2: -1 and -2 (order is not guaranteed, so check the +# SET via the symmetric functions below as well). +rr is complex.poly_roots of [3, 2] +assert_eq of [len of rr, 2, "poly_roots returns n roots"] + +# VIETA, which checks the whole root set at once rather than one root: +# for a monic polynomial, sum of roots = -c1 and product = (-1)^n * cn. +sum is [0, 0] +prod is [1, 0] +for k in range of (len of rr): + sum is complex.add of [sum, rr[k]] + prod is complex.mul of [prod, rr[k]] +assert_near of [sum[0], -3, 1e-9, "Vieta: sum of roots = -c1"] +assert_near of [prod[0], 2, 1e-9, "Vieta: product of roots = (-1)^n cn"] + +# RESIDUAL: every returned root must actually be a root. +for k in range of (len of rr): + res is complex.poly_eval of [[3, 2], rr[k]] + assert_near of [complex.mag of res, 0, 1e-9, "residual at returned root"] + +# A COMPLEX conjugate pair, which is the case a real-only solver cannot +# express at all: z^2 + 1 has roots +i and -i. +ci is complex.poly_roots of [0, 1] +assert_near of [complex.mag of ci[0], 1, 1e-9, "z^2+1 root has modulus 1"] +assert_near of [ci[0][0], 0, 1e-9, "z^2+1 root is purely imaginary"] +assert_near of [(complex.add of [ci[0], ci[1]])[1], 0, 1e-9, "the two roots are conjugates"] + +# A REPEATED root: (z-1)^2 = z^2 - 2z + 1. The Durand-Kerner denominator +# goes to zero here, and an unguarded division would fling an iterate to +# 1e308 instead of holding it at the multiple root. +rep is complex.poly_roots of [-2, 1] +assert_near of [rep[0][0], 1, 1e-6, "repeated root 0"] +assert_near of [rep[1][0], 1, 1e-6, "repeated root 1"] + +# A CUBIC, so the degree is not hard-coded at 2: z^3 - 6z^2 + 11z - 6 has +# roots 1, 2, 3. +cub is complex.poly_roots of [-6, 11, -6] +assert_eq of [len of cub, 3, "cubic returns 3 roots"] +csum is [0, 0] +for k in range of 3: + csum is complex.add of [csum, cub[k]] +assert_near of [csum[0], 6, 1e-8, "cubic Vieta: roots sum to 6"] + +# THE ITERATION MUST BE DOING THE WORK. With iters = 0 the starting +# guesses are returned untouched, and they are NOT the answer — without +# this, a solver that returned its initial spiral would pass any check +# loose enough to accept it. +g0 is complex.poly_roots_tuned of [[3, 2], 0, 1e-14] +assert_eq of [g0[0][0], 1.0, "iters=0 returns the untouched first guess"] +assert_true of [(complex.mag of (complex.poly_eval of [[3, 2], g0[1]])) > 0.5, "the starting guesses are not roots"] + +test_summary of null diff --git a/tests/test_engineering.eigs b/tests/test_engineering.eigs index bf0e09f0..d5cadcb2 100644 --- a/tests/test_engineering.eigs +++ b/tests/test_engineering.eigs @@ -137,6 +137,28 @@ assert_near of [mag[2], 1, 0.001, "magnitude [1,0,0,0] k=2"] pwr is power_spectrum of spec assert_near of [pwr[0], 1, 0.001, "power [1,0,0,0] k=0"] +# phase_spectrum (#1043): the missing sibling of the two above, available +# once lib/complex.eigs made the argument expressible. A real even signal +# has zero phase in every bin; a shifted impulse does not, which is what +# separates a real phase calculation from a stub returning zeros. +ph is phase_spectrum of spec +assert_near of [ph[0], 0, 1e-9, "phase of [1,0,0,0] k=0 is 0"] +assert_near of [ph[1], 0, 1e-9, "phase of [1,0,0,0] k=1 is 0"] + +# A one-sample shift puts -pi/2 in bin 1: X[k] = e^(-2 pi i k / 4), so +# k=1 has argument -pi/2. atan2 is what gets the quadrant right here. +shifted is dft of [0, 1, 0, 0] +phs is phase_spectrum of shifted +assert_near of [phs[1], -1.5707963267948966, 1e-9, "phase of a shifted impulse k=1 is -pi/2"] +# The Nyquist bin is real and negative, so its phase is pi — but its +# imaginary part is a signed ZERO, and atan2 returns -pi for -0. The sign +# there is not a meaningful distinction, so grade the magnitude of the +# angle rather than pretending one branch is correct. +absph is phs[2] +if absph < 0: + absph is 0 - absph +assert_near of [absph, 3.141592653589793, 1e-9, "phase of a shifted impulse k=2 is pi up to the sign of zero"] + # IDFT roundtrip: IDFT(DFT([1,0,0,0])) = [1,0,0,0] recovered is idft of spec assert_near of [recovered[0], 1, 0.001, "IDFT roundtrip x[0]"] diff --git a/tests/test_linalg.eigs b/tests/test_linalg.eigs index 7e29f6f5..4de34732 100644 --- a/tests/test_linalg.eigs +++ b/tests/test_linalg.eigs @@ -3,6 +3,7 @@ # ============================================================ load_file of "lib/test.eigs" import linalg +import complex print of "=== LINEAR ALGEBRA TESTS ===" @@ -209,5 +210,67 @@ c is linalg.mat_col of [M, 2] assert_eq of [c[0], 3, "col 2 [0]"] assert_eq of [c[1], 6, "col 2 [1]"] +# ---- Characteristic polynomial and general eigenvalues (#1042) ---- +print of "--- charpoly / eigenvalues ---" + +# diag(1,2,3): (L-1)(L-2)(L-3) = L^3 - 6L^2 + 11L - 6 +cp is linalg.charpoly of [[1, 0, 0], [0, 2, 0], [0, 0, 3]] +assert_eq of [len of cp, 3, "charpoly returns n coefficients"] +assert_near of [cp[0], -6, 1e-12, "charpoly diag c1"] +assert_near of [cp[1], 11, 1e-12, "charpoly diag c2"] +assert_near of [cp[2], -6, 1e-12, "charpoly diag c3"] + +# The COMPANION matrix of L^2 + 3L + 2 must give back its own coefficients. +cc is linalg.charpoly of [[0, 1], [-2, -3]] +assert_near of [cc[0], 3, 1e-12, "charpoly companion c1"] +assert_near of [cc[1], 2, 1e-12, "charpoly companion c2"] + +# Non-square and empty are refused, not indexed off the end of a row. +assert_true of [(linalg.charpoly of [[1, 2, 3], [4, 5, 6]]) == null, "charpoly refuses non-square"] +assert_true of [(linalg.charpoly of []) == null, "charpoly refuses empty"] + +# Real spectrum, recovered as complex pairs with zero imaginary part. +ed is linalg.eigenvalues of [[1, 0, 0], [0, 2, 0], [0, 0, 3]] +assert_eq of [len of ed, 3, "eigenvalues returns n values"] +esum is 0 +for k in range of 3: + esum is esum + ed[k][0] +assert_near of [esum, 6, 1e-8, "eigenvalues of diag sum to the trace"] + +# THE CASE eigenvalues_2x2 CANNOT EXPRESS. A rotation matrix has +# eigenvalues +-i; eigenvalues_2x2 returns null for it because the +# discriminant is negative, and this is why the general path exists. +assert_true of [(linalg.eigenvalues_2x2 of [[0, -1], [1, 0]]) == null, "eigenvalues_2x2 refuses a complex spectrum"] +er is linalg.eigenvalues of [[0, -1], [1, 0]] +assert_near of [er[0][0], 0, 1e-9, "rotation eigenvalue is purely imaginary"] +assert_near of [(er[0][1]) * (er[1][1]), -1, 1e-9, "rotation eigenvalues are +i and -i"] + +# A 4x4 with a KNOWN complex spectrum: block-diagonal with a rotation +# block (+-i) and a damped block (-1 +- 2i). This is the shape the issue +# was raised for — a 4x4 state matrix whose modes are conjugate pairs. +B4 is [[0, -1, 0, 0], [1, 0, 0, 0], [0, 0, -1, -2], [0, 0, 2, -1]] +e4 is linalg.eigenvalues of B4 +assert_eq of [len of e4, 4, "4x4 returns 4 eigenvalues"] + +# Graded by the two identities that tie the spectrum back to the MATRIX +# rather than to the solver: sum of eigenvalues = trace, product = det. +# mat_trace and mat_det are tested independently above, so this is a real +# cross-check and not a restatement of charpoly. +s4 is [0, 0] +p4 is [1, 0] +for k in range of 4: + s4 is complex.add of [s4, e4[k]] + p4 is complex.mul of [p4, e4[k]] +assert_near of [s4[0], linalg.mat_trace of B4, 1e-8, "sum of eigenvalues = trace"] +assert_near of [s4[1], 0, 1e-8, "the spectrum is closed under conjugation (imag sum 0)"] +assert_near of [p4[0], linalg.mat_det of B4, 1e-7, "product of eigenvalues = determinant"] + +# ...and every eigenvalue must actually annihilate the characteristic +# polynomial, which catches a solver that returns plausible-looking +# numbers with the right sum. +cp4 is linalg.charpoly of B4 +for k in range of 4: + assert_near of [complex.mag of (complex.poly_eval of [cp4, e4[k]]), 0, 1e-8, "eigenvalue is a root of charpoly"] + print of "" test_summary of null