Refactor numerical tower and fix critical bugs - #1
Open
longmathemagician wants to merge 14 commits into
Open
Conversation
…n proposals Documents verified defects (matrix dispatch bitwise-AND bug, Complex magnitude/recip/parse bugs, failing cubic roots path), proposes a Real/Scalar rework of the numerical tower, a migration off generic_const_exprs onto per-size inherent impls, an audit of stale feature gates on current nightly, and a min_specialization-based matrix kernel architecture verified against 1.98.0-nightly. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_012bGPMXfU7c6jnHVSX33DJW
Implements REVIEW.md Phase 1 (sections 1.1-1.6 and 5): - matrix: kernel dispatch and determinant guards used bitwise & instead of &&, so e.g. Matrix<f64,2,3> * 2x2 identity ran the Strassen kernel and zeroed the third row; determinant returned numbers for non-square matrices. Now M == 2 && N == 2 && O == 2 etc.; non-square determinant hits the todo!() arm. Added regression tests for non-square products and a should_panic determinant test. - complex: magnitude computed (re^2+im^2)^2 == |z|^4; final powi(2) is now sqrt(). Added magnitude/normalize tests. - algebra/complex: impl_field! generated recip(self) = Self::recip(self), an infinite recursion for types without an inherent recip (stack overflow for Complex). The macro now calls <T>::recip unambiguously, and c32/c64 get a hand-written Field impl with recip = conj(z)/|z|^2. - complex: replaced From<&str> (which silently parsed "3 - 4i" as 0+0i) with core::str::FromStr and a public ParseComplexError. Handles "a + bi", "a - bi", "a", "bi", "-a - bi" with tolerant whitespace and round-trips Display output. Updated demo example and README. - polynomial: wired the cubic arm of roots() to the generic Blinn solver (previously fell through to NaN; roots_3_generic failed). Documented and enforced the ordering contract: finite roots ascending, non-finite entries last; solver modules keep their native order. - natural/integer: powi(power as u32) turned negative exponents into huge ones; now asserts power >= 0 with a clear message (should_panic tests added). - lib: removed stale feature gates (associated_type_bounds, const_float_bits_conv, generic_arg_infer stable; more_qualified_paths and min_specialization unused) and the removed soft_unstable lint. generic_const_exprs + incomplete_features stay until Phase 3. - solvers: scoped allow(non_snake_case) for paper-notation variables, dropped unneeded mut/unused params, removed superseded commented-out code. cargo build and cargo test are now warning-free. - toolchain: pinned rust-toolchain.toml to nightly-2026-07-02. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_012bGPMXfU7c6jnHVSX33DJW
Numerical tower (REVIEW.md §2): - New Real trait (src/real.rs, replaces Float): Field + PartialOrd with float consts (EPSILON/INFINITY/NAN/MIN/MAX, DIGITS/MANTISSA_DIGITS/ RADIX as u32, MIN_EXP/MAX_EXP as i32), by-value math methods, and from_u32 for building small solver constants. Implemented for f32/f64. - New Scalar trait (src/scalar.rs): Field with a real-valued norm (norm_sqr/norm -> Self::Real), conj, mul_add (fused for reals), is_nan/is_finite, and Mul/Div<Self::Real> supertraits. Blanket impl for Real types plus impl for Complex<F>. - Natural is now unsigned machine integers only (Semiring + Ord + Eq, BITS: u32, pow(u32)); Integer is signed machine integers only (Ring + Ord + Eq, plus abs). Floats and Complex no longer implement them; Complex loses PartialOrd and all todo!() Float methods. - Deleted the Number identity trait; every T: Number<Type = T> + Float bound is now T: Scalar (or T: Real where genuinely real-only). - Element drops the Display supertrait; Display bounds moved to the formatting impls. Group loses the redundant Neg method. Monoid gains iter::Sum. Each impl_* macro implements exactly one trait; the stack_* family macros compose them. Complex algebra impls are now generic over F: Real. Const generics (REVIEW.md §3): - Removed #![feature(generic_const_exprs)] and incomplete_features; the crate now builds and tests clean on stable Rust. - New counted Roots<T, MAX> (src/roots.rs) with as_slice/len/Deref/ IntoIterator/PartialEq on the live prefix; no NaN sentinels in the public roots API. - Polynomial::roots is now per-size inherent impls for N = 2..=5 (linear/quadratic/cubic/quartic) with a documented contract: ascending order, len == number of real roots found. The N + 0^(N-1) - 1 arrays, [(); N]: bounds, and NaN-padded returns are gone. - Polynomial::eval is a single Horner fold over T: Scalar (complex coefficients work); per-degree eval helpers deleted. new is the only constructor plus a real From<[T; N]> impl. - Yuksel solvers genericized from f64-only to T: Real; the quartic path is reachable for the first time and tested. Blinn is now plain generic free functions (PhantomData struct deleted). - Matrix<T, M, N> is now conventionally M rows x N cols with row-major storage [[T; N]; M]; Mul has the standard (M,K) x (K,N) -> (M,N) signature. IDENTITY/determinant/inverse live on Matrix<T, N, N>, so non-square usage is a compile error (compile_fail doctest included). determinant uses closed forms for N <= 3 and LU with norm_sqr partial pivoting above; inverse is Gauss-Jordan returning Option; transpose is implemented for all shapes. - Vector is generic over Scalar: magnitude returns T::Real (a complex vector's magnitude is an f64), normalize divides by the real magnitude via Div<T::Real>. Tests: 44 -> 70 unit tests plus a compile_fail doctest, all green on both the pinned nightly and stable; demo example updated and passing. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_012bGPMXfU7c6jnHVSX33DJW
…tmul kernels
Task 1 — mul_add_fast (measured 5.7x perf trap):
- Scalar gains a provided `mul_add_fast`: hardware FMA on targets that have
it (cfg target_feature=fma or aarch64), plain `a*x + b` elsewhere, so hot
loops never hit the indirect libm software-fma call (8.19 ns vs 1.45 ns
per Horner step on non-FMA x86-64).
- `mul_add` stays always-fused as the *precision* tool (Blinn's discriminant
depends on it) and its docs now state the non-FMA libm cost.
- Polynomial::eval's Horner fold and the naive matrix-multiply accumulation
now use mul_add_fast; a new test pins mul_add_fast == mul_add == plain on
exactly-representable values for f32/f64/c64.
Task 2 — cargo features and no_std (REVIEW.md #6):
- Features: default = ["std"]; `libm` (optional dep, fetched as v0.2.16)
supplies the float math backend for no_std builds; std wins if both are
enabled; compile_error! if neither backend is selected.
- lib.rs is `#![no_std]` without the std feature; the Real impl macro picks
std intrinsics or libm functions per cfg (powi's no_std fallback is
exponentiation by squaring, matching powi's multiply-sequence semantics
instead of a transcendental libm::pow).
- Display for Vector/Matrix/Polynomial writes straight to the Formatter
(no String/format!/pop, no alloc); Complex FromStr reparsed core-only
(trim + slice scanning, no String); ParseComplexError now implements
core::error::Error; deleted the dead display.rs (format_f64).
- Verified: `--no-default-features --features libm` builds warning-free on
nightly and stable, and cross-builds for thumbv7em-none-eabihf.
Task 3 — matmul kernel architecture (REVIEW.md #4.1):
- Default builds: Mul is only the naive triple loop with mul_add_fast
accumulation; the crate stays stable-Rust-compatible.
- `--features specialization` (nightly): min_specialization is enabled and
Mul delegates to an internal Gemm<Rhs> trait in src/matrix/kernels.rs with
a `default fn` naive impl plus specialized impls for the concrete square
sizes: 2x2 Strassen, 3x3 Laderman, 4x4 AlphaTensor-style. Dimension
dispatch is type-checked (a 3x2 * 2x2 product can only match the default
impl); the module is out-of-line so stable builds never parse the
unstable `default fn` syntax. Kernel bodies moved unchanged; module docs
record the stability/perf caveats.
- New test compares Mul against a plain reference loop for 2x2/3x3/4x4 and
non-square/degenerate shapes on small-integer matrices (exact equality),
in both feature configurations.
Acceptance (all warning-free): cargo test on nightly-2026-07-02 (default and
--features specialization) and stable 1.94.1 — 72+1 tests green each;
lib builds for {libm}, {libm,specialization}, and thumbv7em-none-eabihf+libm;
examples/demo passes.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_012bGPMXfU7c6jnHVSX33DJW
Design spec for building toward a polyhedral homotopy continuation solver: ascending coefficient order, Algebra<T> evaluation trait, Dual/DualN forward AD, Real/Complex completeness with Smith division, LU factorization, Ring-relaxed containers, lattice (SNF/HNF) module, running-error-bound evaluation, and multivariate sparse polynomials. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_012bGPMXfU7c6jnHVSX33DJW
…ce, mvpoly Implements HOMOTOPY.md §5.1–§5.8, the tower refactor toward a polyhedral homotopy continuation solver: - §5.1 Polynomial flipped to ascending coefficient order (c[i]·x^i): Horner folds from the top via iter().rev(), Display still prints highest-degree first, Blinn/Yuksel indexing and all tests/demo updated. - §5.2 Algebra<T> marker trait (Ring + Mul<T> + Add<T> + From<T>) with the blanket self-instance and Complex<F>: Algebra<F>; Polynomial::eval_at evaluates one Horner body over any Algebra (scalar, complex, dual points); monomorphic eval keeps its mul_add_fast fast path. - §5.3 src/dual.rs: Dual<T> and DualN<T, K> forward AD as ring elements — full tower impls (Field recip = (1/v, −d/v²)), mixed-scalar ops, Algebra instances, nested Dual<Complex<F>>/DualN<Complex<F>, K>: Algebra<F>, and chain-rule lifts (sqrt/sin/cos/sin_cos/tan/exp/ln/powi/mul_add) on Real. - §5.4 Real gains PI/TAU/E and exp/ln/powf/signum/min/max/clamp in both std and libm backends; Complex gains exp/ln/arg/powf/from_polar, by-value powi via exponentiation-by-squaring, nth_root_of_unity, and Smith's algorithm replaces textbook complex division (Div/DivAssign, f32/f64 LHS, recip) so 1e±300-magnitude quotients survive. - §5.5 Lu<T, N> packed factorization (partial pivoting by norm_sqr) as the linear-solve primitive: Matrix::lu/solve, Lu::solve/determinant; the N>3 determinant arm and inverse now share this one pivoting code path. - §5.6 Matrix/Vector structural ops (construction, ZERO/IDENTITY, add/sub/ neg, scalar mul, transpose, matmul) relaxed to T: Ring — integer exponent matrices are first-class; norms/determinant/inverse/lu/solve stay Scalar. New src/lattice.rs: Smith and Hermite normal forms over Matrix<i64, M, N> by gcd reduction, returning unimodular transforms. - §5.7 Polynomial::eval_with_error: Horner with Higham's running error bound, turning stopping rules into "below evaluation noise". - §5.8 src/mvpoly.rs: Monomial/MPoly/MSystem sparse multivariate types with eval_at over any Algebra, structural partial, eval_grad via DualN, MSystem::eval/eval_jacobian, and allocation-free Display. Deviations from the notes, both forced by min_specialization (it rejects specializing impls that add non-marker bounds like Scalar over a Ring default): the Gemm kernels' bounds are relaxed Scalar→Ring (algorithms untouched — they only use +,−,*,ZERO) so Ring matmul exists under the specialization feature too, and matmul accumulation is plain mul+add in all configurations (the §5.6-accepted trade; LLVM still contracts to FMA). Tests grow 72 → 132 (+1 doctest). Full feature matrix green with zero warnings: nightly + stable cargo test, --features specialization, --no-default-features --features libm (also with specialization, and for thumbv7em-none-eabihf), and the demo example gains an eval_at + Dual derivative + DualN Jacobian + Newton-corrector showcase. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_012bGPMXfU7c6jnHVSX33DJW
…s, CI - Containers: hermitian Vector::dot, 3-D cross, Index/IndexMut, Default, From<array>, scalar Div/Neg; Matrix (row,col) indexing, Default, From<array>, Matrix x Vector; scalar-left multiplication macro-stamped for f32/f64/c32/c64 on both Vector and Matrix. - tests/laws.rs: semiring/group/field law macros instantiated for u32, i32, f32, f64, c32, c64, Dual<f64> on exactly-representable samples. - tests/oracle.rs: proptest oracles vs num-complex (ring ops, Smith division at 2^+/-950 magnitudes via exact power-of-two rescaling, exp, ln) and nalgebra (matmul/det/inverse, 4x4 and 5x5); root-finder residuals, LU solve residuals, Dual vs central differences, SNF invariants. Dev-deps pinned (=1.11.0 / =0.4.6 / =0.35.0). - tools/check_codegen.sh + tools/codegen-probe: fails if any call lands in Polynomial::<f64,4>::eval / eval_at bodies, default and +fma configs (follows MergeFunctions symbol aliases); workspace-excluded. - Docs: #![warn(missing_docs)] + #![forbid(unsafe_code)], every public item documented, crate-level homotopy example; cargo doc warning-free. - CI: stable fmt/clippy/test, pinned-nightly test + specialization lane, no_std libm builds (host + thumbv7em-none-eabihf), advisory codegen guard. - README rewritten around the tower/AD/LU/lattice/mvpoly APIs with the ascending-coefficient convention called out; snippets sourced from examples/demo.rs. examples/bench.rs micro-benchmarks (std::time only). - Clippy-clean on stable 1.94 and pinned nightly with -D warnings (LU/lattice row ops rewritten iterator-style, test lint fixes). - Version 0.2.0, rust-version 1.94. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_012bGPMXfU7c6jnHVSX33DJW
The stable and no_std jobs invoked bare cargo, which resolves the toolchain through rust-toolchain.toml's nightly pin rather than the action-installed default — so cargo fmt ran on a nightly without the rustfmt component, and the thumbv7em target was added to stable while the build used the pinned nightly. Set RUSTUP_TOOLCHAIN at the job level (the env var outranks the toolchain file), correct the header comment that claimed the action's default takes precedence, and bump actions/checkout to v5 for the Node 20 deprecation. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_012bGPMXfU7c6jnHVSX33DJW
…al starts
New std-gated module tree src/solvers/homotopy/ implementing the offline
half of a Huber–Sturmfels polyhedral homotopy on top of the Phase 5–6
tower (lattice SNF, Ring matrices, LU solve, mvpoly, Complex polar ops):
- support.rs: Support<NV> (deduplicated exponent-vector sets, extracted
from MPoly/MSystem by skipping zero-coefficient padding terms) and
Lifting<F> (per-point lift values; deterministic Knuth-MMIX LCG mapping
the top 24 state bits to [0,1) — same seed, same subdivision), plus
random_liftings drawing one stream across all supports.
- cells.rs: fine mixed-cell enumeration by naive Π C(|A_i|,2) edge-tuple
search — level system solved per tuple via Matrix::lu/solve, strict
minimality checked against a documented 1e-9 relative tolerance band.
Ambiguous margins on tuples not otherwise rejected abort with the
first-class GenericityError ("re-lift with a new seed") instead of
guessing. MixedCell stores edges, normal α, and the i64 edge matrix V;
volume() computes |det V| exactly as the product of the SNF diagonal.
mixed_volume(supports, seed) does lift→cells→sum.
- start.rs: binomial start systems x^V = β with β_i = −c_{a_i}/c_{b_i},
solved exactly through (U, S, W) = smith_normal_form(V): raise to U
(γ = β^U), substitute x = y^W to reach y_i^{s_i} = γ_i, enumerate the
s_i-th roots per coordinate, map back through x = y^W; |det V|
pairwise-distinct solutions per cell.
- complex.rs: additive Complex::nth_roots(n) iterator (radius through
the overflow-safe ln; no_std-friendly, ExactSizeIterator) with tests.
Complex::powi already handled negative exponents — no change needed.
Anchor results: unit-simplex pair MV 1 (one cell), dense conic pair MV 4,
sparse trinomial pair MV 2 (< Bézout 4), three 3-var linear supports MV 1;
same seed reproduces identical cells, different seeds the same volume;
end-to-end lift→cells→starts residuals < 1e-9. 24 new tests (200 total);
fmt/clippy/test matrix, libm/thumbv7em builds, docs, and codegen guard
all clean on stable and the pinned nightly.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_012bGPMXfU7c6jnHVSX33DJW
The online half of the Huber-Sturmfels pipeline, completing the solver:
- CellHomotopy<F, NV, MAXT>: precomputed per-term coefficients and shifted
levels e_{i,a} = <a,alpha> + w_i(a) - m_i for one mixed cell, with
allocation-free eval / eval_jacobian / dt and a write-into-buffer
system_at variant so the tracker reuses one MSystem. Edge-term levels are
pinned to exactly zero, so H(.,0) is the binomial start system and
H(.,1) is the target bit-for-bit. Levels are globally normalized by
t -> t^(1/e_min) (same path set, different speed): raw lifted levels can
be ~0.07, which parks the binomial regime below t ~ 1e-14 and made the
dense-conic paths untrackable before normalization.
- track_path: Euler predictor (J*ydot = -dH/dt) + Newton corrector at
fixed t with step doubling/halving, relative update-norm convergence,
a corrector-only first step (dH/dt ~ t^(e-1) is singular at t = 0), an
exact landing on t = 1, and a terminal Newton polish against F itself.
Honest statuses: Converged / MinStepReached / MaxStepsReached /
SingularJacobian / Diverged, with t_reached, steps, and Newton counts in
PathResult. Allocation-free by construction for future no_std exposure.
- solve(system, seed, options) -> SolveReport: supports -> seeded lifting
-> mixed cells -> binomial starts -> one tracked path per unit of mixed
volume; re-lifts once with seed+1 on GenericityError before surfacing
it. SolveReport keeps every raw path and offers solutions() /
distinct_solutions(tol) as views.
Validation: linear pair vs direct LU solve (1e-10), {x^2=2, y^2=3} smoke
(4 paths, +-sqrt2 x +-sqrt3), trinomial pair vs its eliminated quadratic
through the univariate solvers, dense conic pair (4 distinct verified
roots), proportional-equations failure honesty (no panic, non-Converged),
and bit-for-bit same-seed determinism.
Also: examples/homotopy.rs (end-to-end runnable proof), two bench rows
(tracker step ~1.8us, full trinomial solve ~65us), README "Polyhedral
homotopy (experimental)" section replacing the Roadmap future-work
framing, HOMOTOPY.md status appendix, crate-doc feature note.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_012bGPMXfU7c6jnHVSX33DJW
…nostics Predictors: TrackOptions::predictor selects Euler / Rk2 (midpoint) / Rk4 (classical), each stage a fresh Jacobian + LU tangent solve of the Davidenko ODE. Rk4 is the benchmarked default (~20x fewer steps, ~8x less wall time than Euler on the conic pair and cyclic-3; table recorded on Predictor::default). TrackOptions stays exhaustive on purpose (documented): pre-1.0 field additions are an accepted breaking change and callers keep `..Default::default()` construction. Step control: accepted steps now adapt dt by observed Newton effort (1 iteration -> *grow, 2 -> hold, converged on the max_newton-th -> *0.8); rejections still halve and the exact t=1 landing clamp stays. Documented on TrackOptions, including the coupling to predictor order; regression test pins total conic steps at <= the Phase 8 rule's measured 419 (now 186). The γ-twist (the headline finding): cyclic-3 choked — the textbook coefficient paths c·t^e never leave the real slice, so its real symmetric (maximally non-generic) coefficients made all six paths fold pairwise on the discriminant at one interior t (conjugate collisions, pivot ratios ~1e-8) for every seed. Fix: rotate every non-edge term by the endpoint-preserving phase exp(iγe(1−t)) — a homotopy-level gamma trick (γ = ln 2 fixed for reproducibility; CellHomotopy::with_gamma for explicit control, γ = 0 = textbook). Both endpoints stay bit-for-bit; cyclic-3 now tracks 6/6 on every seed and predictor. Diagnostics: Lu::pivot_ratio() (min/max pivot-norm ratio; documented as a cheap singularity-proximity signal, not a condition number) surfaces as PathResult::pivot_ratio from the final polished Newton solve; SolveReport gains converged_count(), failed_paths(), real_solutions(tol) (filters, never zeroes imaginary parts) and allocation-free Display impls for PathStatus and SolveReport; examples/homotopy.rs prints the report Display. Validation: cyclic-3 end-to-end with the structural oracle (permutations of (1, ω, ω̄): |coords| = 1, sum = 0, product = 1, MV = 6 = 6 distinct verified roots); the trinomial pair end-to-end over Complex<f32>; a proptest lane (tests/homotopy_prop.rs, 32 cases) with random complex coefficients on the fixed trinomial supports. Bench gains a cyclic-3 solve() row and the Euler/RK2/RK4 comparison table (wall time + steps). README and HOMOTOPY.md updated. 218 tests green on nightly, stable, and specialization; fmt/clippy clean on both toolchains; no-default+libm builds (host + thumbv7em) and the codegen guard unchanged. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_012bGPMXfU7c6jnHVSX33DJW
Tier A — examples/bench_suite.rs (also a test target): runs cyclic-3/4/5, katsura-3/4, noon-3, eco-4/5 and the trinomial/conic calibration rows end-to-end, printing mixed volume, cell/path counts, converged paths, offline vs tracking wall time, us/path, verified max residual, and honest failure tallies (--csv for machine-readable output). All 205 tracked paths across the 9 zero-dimensional systems converge with residuals <= 6e-15; cyclic-4 (proven positive-dimensional) fails all 16 paths honestly and stays in the table as a negative control. A pre-run gate excludes systems whose naive Pi C(|A_i|,2) cell enumeration projects past 60 s (cyclic-7: ~8.6e7 tuples, ~294 s) — the enumeration frontier, documented as a finding. Tier C — tools/oracle-sympy: an exact sympy oracle (Groebner over Q/Q(i)) recording quotient dimensions, Seidenberg radicality, and torus counts for every suite system. It settled: katsura-n (n+1-unknown convention) has 2^n distinct roots of which exactly MV lie on the torus (katsura-3: 6 of 8, katsura-4: 12 of 16 — matching the converged path counts exactly); cyclic-4 is positive-dimensional (two curves, verified by substitution); eco-4/5 have 4/8 roots, all on the torus. sympy's solve_poly_system was caught silently dropping 6 of katsura-3's 8 roots, so all counts come from quotient dimensions, never from solve(). Bug found by the suite and fixed: mixed_cells accepted exactly integer-singular candidate edge tuples (f64 LU leaves a ~1e-16 pivot where the exact pivot is 0), producing garbage normals that falsely tripped the genericity check on every katsura seed. Tuple singularity is now decided exactly over Z via Smith normal form, with a katsura-3 regression test. Tier B — BENCHMARKS.md (methodology incl. the fairness asymmetry stated both ways, recorded results, findings, literature context without any restated numbers) and tools/bench-external/hc_bench.jl + README: a same-schema HomotopyContinuation.jl harness for unrestricted machines (the dev container's proxy 403-blocks the Julia CDN, so it is untested here by construction). Tests: 220 per lane (was 218), all three lanes green; fmt/clippy (-D warnings, stable + nightly, --all-targets), no-default+libm (host + thumbv7em), doc, codegen guard, demo/homotopy/bench examples all clean. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_012bGPMXfU7c6jnHVSX33DJW
Replace the naive Pi C(|A_i|,2) tuple scan with an LP-pruned depth-first
tree search (a simplified variant of DEMiCs, Mizutani-Takeda-Kojima, DCG
2007), moving the polyhedral solver's enumeration frontier from cyclic-5
to cyclic-7 end-to-end and cyclic-8 enumeration-only.
- src/solvers/homotopy/lp.rs (new, module-private): a dense two-phase
simplex under Bland's rule solving the maximize-delta feasibility LP
{A_eq*alpha = b_eq, A_in*alpha - delta*1 >= b_in}; delta* > band means
strictly feasible, below -band prunes, inside the 1e-9 band surfaces
GenericityError (same re-lift contract as before). Iteration-cap safety
valve (Stalled) is never read as a pruning verdict; tolerances are
documented and sit far below the genericity band.
- src/solvers/homotopy/cells.rs: mixed_cells is now the tree search
(per-support edge pre-filter, static ascending-viable-count support
ordering, per-node feasibility LPs); full-depth tuples go through the
exact shared decision path (Z-SNF singularity gate, LU normal solve,
relative-band minimality check), so accepted cells match the naive
enumerator bit for bit. The old scan survives as mixed_cells_naive,
the documented reference implementation.
- tests/cells_oracle.rs (new): asserts identical cell sets (count, edge
tuples, edge matrices, exact volumes, normals to 1e-12) between both
enumerators on trinomial/conic, cyclic-3/4/5, katsura-3/4, noon-3,
eco-4/5, and seeded randomized supports; degenerate-lifting and
Phase-10 singular-tuple regressions pass unchanged.
- examples/bench_suite.rs: cyclic-6/7 rows with published mixed volumes
asserted (156, 924) - all 156/924 paths converge; the naive-projection
exclusion gate is gone; new --enum mode measures naive vs DEMiCs side
by side out to cyclic-8 (MV 2560 asserted, 10.4 s vs ~19 h projected).
- BENCHMARKS.md / HOMOTOPY.md: before/after tables (katsura-4 offline
113 ms -> ~8 ms, cyclic-6 enumeration 558 ms -> 36 ms, cyclic-7
~294 s projected -> 0.70 s measured), frontier section rewritten,
Phase 11 status recorded.
234 tests across 3 lanes; fmt/clippy (nightly + stable), no_std libm
(host + thumbv7em), doc, and the codegen guard all clean.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_012bGPMXfU7c6jnHVSX33DJW
Paths ending at singular isolated roots no longer die MinStepReached:
inside the endgame zone (t >= t_endgame, default 0.99) a min-step
failure — or an accepted corrector whose pivot ratio has collapsed
below endgame_pivot_threshold while dt sits below endgame_dt_threshold,
including the terminal accept at t = 1 where the ~sqrt(eps)-wide
cancellation zone of H around a multiple root lets Newton "converge" at
a point far less accurate than newton_tol suggests — hands the path to
a Cauchy endgame that computes the endpoint and its winding number.
- Complex-t evaluation on CellHomotopy: write_system_at_ct / eval_ct /
eval_jacobian_ct / dt_ct continue c·t^e through the principal branch
(single-valued for Re t > 0; the circle |1-t| <= r < 1 stays there)
and the gamma-twist exp(i·gamma·e·(1-t)) verbatim (entire in t). The
real-t methods remain the tracker's hot path (one real powf per term
beats complex ln/exp; endpoint exactness is easiest there); axis
agreement is pinned to ~1e-14 and dt_ct against complex central
differences in both the real and imaginary directions.
- The endgame (track.rs, private Endgame): walk-out to
r = max(1 - t_entry, endgame_radius) by radius-doubling Newton hops
at real t (at the raw stall radius ~1e-13 the Newton noise floor
eps/sigma_min exceeds the closure tolerance for windings >= 3), then
circle tracking of t(theta) = 1 - r·e^{i·theta} with an
Euler-in-theta predictor (dt/dtheta = -i·r·e^{i·theta}, sign pinned
by unit test — a sign error makes closure impossible) and a Newton
corrector at fixed complex t, bisecting failed theta-steps up to four
times. Closure after each full loop finds the winding (<= 8 default);
the endpoint is the mean of the equally spaced samples over the
closed cycle (trapezoid on a periodic function = the Cauchy integral;
kills every fractional Puiseux power exactly), gated by
||H(y,1)||_inf <= sqrt(newton_tol)·max(1, ||y||_inf) — singular roots
cannot be Newton-polished to the regular tolerance.
- New statuses ConvergedSingular { winding } / EndgameFailed (Display:
conv-singular / endgame-fail, width-aware); PathResult grows an
endgame_entered flag so tests can *prove* healthy paths never enter.
- Report surface: SolveReport::singular_count(), multiplicity_of(),
solutions() documented to include gated singular endpoints,
failed_paths() excludes them, Display gains a (+n singular) tally.
- Validation: double root {x²-2x+1, y-x} — both paths winding 2,
endpoints ~1e-15 from (1,1) (vs ~1e-8 attainable by plain Newton);
triple root {(x-1)³, y-1} — three paths winding 3, ~2e-12; the same
supports with simple roots stay endgame-free; f32 double root with
f32-scaled trigger constants — winding 2, ~1e-7. The full regular
suite (trinomial, conic, cyclic-3/5/6/7, katsura, noon, eco — 1000+
paths) is asserted endgame-free, in cargo test and at release
bench_suite scale.
- cyclic-4, the honest caveat: the endgame closes its 16 paths pairwise
at winding 2 on points that genuinely lie on the positive-dimensional
curves ((a,b,-a,-b), ab = ±1 — structurally verified, residuals
~1e-15). Nothing is fabricated, but winding certifies local branch
structure, not isolatedness: a winding-2 landing on a curve is
locally indistinguishable from an isolated double root, and telling
them apart needs witness sets (deferred, documented on the status,
in the pinned test, and in HOMOTOPY.md/BENCHMARKS.md).
- bench_suite: dbl-root row (endgame cost visibility), root-bearing
conv column, endgame-fail tally, per-row no-endgame assertions on
regular systems; README/HOMOTOPY.md/BENCHMARKS.md updated.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_012bGPMXfU7c6jnHVSX33DJW
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
This is a major architectural refactor of talrost's numerical tower and algebraic abstractions, fixing three critical correctness bugs and establishing the foundation for a polyhedral homotopy solver. The crate is reorganized around a proper trait hierarchy (
Monoid→Group→Semiring→Ring→Field), with newRealandScalartraits replacing the invertedFloat/Numberdesign. All container types (Matrix,Vector,Polynomial) are now generic over the algebraic tower, enabling integer matrices and dual-number automatic differentiation as first-class citizens.Key Changes
Critical bug fixes:
Matrixkernel dispatch using bitwise&instead of&&(e.g.,M == 2 & N & OnowM == 2 && N == 2 && O == 2), which silently corrupted results for non-matching dimensionsComplex::magnitude()returning|z|⁴instead of|z|(was.powi(2).powi(2), now.sqrt())Complex::normalize()and related operations that depended on the broken magnitudeNumerical tower redesign:
Realtrait for IEEE-754 ordered fields (f32,f64), replacing the invertedFloattraitScalartrait for fields with real-valued norms (implemented byRealtypes andComplex<F>)Number,Natural, andFloattraits; replaced with proper algebraic hierarchyComplex<F>now requiresF: Real(notFloat) and no longer implementsPartialOrd(mathematically correct)Complexto maintain accuracy at extreme magnitudes (1e±300)Container generalization:
Matrix<T, M, N>now requires onlyT: Ringfor structural operations (add/sub/mul/transpose),T: Scalarfor numerics (determinant/inverse/solve)Vector<T, N>similarly split:Ringfor algebra,Scalarfor normsPolynomial<T, N>generalized; root finders remainReal-boundMatrixstorage from column-major[[T; M]; N]to row-major[[T; N]; M](conventional)Index<(usize, usize)>andIndexMutfor matrix element accessNew modules:
src/real.rs: TheRealtrait with IEEE-754 constants and transcendental functionssrc/scalar.rs: TheScalartrait unifying real and complex numeric operationssrc/dual.rs: Dual numbers (Dual<T>andDualN<T, K>) for forward-mode automatic differentiationsrc/mvpoly.rs: Multivariate sparse polynomials (MPoly,MSystem) with monomial supportsrc/lattice.rs: Smith and Hermite normal forms for integer matricessrc/matrix/lu.rs: LU factorization with partial pivoting (extracted from matrix.rs)src/matrix/kernels.rs: Specialized multiply kernels (Strassen 2×2, Laderman 3×3, AlphaTensor 4×4) gated bymin_specializationsrc/roots.rs:Rootstype for bounded root collections (replaces NaN sentinels)Documentation and testing:
REVIEW.md: Architecture review identifying defects and design decisionsHOMOTOPY.md: Implementation notes for the homotopy solver pipelinetests/oracle.rs: Property-based tests againstnalgebraandnum-complextests/laws.rs: Algebraic law verification for the trait towerexamples/bench.rs: Micro-benchmarks for matrix operationstools/check_codegen.sh: Assembly-level codegen guard ensuring hot paths inline**Feature and toolchain
https://claude.ai/code/session_012bGPMXfU7c6jnHVSX33DJW