Part 4 sessions 4-5: symbolic algebra, primes, and number theory - #4
Merged
Conversation
Add exact/symbolic.rs: expression trees with a precedence-climbing parser, Display and LaTeX rendering, evaluation, exact symbolic differentiation, simplification, expansion, substitution, polynomial extraction, Taylor series, a stack-machine compiler, table-driven integration, numeric limits, root finding and critical points, gradients and Hessians. Simplification is a normaliser, not a prover. It folds constants, flattens nested sums and products, collects like terms by their non-numeric part, groups repeated bases into powers, and applies the power, exponential and logarithm identities. That is enough for the roadmap's headline property to fall out of arithmetic rather than a special case: differentiating sin(x)^2 + cos(x)^2 produces +2*cos*sin and -2*cos*sin, and the two collect to exactly zero. Two rules were missing until the tests found them. split_coeff did not recognise a negation, so its -1 stayed an opaque factor and x - x hashed under two different keys instead of cancelling. And a product of exponentials did not combine its arguments, so exp(x)*exp(-x) could not reach 1; exponential factors are now gathered into a single argument sum, which is the companion to the ln(exp(x)) rule that was already there. Sums and products need opposite operand orders. Constants sort last in a sum so a polynomial reads x^2 - 1, and first in a product so a term reads 5*x. Also backfills atan, sinh, cosh and tanh on core::dual::Dual. They were absent, which blocked cross-checking the symbolic derivative against forward-mode automatic differentiation for those functions -- the roadmap's stated property for this session. The new rules are tested against their closed forms and against cosh^2 - sinh^2 = 1. Verified by extracting the staged tree into a clean checkout: 2935 lib tests, 107 property tests, and clippy --all-targets -D warnings pass there, and the committed tree hash matches the one tested. Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Add discrete/primes.rs: three sieves that cross-check each other, deterministic Miller-Rabin for every u64, BPSW for BigInt, Pollard rho and p-1, Fermat's method, complete factorization, and prime counting. prime_count_meissel uses the Lucy_Hedgehog recurrence over the distinct values of n/i rather than Meissel's own formula. The state holds O(sqrt n) partial counts and each prime up to sqrt(n) sieves all of them at once, which is what the name promises -- pi(n) without a sieve to n. pi(10^9) is 50847534 in 87 ms in a debug build, so the roadmap's headline value is an ordinary test rather than something skipped for cost. The Lucas half of BPSW needed testing on its own terms. is_prime_bigint short-circuits below 2^62 and otherwise runs Miller-Rabin on random bases first, which rejects essentially every composite before Lucas is reached: replacing the whole strong Lucas test with `return true` left the entire suite green. It is now tested directly against the odd numbers below 20000, where it must accept every prime and exactly the five strong Lucas pseudoprimes for Selfridge's parameters. The implementation reproduces that set, 5459, 5777, 10877, 16109 and 18971, which is a much stronger statement than not crashing. The test also pins the complementarity BPSW rests on: none of the six base-2 strong pseudoprimes below 20000 is a Lucas pseudoprime, so each test catches what the other misses. Factorization is checked by reconstruction on 2000 random values up to 10^12 plus the shapes that defeat a single method: a semiprime of two near-equal factors, prime powers, and a prime beyond trial division. number_theory.rs is committed empty; it lands in its own commit. Verified by extracting the staged tree into a clean checkout: 2942 lib tests, 107 property tests, and clippy --all-targets -D warnings pass there, and the committed tree hash matches the one tested. Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Add discrete/number_theory.rs: modular arithmetic and the Chinese remainder theorem for general moduli, the multiplicative functions and their Dirichlet relations, multiplicative order and primitive roots, discrete logarithms by baby-step giant-step and Pohlig-Hellman, Legendre and Jacobi symbols, Tonelli-Shanks, the Carmichael function, digit and Collatz utilities, sums of two and four squares, primitive Pythagorean triples by the Berggren tree, Gaussian integer factorization, the Frobenius number, Egyptian fractions, Zeckendorf representations, Lucas sequences, Diophantine solving, and Stern-Brocot and Farey navigation. The tests assert the defining identities rather than sampled values: sum of phi over the divisors of n is n, sum of mobius is one exactly at n = 1, and mobius and phi are a Dirichlet inverse pair. Carmichael numbers below 10^4 are exactly the known seven and each is verified composite yet Fermat-pseudoprime to every coprime base. Tonelli-Shanks is checked by squaring the root back, and None is checked to mean a genuine non-residue rather than a failure to find one. The two discrete logarithm routines agree with each other and with brute force. Three notes on the specification. quadratic_diophantine_solve is read as a x^2 + b y^2 = c, so an indefinite form returns an empty vector rather than enumerating an infinite Pell family. frobenius_number returns zero when a unit coin is present and None when the coins share a factor. stern_brocot_nth indexes the tree breadth-first from the root 1/1. Verified by extracting the staged tree into a clean checkout: 2981 lib tests, 107 property tests, and clippy --all-targets -D warnings pass there, and the committed tree hash matches the one tested. Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
There was a problem hiding this comment.
🟡 Changes recommended
There are correctness/contract issues in the new discrete-math code paths (notably overflow in segmented sieving and factorization failure handling) that should be addressed before approval.
Once you've addressed the issues Copilot identified, you can request another Copilot review.
Pull request overview
This PR completes Phase A and starts Phase B of docs/ROADMAP_PART4.md by adding new discrete-math functionality (primes + number theory) and a new symbolic algebra subsystem, with tests that directly exercise previously untested rare branches (notably the Lucas half of BPSW).
Changes:
- Add
exact::symbolicexpression trees with parsing/printing, simplification, differentiation, Taylor series, compilation, and numeric limit/solve helpers. - Add
discrete::primes(sieves, deterministic u64 primality, BPSW forBigInt, factorization, prime counting) plusdiscrete::number_theory(CRT, arithmetic functions, discrete logs, residues, Diophantine/digit problems, etc.). - Extend
core::Dualwithatan/sinh/cosh/tanhto enable AD cross-checks for the new symbolic derivatives.
File summaries
| File | Description |
|---|---|
| src/lib.rs | Exposes the new discrete module at the crate root. |
| src/discrete/mod.rs | Adds the discrete-math module namespace. |
| src/discrete/primes.rs | Implements sieves, primality (u64 + BigInt/BPSW), factorization, and prime counting with dedicated Lucas-branch tests. |
| src/discrete/number_theory.rs | Implements CRT, arithmetic functions, discrete logs, residue tests/roots, and related number theory utilities with extensive tests. |
| src/exact/mod.rs | Exposes the new symbolic module and re-exports Expr. |
| src/exact/symbolic.rs | Adds the symbolic expression system (parse/print/simplify/diff/compile/etc.) and tests including AD cross-checks. |
| src/core/dual.rs | Adds missing transcendental/hyperbolic Dual operations plus derivative/identity tests. |
Review details
Suppressed comments (1)
src/discrete/primes.rs:536
factorize_bigintis documented as a complete factorization, but onpollard_rho_bigintfailure it pushesminto the output even though it is known composite (it already failedis_prime_bigint). This makes the API contract unclear and can also make tests/callers flaky if rho occasionally fails to split a composite.
/// The complete factorization of a `BigInt`.
///
/// # Panics
/// Panics if `n` is not positive.
#[must_use]
pub fn factorize_bigint(n: &BigInt, rng: &mut Rng) -> Vec<(BigInt, u32)> {
assert!(!n.is_negative() && !n.is_zero(), "factorization needs a positive integer");
- Files reviewed: 7/7 changed files
- Comments generated: 5
- Review effort level: Lite
💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.
Comment on lines
+48
to
+54
| // First multiple of p at or above lo, never below p^2. | ||
| let start = (lo.div_ceil(p) * p).max(p * p); | ||
| let mut m = start; | ||
| while m < hi { | ||
| is_p[(m - lo) as usize] = false; | ||
| m += p; | ||
| } |
Owner
Author
There was a problem hiding this comment.
Confirmed and fixed. lo.div_ceil(p) * p overflows at lo = u64::MAX for p = 7, 11 and 13, which in release wraps below lo and then underflows the m - lo index. It now rounds up by the remainder with checked_add, so every intermediate stays at or below the answer and a first multiple that does not fit means the window holds no multiple at all.
One note on testing it: the overflow is not reachable through sieve_segmented itself, because a window that high needs a base sieve to sqrt(u64::MAX) — four billion booleans. So the regression test pins the arithmetic directly, against the naive form wherever the naive form is valid and against the top of the range where it is not.
Generated by Claude Code
Comment on lines
+510
to
+516
| match pollard_rho(m) { | ||
| Some(d) => { | ||
| stack.push(d); | ||
| stack.push(m / d); | ||
| } | ||
| None => found.push(m), | ||
| } |
Owner
Author
There was a problem hiding this comment.
Fixed for factorize, and the u64 contract is now provable: trial division to the square root is added as a guaranteed-terminating fallback when Pollard rho gives up. (Rho splits every composite below two million, so nothing observed reaches it, but the path existed.)
factorize_bigint at line 530 is different and I have not made it unconditional. Splitting a large composite has no cheap certain method, so there is no BigInt analogue of the trial-division fallback. It now tries rho three times — forty-eight independent polynomials — and the doc says plainly that an unsplit cofactor comes back as a single entry despite being known composite, and that a caller needing certainty should test each returned base with is_prime_bigint. Claiming completeness there would be false, so the contract was narrowed rather than the code made to promise something it cannot.
Generated by Claude Code
Comment on lines
+1357
to
+1365
| if q == 1 { | ||
| return Rational::from_i64(p * order + 1, order); | ||
| } | ||
| let inv = mod_inverse_u64(p.rem_euclid(q) as u64, q as u64) | ||
| .expect("a reduced fraction has coprime parts"); | ||
| let s0 = ((q as u64 - inv % q as u64) % q as u64) as i64; | ||
| let s = s0 + q * ((order - s0) / q); | ||
| let r = (1 + p * s) / q; | ||
| Rational::from_i64(r, s) |
Comment on lines
+161
to
+176
| pub fn variables(&self) -> Vec<String> { | ||
| let mut out = Vec::new(); | ||
| fn walk(e: &Expr, out: &mut Vec<String>) { | ||
| if let Expr::Var(n) = e { | ||
| if !out.contains(n) { | ||
| out.push(n.clone()); | ||
| } | ||
| } | ||
| for c in e.children() { | ||
| walk(c, out); | ||
| } | ||
| } | ||
| walk(self, &mut out); | ||
| out.sort(); | ||
| out | ||
| } |
The Verify workflow has never run. Every one of its 18 runs since it was
added failed instantly with zero jobs, which is what a workflow that fails
to parse looks like from the API.
Two defects:
1. The Miri step's command was unquoted:
run: cargo miri test --lib -- core:: linalg:: spatial::
`core:: linalg::` is a colon followed by a space, so YAML parses it as a
nested mapping and rejects the whole file. Locally:
mapping values are not allowed here
in ".github/workflows/verify.yml", line 50, column 44
Every `run:` in the file is now quoted so a command containing a colon
cannot recur as a parse error.
2. The clippy step combined `-D warnings` with `-W clippy::float_cmp` and
three other float-accuracy lints, intending them as advisory. `-D
warnings` denies the whole warning level, so naming an allow-by-default
lint with `-W` promotes it to an error rather than softening it. That
combination produces 3429 errors, mostly suboptimal_flops. Rewriting
those expressions as mul_add changes rounding, so they need review one
at a time rather than a blanket gate; the flags are removed and the
reasoning recorded in the file.
Miri is also narrowed from three module filters to `core::` (27 tests).
Miri interprets at roughly a hundredth of native speed and the crate has no
`unsafe`, so it is a backstop, not the primary check.
Verified locally: the file parses, `cargo test --release --test properties`
passes 107 tests, and `cargo clippy --all-targets -- -D warnings` finishes
clean. The kani and miri jobs remain unverified — neither tool is installed
here, so CI will be their first real execution.
An earlier commit message claimed the strict lint job was "green rather
than red on arrival". The lint command did pass locally, but the workflow
containing it never parsed, so no job ever ran.
Co-Authored-By: Claude <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Add discrete/combinatorics.rs, partitions.rs, sequences.rs and disjoint_set.rs. combinatorics.rs covers binomials and multinomials, permutation and combination enumeration, the permutation group, the classical counting numbers, Burnside and Polya, and the named puzzles. Most tests compare a formula against exhaustive enumeration of the objects it counts rather than against a published table: Eulerian numbers against permutations sorted by ascent count, Narayana numbers against Dyck paths sorted by peak count, ballot numbers against every vote sequence, necklaces and bracelets against brute-force orbits under rotation and reflection, and all twelve entries of the twelvefold way against enumerated maps quotiented four ways. The Hanoi move list is simulated on three real pegs, the shuffle order is checked by shuffling until the deck returns, and the magic squares are checked on every row, column and diagonal in all three residue classes. partitions.rs covers the pentagonal recurrence, Young diagrams and hook lengths, and RSK. The hook length formula is checked against a direct fill of the diagram; RSK is checked to be injective on S_n, to satisfy Schuetzenberger's theorem that inverting the permutation swaps the two tableaux, and to satisfy Schensted's, that the first row is the longest increasing subsequence and the row count the longest decreasing one. sequences.rs covers Taylor coefficients by Cauchy's integral, linear recurrences by iteration and by matrix power, Berlekamp-Massey over Q and over GF(2), and the named integer sequences. The recovered recurrences are checked by regenerating the input rather than by comparing coefficients, and the OGF error is checked against the exact aliasing term r^N/(1-r^N) rather than a tolerance. Four defects the tests found while writing them. nth_permutation read the factoradic digits from the wrong end: it divided by the position radix and used the remainder, which walks the positions in reverse. It now divides by (n-1-j)! at position j. The restricted growth string successor never reset the suffix, so set_partitions_iter emitted one string per n rather than Bell(n). binomial_u64 reported overflow for results that fit. The running product before the division is C(n,k+1)*(k+1), up to k times the answer, so a u64 accumulator overflows first. It accumulates in u128 and tests the coefficient itself against u64::MAX, which makes None mean exactly "does not fit". random_permutation drew its swap index with next_u64() % (i+1). The generator is an LCG modulo 2^64, where bit b has period 2^(b+1), so a small modulus reads the shortest-period bits: the lowest merely alternates. Shuffling six elements that way cycles through a handful of arrangements instead of sampling the 720. It now takes the high half of a widening multiply. The test that caught this counts how often each symbol lands in each position over 20000 draws, which a validity-only check would have passed. Also fixes seven clippy errors that only the current stable raises, none of which had ever been seen: verify.yml has never had a passing run, and this environment's toolchain was four releases behind. Substantive rather than suppressed -- an explicit counter loop becomes a range, a comparator becomes sort_by_key, a manual checked division becomes checked_div, chunks_exact(2) becomes as_chunks, a loop becomes while let, and two no-op .max(0) calls on unsigned values are dropped. One of those two was hiding a real fault: additive_evolving computed (env.len() - 1) on an envelope it never checked was non-empty, which underflows. Empty envelopes now contribute nothing. Verified by extracting the staged tree into a clean checkout: 3065 lib tests, 117 property tests, and clippy --all-targets -D warnings pass there under rustc 1.98, the same version CI uses. Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
The Verify workflow ran for the first time on the previous commit. Two of
its four jobs still failed; this addresses Miri.
Miri reported five failures in core::dual and core::interval, all of them
exact-value float assertions:
assertion `left == right` failed: p' at -3
left: -70.99999999999997
right: -71.0
None of them is a defect. Miri does not call the host's sin, exp or powi;
it evaluates them itself, and is deliberately non-deterministic within the
slack the language allows for those operations. A test asserting an exact
double therefore fails under Miri whatever the code does. The five are
marked #[cfg_attr(miri, ignore)] with that reason recorded at each site.
They are unaffected everywhere else -- core:: still runs 27 tests, none
ignored, under a normal cargo test.
Also bounds the two slow jobs. Kani had been running for over fifty
minutes and was still going when this was written: the harnesses quantify
over whole f64 domains and several take a square root, which the solver
bit-blasts. With no bound a single slow harness holds a runner until the
six-hour job default. 90 minutes for Kani and 45 for Miri turn that into a
failure that names the problem instead. Miri's own 27 tests took twelve
minutes, which is recorded in the file so the scope choice is legible.
Co-Authored-By: Claude <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Add graph/core.rs and graph/paths.rs. core.rs holds the representation and the structural queries: connected and strongly connected components by an iterative Tarjan, condensation, bipartiteness, topological order, bridges and articulation points, Eulerian circuits and paths by Hierholzer, Hamiltonian paths by bitmask dynamic programming, girth, the distance metrics, clustering and transitivity, k-cores, the named and random generators, line and product graphs, canonical forms and isomorphism, graph6, and the matrix-tree theorem over the integers. paths.rs holds Dijkstra, Bellman-Ford, Floyd-Warshall, Johnson, A*, bidirectional search, Yen's k shortest paths, widest and minimax paths, the DAG routines, three minimum spanning tree algorithms, the second-best spanning tree, Dreyfus-Wagner Steiner trees, Held-Karp, the tour heuristics, Christofides, and the Chinese postman. The tests compare each algorithm against the definition it implements rather than against stored answers. Strongly connected components are checked against mutual reachability computed by transitive closure, one pair at a time. Bridges and articulation points are checked by actually removing each edge or vertex and recounting the components. The girth is checked against exhaustive cycle search, the k-core against direct peeling, and Hamiltonian path existence against every permutation. Eulerian walks are checked to consume each edge exactly once with a tally, not merely to have the right length. Held-Karp is checked against every tour, Christofides against Held-Karp for its 1.5 bound, and the Steiner tree against a brute force over every subset of Steiner points. Cayley's formula is checked to n = 14, where n^(n-2) is past 2^53 and an f64 determinant could not be exact. One defect the tests found. The heap key ordered `f64` with `partial_cmp(..).unwrap_or(Equal)`, which is the obvious thing to write and is wrong: it makes a NaN key compare equal to every other key, so the ordering is not transitive, `Ord`'s contract is broken, and the heap can return items out of order. A NaN pushed among 1, 2 and 3 came back second. `total_cmp` is a genuine total order and puts a positive NaN above infinity, so a NaN weight now settles last instead of corrupting the search. Every weight comparison in the module uses it. Also corrects the Chinese postman's documented contract: an edgeless graph is disconnected but has nothing to cross, so it returns the empty route rather than failing. Only edges spanning more than one component make the problem unsolvable. Verified by extracting the staged tree into a clean checkout: 3109 lib tests, 125 property tests, and clippy --all-targets -D warnings pass there under rustc 1.98, and the committed tree hash matches the one tested. Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Five findings from the automated review on PR #4. Four are real; the fifth is a contract that could not be made unconditional, and is now documented accurately instead. sieve_segmented computed the first multiple of p at or above lo as `lo.div_ceil(p) * p`. That rounds lo up past u64::MAX when lo is within p of the top, which wraps to a value below lo and then underflows the `m - lo` index. Confirmed directly: it overflows at lo = u64::MAX for p = 7, 11 and 13. It is not reachable through sieve_segmented itself, because a window that high needs a base sieve to sqrt(u64::MAX), four billion booleans, so no test can drive it; the arithmetic is pinned directly instead. Rounding up by the remainder keeps every intermediate at or below the answer, and a first multiple that still does not fit means the window holds no multiple at all. The `p * p >= hi` guard is also written with checked_mul now -- it cannot overflow for any prime the base sieve produces, and saying so beats leaving the reader to check. factorize is documented as a complete prime factorization but pushed the cofactor into the result when Pollard rho gave up, even though it had already failed is_prime_u64 and was therefore known composite. Rho splits every composite below two million, so nothing observed reaches it, but the path existed. Trial division to the square root is added as a guaranteed-terminating fallback, which makes the guarantee provable rather than overwhelmingly likely. factorize_bigint has the same shape and no such fallback: splitting a large composite has no cheap certain method. It now tries rho three times, forty-eight independent polynomials in total, and the doc says plainly that an unsplit cofactor comes back as a single entry and that a caller needing certainty should test each base with is_prime_bigint. Claiming completeness there would be false. farey_next formed `p * order` and `p * s` as i64 products, which wrap silently in release. They are formed in i128 and converted back with a checked cast, so an out-of-range successor panics with a message rather than returning a wrong fraction. Expr::variables deduplicated by scanning a growing vector, which is O(v^2) string comparisons in the number of distinct variables. It collects into a BTreeSet, which also supplies the sort. Also removes a dead `a = 0; let _ = a;` at the end of jacobi_bigint and the `mut` it existed to justify. Each fix carries a regression test that pins the behaviour: the first-multiple arithmetic against the naive form wherever the naive form is valid and against the top of the range where it is not, factorize against primality of every returned base on prime powers and near-equal semiprimes, farey_next against the neighbour identity r*q - p*s = 1 across F_1 to F_40 plus a should_panic for the overflow, and variables against a 300-variable expression repeated twice. 3114 lib tests and 125 property tests pass, and clippy --all-targets -D warnings is clean under rustc 1.98. Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
lu_decompose_3x3_finite_or_err carried #[kani::unwind(6)] on a loop that
runs nine times, so the unwinding assertion could never discharge:
Check 655: verification::linalg::lu_decompose_3x3_finite_or_err.unwind.0
- Status: FAILURE
- Description: "unwinding assertion loop 0"
That is not a property violation. Kani reports it when the bound is too
small to cover the loop, and the bound must exceed the trip count, so nine
iterations need ten. At ten the harness verifies: 655 checks, none failed,
68 seconds.
The harness has been wrong since it was written. Nobody could have seen
it, because verify.yml never parsed and so the Kani job never ran once.
This was found by installing Kani locally and timing every harness rather
than pushing another guess at CI.
That timing sweep also shows three harnesses exceeding a five-minute
budget -- bisection_result_is_inside_bracket, interval_mul_contains_corner
_products and mat3_inverse_never_divides_by_zero -- against 25 to 32
seconds for the ones that finish. Those are handled separately once the
full table is in; they quantify over whole f64 domains with symbolic
multiplication, division, or a 34-iteration loop, which CBMC has to
bit-blast.
Co-Authored-By: Claude <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Add graph/flow.rs and graph/matching.rs.
flow.rs holds Dinic and push-relabel, minimum cut, the global minimum cut
by Stoer-Wagner, minimum-cost maximum flow, circulations with demands and
lower bounds, both forms of Menger's theorem, Gomory-Hu trees by
Gusfield's construction, and the maximum-weight closure with its project
selection specialisation. matching.rs holds Hopcroft-Karp, the Hungarian
algorithm, the auction algorithm, Edmonds' blossom algorithm,
Gale-Shapley, Irving's stable roommates, Konig vertex covers and Hall's
condition.
The tests check each result against the definition rather than against a
stored answer. Max-flow is compared with min-cut and with a second flow
algorithm that never holds a valid flow until it finishes, and the flow
itself is checked for conservation at every interior vertex and capacity
on every arc. Stoer-Wagner is compared with the best over all C(n,2) s-t
cuts. Menger's theorem is checked by removing edges and vertices and
recounting. The Gomory-Hu tree is checked to encode every pairwise
minimum cut as the lightest edge on its tree path. Blossom matching is
checked by Berge's lemma -- that no augmenting path remains -- which is a
different statement from the algorithm's own search.
Four defects the tests found.
blossom_max_matching contracted blossoms but never rewired the parent
pointers through them. That rewiring is the lifting step and is the whole
difficulty of Edmonds' algorithm: without it an augmenting path traces
back out of the tree by the wrong edge, and the result was an asymmetric
pairing with m[1] = 4 but m[4] != 1. Rewritten with a proper lowest
common ancestor over blossom bases and a path-marking pass that rewires
parents along the odd cycle.
stable_roommates never terminated. Irving's rotation elimination has
y_i reject x_{i+1} and everyone below, but x_{i+1} is *defined* as the
last entry of y_i's list, so a non-strict comparison rejects nobody: no
list shrinks, no rotation is consumed, and the loop spins forever. The
test hung rather than failed. With the strict comparison the same test
finishes in 0.01 seconds.
circulation_with_demands had the super-source and super-sink inverted. A
vertex that must receive was wired to the source rather than the sink, so
every feasible instance came back None.
vertex_disjoint_paths gave each original edge capacity n. Two vertex-
disjoint paths cannot share an edge either, since sharing one means
sharing both its endpoints, so a single edge between adjacent vertices
reported two paths instead of one.
The auction algorithm was correct but reset its prices between epsilon
rounds, which discards the entire mechanism of epsilon scaling: each
round became a fresh auction and the last one, at the smallest epsilon,
ran the slowest variant there is, on the order of n^2 (cost range) / eps
bids. Its test took minutes. Carrying the prices over brings it to 0.04
seconds.
Also drops three n == 0 guards that were dead code. Matrix::zeros asserts
positive dimensions, so an empty cost matrix cannot be constructed.
Verified by extracting the staged tree into a clean checkout: 3133 lib
tests, 132 property tests, and clippy --all-targets -D warnings pass
there under rustc 1.98, and the committed tree hash matches the one
tested.
Co-Authored-By: Claude <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
The Kani job had never finished. It ran for over 100 minutes on one PR
and was still going. Rather than guess at it again from CI, Kani was
installed locally and all twenty harnesses were timed individually
against a five-minute budget. The result splits cleanly along what each
harness asserts:
24s normalize_angle_never_panics
25s displacement_is_finite_on_bounded_inputs
25s mat3_identity_inverse_is_identity
25s mean_panics_on_empty
25s projectile_range_panics_on_nonpositive_g
25s vec3_dot_with_self_nonnegative
27s kinetic_energy_nonnegative_for_nonnegative_mass
28s factorial_is_monotone_and_finite_below_171
28s projectile_range_never_panics_with_positive_g
32s escape_velocity_finite_nonnegative
34s shannon_entropy_never_panics_on_nonempty
52s vec3_normalized_never_produces_nan
68s lu_decompose_3x3_finite_or_err
---
>300s bisection_result_is_inside_bracket
>300s interval_mul_contains_corner_products
>300s mat3_inverse_never_divides_by_zero
>300s mean_of_bounded_slice_is_bounded
>300s orbital_velocity_below_escape_velocity
>300s ray_aabb_interval_ordered
>300s variance_is_nonnegative
Everything above the line asserts panic-freedom, finiteness or a sign.
Everything below asserts a numeric relation between symbolic float
expressions. CBMC decides floating point by bit-blasting it into SAT:
proving a result is finite constrains few bits, while proving that one
symbolic product or quotient bounds another constrains the whole 53-bit
mantissa of every intermediate, and the instance stops being tractable.
That is the distinction the module's own documentation already drew --
transcendentals are modelled as unconstrained finite values, so those
harnesses prove panic-freedom rather than numeric bounds. The slow seven
ask for numeric bounds anyway.
They are kept rather than deleted: each states something true, and a
future Kani or solver may decide them. They now sit behind a `kani-slow`
cargo feature, off by default, and are runnable with
cargo kani --features kani-slow --harness <name>
`cargo kani list` reports 13 harnesses by default and 20 with the
feature. The full default run was executed locally on Kani 0.67, the
version the action installs: 13 successfully verified, 0 failures, 67
seconds wall-clock with -j. The workflow now passes -j and its timeout
drops from 90 minutes to 45.
The measured table is recorded in src/verification/mod.rs so the split is
legible and revisitable rather than folklore.
Co-Authored-By: Claude <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Part 4 session 9, first half: src/graph/spectral.rs.
Laplacians (combinatorial and normalized) and their spectra, algebraic
connectivity and the Fiedler vector, spectral bisection and k-way
spectral clustering, Kirchhoff's spanning-tree count, PageRank, HITS,
eigenvector and Katz centrality, Brandes betweenness, closeness and
harmonic centrality, effective resistance and commute time, random-walk
stationary distributions and mixing time, Cheeger bounds and expander
testing, graph energy, the Estrada index, cospectrality, modularity,
Louvain, and label propagation.
Twenty tests, each against a closed form, an independent algorithm, or
the definition: the Laplacian's zero-eigenvalue multiplicity against a
component count, the spectra of K_n, C_n, P_n, the star and the
hypercube against their formulas, the normalized spectrum reaching two
exactly on a bipartite component, the matrix-tree count against exact
enumeration, betweenness of a star centre against (n-1)(n-2)/2, Foster's
theorem, commute time against 2m R, Cheeger's inequality against
brute-forced conductance, the K_{1,4} / C_4+K_1 cospectral pair, and
planted-partition recovery.
Three defects the tests found, all in the implementation:
- eigenvector_centrality did not converge on a bipartite graph. The
adjacency spectrum is symmetric about zero there, so the extreme
eigenvalues tie in magnitude and power iteration flips between their
eigenvectors forever, returning whichever phase the loop bound
happened to stop in. Iterating on A + cI for a Gershgorin bound c
breaks the tie and moves no eigenvector.
- mixing_time_estimate filtered every transition eigenvalue of
magnitude one out of its gap computation, which discarded exactly
the evidence that the walk does not converge. It reported a finite
mixing time for bipartite and disconnected graphs. It now drops one
zero -- the Perron eigenvalue -- and takes the largest magnitude
among the rest.
- modularity counted degrees with weighted_degrees, which drops
self-loops because the Laplacian does. Modularity does not: a loop
adds two to a vertex's degree with nothing to cancel it. Louvain's
contraction turns each community's internal edges into exactly such
a loop, so the total edge weight shrank at every level of the
recursion and the gain formula compared against the wrong 2m. On
four cliques joined in a path it merged everything into one
community and scored modularity zero. Now pinned by a test that
contracts a random partition and requires modularity to be
unchanged.
Co-Authored-By: Claude <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Part 4 session 9, second half: src/graph/coloring.rs.
Greedy colouring in four orders (natural, largest-first, degeneracy,
DSATUR), Welsh-Powell and its bound, the exact chromatic number by
bracket-and-search, the chromatic polynomial by memoised
deletion-contraction, Vizing edge colouring by the Misra-Gries
construction, a time-limited constraint search for k-colourability,
optimal interval-graph colouring, map colouring from adjacency lists,
Bron-Kerbosch with pivoting for maximal and maximum cliques, greedy and
exact independent sets, the vertex-cover two-approximation and its exact
counterpart, a greedy dominating set, and the Eades-Lin-Smyth feedback
arc set.
Eleven tests, each against a closed form, an independent algorithm, a
theorem, or exhaustive search:
- Welsh-Powell's colour-class sweep is asserted equal to largest-first
greedy, which is the theorem that the two procedures are one.
- The chromatic polynomial is checked to count what it claims: its
value at every k from zero to five equals the exhaustively counted
proper k-colourings, on random graphs, plus the closed forms for
cycles, trees and complete graphs, P(C5, 3) = 30, degree n, monic,
and an x^(n-1) coefficient of minus the edge count.
- Gallai's identity alpha + tau = n, the Caro-Wei bound on the greedy
independent set, the factor of two on the cover approximation, and
set cover's ln(n) + 1 on the dominating set.
- Vizing's bound both ways: never more than Delta + 1, never fewer
than Delta, and exactly Delta + 1 on every odd cycle. Separately
stressed over twenty thousand random graphs of up to thirty-one
vertices.
- The crown graph, where the natural order takes four colours on a
two-chromatic graph and DSATUR takes two, so the order really is the
algorithm.
Three defects the tests found:
- The Misra-Gries rotation chose its fan vertex without rechecking the
fan. Inverting the alternating path recolours an edge at the hub and
can occupy, at a fan vertex, the colour the fan property needed free
there, so a prefix that was a fan before the inversion need not be
one after. The prefix is now re-established first.
- The path inversion wrote each edge back into the incidence table as
it went. Two consecutive path edges meet at a vertex and exchange
colours there, so the second erased the entry the first had just
made and the table drifted out of step with the colouring, which
then went silently improper. Clearing the whole path before writing
any of it back fixes it.
- The interval sweep spent a colour on an empty interval. A half-open
interval whose ends coincide meets nothing, so charging it a colour
of its own pushed the total past the maximum overlap, which is the
one guarantee the sweep exists to make.
Co-Authored-By: Claude <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Part 4 session 9, third half: src/graph/layout.rs. Completes roadmap
item 6c.
Metric layouts (Kamada-Kawai, stress majorization by SMACOF from a
classical-scaling start, Fruchterman-Reingold, spectral), structural
layouts (circular, shell, Reingold-Tilford, Sugiyama), the stress
functional itself in two dimensions and in n, straight-line crossing
counting, biconnected decomposition, and planarity by Demoucron's path
addition -- which returns the faces, so the embedding comes with the
answer.
Thirteen tests. What each of them pins:
- Stress majorization is required to be monotone at every round on
every random graph, which is the only property that distinguishes
majorization from gradient descent on the same objective. Kamada-
Kawai is held to the weaker statement it can actually make: never
worse than the drawing it started from. That is not free either --
Newton on a non-convex energy steps uphill happily, so the step is
accepted only when the energy falls.
- A path laid out in one dimension must come out with consecutive
vertices exactly one apart, and a nine-cycle in two dimensions must
come out with every edge the same length.
- The spectral layout's two coordinates are checked to satisfy L x =
lambda x entry by entry for the right two eigenvalues, to be
centred, and to be orthogonal to each other.
- Reingold-Tilford is held to all three of its defining properties at
once -- depth is the height, nothing at a depth overlaps, every
parent is centred over its outermost children -- plus zero crossings
and, on a complete binary tree, exact symmetry and leaves packed at
exactly the separation.
- Crossing counts against the closed form: a complete graph drawn in
convex position has one crossing per four vertices, so K_n gives
n choose 4, checked for n from three to nine.
- Planarity against Kuratowski's graphs and their subdivisions, the
planar families, and the sharp boundary that K5 and K3,3 become
planar on the removal of any single edge. Then against an
independent algorithm: enumerate every rotation system, trace its
faces, and read off the genus. Demoucron and the genus computation
agree on every graph small enough to enumerate.
- The returned embedding is checked to be one: Euler's formula holds,
every face is a closed walk in the graph, and every edge borders
exactly two faces.
Co-Authored-By: Claude <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Five properties that hold between graph/spectral, graph/coloring,
graph/layout and graph/matching, which no one module's own tests can
state:
- Koenig's theorem, with the matching from the blossom algorithm and
the cover from maximum independent sets, on bipartite graphs; plus
the odd cycles that show why the theorem needs bipartiteness, where
the cover is exactly one larger than the matching.
- Hoffman below and Wilf above: the adjacency spectrum brackets the
chromatic number, computed by Jacobi rotations on one side and an
exhaustive colouring search on the other, with both bounds shown
tight on complete graphs.
- The four colour theorem against the planarity test: every graph
called planar must be four-colourable, must obey Euler's bound, and
must have degeneracy at most five, so the smallest-last order never
opens a sixth colour. Also asserts that some drawn planar graph
actually needed four, so the property is not passing vacuously.
- The chromatic number squeezed by the clique number below and by
chi times alpha at least n above.
- Every straight-line drawing of a graph the planarity test rejects
must cross, checked across four layout algorithms -- a planarity
test that wrongly said no would pass its own module's tests and fail
here -- and the tree layout draws every random tree with no
crossings at all.
Co-Authored-By: Claude <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Part 4 session 10, first half: src/codes/checksum.rs, a new codes/
module.
Parity, Fletcher-16 and -32, Adler-32, a parametric CRC covering every
named variant, CRC-32/CRC-16-CCITT/CRC-8 as named instances, the
reflected table form, Luhn, ISBN-10 and -13, Verhoeff, Damm, and Hamming
distance.
Ten tests. The point of a checksum is which errors it catches, so that
is what they assert rather than that the bytes come out the same twice:
- Eight named CRCs against their published check values, and the
table-driven CRC-32 against the bit-at-a-time one on random input.
- Every burst of w bits or fewer is detected by a CRC of width w, on
five parameter sets. This needed the bits numbered in the order the
CRC actually consumes them -- most-significant first within a byte,
or least-significant first when the CRC reflects its input -- since
numbering them the other way scatters a window across up to twice
its span and the theorem stops applying.
- A zero-seeded CRC is linear over GF(2), which is the fact that makes
the burst statement a statement about error patterns at all.
- A generator with an even number of terms is divisible by x + 1 and
so detects every odd number of bit errors. CRC-16/CCITT and
CRC-8/SMBUS qualify; the test records that CRC-32 has fifteen terms
and does not, rather than claiming a guarantee it lacks.
- Fletcher and Adler against transposition. Both weight byte m by the
number of bytes after it, so a swap moves the checksum by
(d_i - d_j)(i - j); the test asserts Fletcher-16 catches the swap
exactly when that survives its modulus of 255, and that Adler-32,
whose modulus is the prime 65521, can never be reached by a product
that small -- so it never misses one.
- Luhn catches every single-digit error and every adjacent
transposition except 09 against 90, which the test requires to
occur rather than working around.
- Verhoeff and Damm catch every single-digit error and every adjacent
transposition, with no exception, over thousands of cases.
- ISBN-10's prime modulus catches every transposition; ISBN-13's
composite one misses exactly those of digits differing by five, and
the test asserts both halves of that.
- Hamming distance satisfies the metric axioms and is translation
invariant, which is why a linear code's minimum distance is its
minimum non-zero weight.
Co-Authored-By: Claude <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Part 4 session 10, second half: src/codes/block.rs. Completes roadmap
item 7a.
A bit-packed GF(2) matrix with elimination, rank, solving and kernel
bases; LinearCode carrying both a generator and a parity check matrix;
the Hamming, extended Hamming, repetition, single parity check,
Golay(23), Golay(24) and Reed-Muller families; syndrome decoding and the
explicit standard array; weight enumerators over exact integers; duals;
Hamming(7,4) in its classical bit layout; the Singleton, sphere-packing,
Gilbert-Varshamov and Plotkin bounds; Gallager's regular LDPC
construction with belief-propagation and bit-flipping decoders.
Nine tests:
- GF(2) linear algebra against its definitions: the product entry by
entry, transposition reversing products, rank-nullity, kernel
vectors actually in the kernel and independent, rref idempotent
with cleared pivot columns, and solve returning None exactly when
the augmented rank exceeds the rank, which is Rouche-Capelli.
- Hamming(7,4) exhaustively: all sixteen nibbles by all seven single
errors corrected, all sixteen codewords pairwise at distance three
or more, and every one of the 336 double errors flagged and
miscorrected rather than passing silently.
- Fourteen named codes against their stated length, dimension and
distance, with G H' zero and H full rank in each.
- The perfect codes: Hamming, Golay(23) and the odd repetition codes
meet the sphere-packing bound with equality and nothing else does,
every coset leader is within the correction radius, and the
extended Golay code's covering radius is one past its correction
radius.
- Syndrome decoding corrects every error up to the radius and, on a
perfect code, provably lands on a different codeword one past it,
since there is nowhere else to land. The incremental search and the
standard array agree.
- MacWilliams's identity: the dual's weight enumerator computed by
enumerating the dual equals the Krawtchouk transform of the
primal's, exactly, for every code and every weight. Plus duality as
an involution, Golay(24) self-dual, and the dual of repetition
being the single parity check code.
- The Golay(24) weight distribution against the classical
1, 759, 2576, 759, 1.
- The four bounds, with repetition meeting Singleton and Plotkin at
once and the perfect codes meeting sphere-packing.
- LDPC regularity, and belief propagation against bit flipping on a
binary symmetric channel at two crossover probabilities chosen to
straddle bit flipping's threshold, so the soft decoder's advantage
shows as a difference in kind.
Two defects the tests found:
- hamming_74_encode used the wrong coverage masks for two of the
three parity bits: 0b0110011 and 0b0001111 rather than 0b1100110
and 0b1111000. The syndrome then named the wrong position and a
single error was "corrected" into a different nibble.
- ldpc_decode_bitflip flipped every bit tied for the most unsatisfied
checks. When the maximum is one, a large fraction of the block ties
for it, so the decoder flipped them all and oscillated: on a
length-240 code at five per cent crossover it ended with more
errors than it started, 684 against 665. It now keeps the iterate
with the fewest unsatisfied checks and returns that, which turns
the oscillation into a plateau: 217 of 495 on the same channel, and
4 of 181 at two per cent.
LinearCode::hamming and extended_hamming now carry their known distance
rather than searching for it, so the family reaches r = 8 instead of
stopping where enumerating 2^k codewords becomes impossible.
Co-Authored-By: Claude <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Part 4 session 11: src/codes/reed_solomon.rs and
src/codes/convolutional.rs. Completes roadmap item 7b.
reed_solomon.rs: GF(256) by logarithm tables, a prime field, a general
GF(2^m) with trace and minimal polynomials, systematic Reed-Solomon with
syndrome decoding through Berlekamp-Massey, Chien search and Forney's
formula, erasure decoding, the CCSDS, QR and disc parameter sets, binary
BCH codes, and an enumeration of the cyclic codes of a given length.
convolutional.rs: rate-1/n convolutional codes with hard and soft
Viterbi, free-distance search, puncturing and depuncturing, block,
random and quadratic-permutation interleavers, recursive systematic
encoders with trellis termination, max-log-MAP BCJR, turbo codes with
iterative decoding, Gaussian and binary symmetric channels, a bit error
rate sweep, and the capacity functions the whole subject is measured
against.
Twenty-five tests. The ones that carry weight:
- GF(256) exhaustively: the powers of the primitive element hit every
non-zero element exactly once, every element has the inverse it
should, and multiplication is commutative across all 65536 pairs.
GF(2^m) for m from two to eight likewise, with the trace shown to
land in GF(2), to be additive, and to split the field exactly in
half -- which is what being a surjective linear map onto GF(2)
means -- and every minimal polynomial shown to vanish at its own
root with a degree dividing m.
- Reed-Solomon corrects every error pattern up to (n-k)/2 on six
parameter sets, including the roadmap's RS(255, 223) against
sixteen random byte errors.
- A burst of 128 flipped bits confined to sixteen bytes is corrected
exactly, which is the property the code is deployed for and which
no bit-level code of that rate could match.
- Erasures cost half what errors do: n-k erasures are recovered,
which is the maximum distance separable property stated
operationally, and n-k+1 is refused.
- Past its capacity the decoder never returns a non-codeword: it
corrects to something valid or reports failure.
- BCH parameters against the classical table for nine (m, t) pairs,
with each generator verified to divide x^n - 1.
- The cyclic codes of length n are exactly the divisors of x^n - 1,
so the count must be two to the power of the number of cyclotomic
cosets. It is, for seven lengths.
- Free distances against the published generator tables: 5, 6, 7, 8
and 10 for constraint lengths three to seven.
- A terminated convolutional code is a block code of minimum distance
dfree, so Viterbi corrects any (dfree-1)/2 errors wherever they
fall. Checked for every count up to that on six codes.
- Soft decisions beat hard ones by the two decibels they are supposed
to be worth.
- Turbo iteration more than halves the error count against a single
pass, at a signal-to-noise chosen to sit in the waterfall.
- The Shannon limit for binary signalling against the values every
coding paper quotes: 0.187 dB at rate 1/2, -0.495 at 1/3, 1.059 at
2/3, and the -1.59 dB floor as the rate falls. The capacity behind
those is integrated numerically, so the agreement is a real check
on the integration.
One defect the tests found: the Reed-Solomon decoder had its position
convention backwards. Symbol zero is the leading coefficient of the
codeword polynomial, so position j carries x^(n-1-j) and its locator
value is alpha^(n-1-j); the Chien search instead read a root at
alpha^-i as naming position i, and Forney's formula carried a spurious
factor of alpha^i that belongs only to a generator whose roots start
elsewhere. The algebra verified against itself and corrected the wrong
symbols, so nothing short of an end-to-end decode would have caught it.
Co-Authored-By: Claude <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Part 4 session 12: src/codes/compression.rs and
src/codes/crypto_math.rs. Completes roadmap item 7c and section 7.
compression.rs: bit-level readers and writers, Huffman with canonical
codes, Shannon-Fano, arithmetic coding, LZ77, LZW, PackBits run lengths,
suffix and longest-common-prefix arrays, the longest repeated substring,
the Burrows-Wheeler transform and its inverse, move-to-front, delta
coding, byte entropy, and the normalized compression distance.
crypto_math.rs: RSA with Chinese-remainder decryption, Diffie-Hellman,
short Weierstrass curves over a prime field with the full group law and
the secp256k1 and P-256 constants, elliptic-curve Diffie-Hellman,
Shamir's secret sharing, the one-time pad, shift registers with the
Berlekamp-Massey attack on them, an avalanche measurement, the birthday
bound, and the classical cipher analyses. The module documents at length
that none of it is safe to deploy: every routine branches and indexes on
secrets, so the timing and the memory trace leak them.
Twenty tests. The ones that carry weight:
- Huffman is checked to be optimal, not merely valid: for alphabets up
to five, every length assignment satisfying Kraft is enumerated and
none beats it. Plus Kraft with equality, prefix-freeness by
construction, and Shannon's bound on both sides.
- Arithmetic coding is required to land within two bytes of the
message's own information content, and to beat Huffman where a whole
bit per symbol is too coarse a unit.
- The suffix array is checked against a naive sort of the suffixes and
the LCP array against naive comparison, on a three-letter alphabet
where ties are everywhere.
- The Burrows-Wheeler transform is required to preserve the histogram,
to halve the run count on repeated text, and to invert on periodic
strings, where the rotations tie and the sort order is ambiguous.
- The elliptic curve group law is checked to be a group -- identity,
inverses, closure, commutativity across every pair, and associativity
on a sample -- for four curves, and scalar multiplication against
repeated addition, with Lagrange and Hasse both holding.
- The published secp256k1 and P-256 generators are verified to lie on
their curves, to have the stated prime order, and to make scalar
multiplication a homomorphism at full 256-bit size.
- Shamir: every k-subset of the shares reconstructs and no (k-1)-subset
does, for every k and n up to five and seven.
- Berlekamp-Massey is required not just to report the right register
length but to predict the rest of the keystream, which is the actual
attack.
- Caesar is broken for all 26 shifts and Vigenere for four keys, with
Kasiski's suggestions required to include a multiple of the true
length.
Two defects the tests found:
- shamir_reconstruct dropped the minus sign in the Lagrange numerator,
computing the product of x_j where it needed the product of -x_j.
The two agree when the threshold is odd, so a three-of-n split
worked and a two-of-n split reconstructed the negation of the
secret.
- The shift register's step map is a bijection only when bit zero is
tapped: without it the outgoing bit does not reach the feedback, two
states share an image, and the register enters a cycle it never
started on. lfsr_period returned zero there with nothing saying why.
Both functions now document the requirement, and the test asserts
the zero rather than treating it as a short period.
Co-Authored-By: Claude <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Part 4 session 13, first half: src/stochastic/markov.rs, a new
stochastic/ module.
MarkovChain: construction with validation, estimation from counts or
from an observed sequence, stationary distributions by linear solve,
n-step transitions, simulation, reachability and state classification,
absorption through the fundamental matrix, hitting times and
probabilities, Kac return times, the mean first passage matrix, total
variation distance, mixing time, the spectral gap through the additive
reversibilisation, detailed balance, entropy rate, exact sampling by
coupling from the past, and the PageRank chain of a graph.
Mcmc: Metropolis-Hastings, an adaptive variant, Gibbs, Hamiltonian Monte
Carlo, the no-U-turn sampler, slice sampling, parallel tempering,
autocorrelation time, effective sample size, the Gelman-Rubin statistic,
and simulated annealing.
Eleven tests:
- The stationary distribution is checked to satisfy pi P = pi entry by
entry and to be the limit of the matrix powers, and separately that
a periodic chain still has one -- which is why it is solved rather
than iterated.
- Absorption against the gambler's ruin in closed form, for three
board sizes and four win probabilities, both the ruin probability
and the expected duration, with the hitting-time and
hitting-probability routines required to agree by their own separate
routes.
- Kac's formula, the mean-first-passage recurrence, and the entropy
rate bounded above by the stationary marginal's entropy.
- A random walk on a graph is shown reversible with stationary
distribution proportional to degree, and a biased directed cycle
shown stationary but not reversible.
- Coupling from the past is checked against the stationary
distribution over thirty thousand exact draws, and the PageRank
chain's stationary distribution against graph::spectral::pagerank.
- Metropolis-Hastings recovers a Gaussian's mean within four standard
errors, where the standard error is built from the effective sample
size rather than the run length.
- Hamiltonian and no-U-turn samplers are held to the mean, variance
and correlation of a target with correlation 0.9, and required to
achieve a higher effective sample size per draw than the
random-walk proposal does.
- Parallel tempering is required to cross between two modes separated
by twelve nats, at least ten times as often as a single cold chain
at the same proposal width, averaged over five starts.
- The diagnostics are checked to separate the two cases they exist
for: independent draws give an autocorrelation time near one and a
Gelman-Rubin statistic within 0.02 of one, while a correlated walk
gives a time above ten and chains started ten apart give a statistic
above two.
- Annealing escapes a local minimum in at least eighteen of twenty
runs, and a frozen schedule provably does not.
One defect the tests found: the first version of the no-U-turn sampler
ran a trajectory until it turned back and took the last point. That is
not reversible -- the trajectory length depends on the state in a way
the acceptance rule does not account for -- and it sampled a
distribution with variance 242 where the target's was one. Replaced with
Hoffman and Gelman's doubling scheme, where the trajectory grows
forwards or backwards at random, every sub-trajectory is checked for the
turn, and the next state is drawn uniformly from what the slice variable
admits.
Co-Authored-By: Claude <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Part 4 session 13, second half: src/stochastic/hmm.rs. Completes roadmap
item 8.
Hmm with scaled forward and backward recursions, Viterbi, posteriors and
posterior decoding, Baum-Welch over multiple sequences, and simulation.
GaussianHmm with the same interface for continuous emissions. A Kalman
filter sequence runner, the Rauch-Tung-Striebel smoother, the lag-one
smoothed cross-covariances, expectation-maximisation for the noise
covariances, and a bootstrap particle filter with systematic resampling.
Eight tests:
- The forward recursion is checked against literal enumeration of every
state path, for sequence lengths one to nine, and Viterbi against the
maximum over the same set -- so the recursions are verified against
the definitions they are shortcuts for, not against each other. The
posteriors are checked the same way, path by path.
- Scaling is checked to do its job: a five-thousand-symbol sequence
still gives a finite log-likelihood, which an unscaled recursion
would not.
- Viterbi recovers the exact state path on a chain whose emissions name
the state, and decodes the occasionally-dishonest casino correctly
more than seven times in ten.
- Baum-Welch is checked round by round, not end to end, for the
monotonicity that is its only guarantee; and required to find the
loaded die's loaded face from a random start.
- The Gaussian model recovers well-separated states more than
ninety-five per cent of the time and learns means of -2 and 3 from a
start at -0.5 and 0.5.
- The smoother's variance is required to be no larger than the
filter's at every step and every component, its total error smaller,
and its final estimate identical to the filter's -- since there is no
future for the last step to borrow from.
- EM recovers a measurement noise of 0.8 from a start of 0.01, keeps
both covariances symmetric and positive semidefinite, and leaves the
truth alone when started there.
- The particle filter is required to track the Kalman filter to within
0.1 over sixty steps on a linear Gaussian model, where Kalman is
exactly optimal, with the effective particle count restored by each
resampling.
One defect the tests found: the expectation-maximisation step estimated
the noise covariances from the residuals of the smoothed states alone,
with no covariance terms. The smoothed states are shrunk towards each
other, so their residuals understate the process noise, and the
measurement residuals ignore the smoother's own uncertainty; from a
deliberately wrong start it converged to a measurement noise of 3.7
where the truth was 0.8. Replaced with the closed-form maximisation in
the three second-moment sums, which needs the lag-one smoothed
cross-covariances -- so those are now computed and exposed rather than
implicitly assumed to be zero.
Co-Authored-By: Claude <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Part 4 session 14: src/stochastic/sde.rs and
src/stochastic/point_process.rs. Completes roadmap item 9.
sde.rs: Brownian motion, bridges and their higher-dimensional forms;
geometric Brownian motion and Ornstein-Uhlenbeck stepped exactly;
Euler-Maruyama in one dimension and in n; Milstein; the Stratonovich
Heun scheme; an order-1.5 scheme for additive noise; measured
convergence orders; Cox-Ingersoll-Ross and Heston by full truncation;
Merton jump diffusion; stable sampling by Chambers-Mallows-Stuck;
fractional Brownian motion by Davies-Harte with a Cholesky fallback;
rescaled-range and detrended-fluctuation Hurst estimators; first passage
by simulation and in closed form, both density and distribution;
Feynman-Kac against Black-Scholes; Ito's isometry; BAOAB Langevin
dynamics; the Chang-Cooper Fokker-Planck solver and the stationary
density it should reach; Kramers' escape rate; stochastic resonance.
point_process.rs: Poisson processes in time, space and with a varying
rate; compound Poisson; Hawkes with intensity, likelihood, simulation
and fitting; renewal and Cox processes; Matern and Thomas cluster
processes; Ripley's K, Besag's L and the pair correlation with edge
correction; the Clark-Evans index; quadrat and Kolmogorov-Smirnov tests;
Galton-Watson branching with its extinction probability; Yule and
birth-death processes.
Twenty-two tests. The ones that carry weight:
- Euler-Maruyama is measured at strong order one half and Milstein at
one, path by path against the exact solution driven by the same
noise. That is the only statement distinguishing a correct Milstein
from an Euler step with a small extra term.
- Heun and Euler-Maruyama are shown to disagree about the same
equation in exactly the way Ito and Stratonovich do: for dX = X dW
the Ito solution is a martingale with mean one and the Stratonovich
one has mean exp(1/2). Both are measured.
- The Fokker-Planck solver is required to conserve probability to
within a part in a billion and to relax pointwise onto the closed-
form stationary density, which is separately checked against the
Gaussian it should be.
- Langevin dynamics is held to equipartition in the velocity and the
Boltzmann distribution in the position, both to within six per cent.
- First passage times are compared to the exact distribution over the
whole sample rather than bin by bin, and the two closed forms --
density and distribution -- are checked against each other by finite
difference. A downward drift is required to reach the barrier with
probability exp(2 mu b / sigma squared) and no more.
- Poisson counts are checked against the mass function term by term
for thirty values of k, not merely on their mean and variance.
- Clustered patterns are required to be detected as clustered by all
four spatial diagnostics, and a random one by none of them.
- Hawkes: the intensity against its definition, the stationary rate
against mu over one minus the branching ratio, the likelihood shown
to fall when any parameter is displaced, and the fit recovering the
branching ratio to fifteen per cent.
- Galton-Watson extinction against the generating function's smallest
fixed point, including the critical case where a population that
replaces itself on average still dies out with probability one.
Four defects the tests found:
- cir_process clipped its state at zero, which is reflection rather
than full truncation. Every reflection injects probability mass the
exact process does not have, and the bias grew as the step was
refined -- 0.03 at dt = 1e-3 and 0.10 at 5e-5 against a true mean of
0.02. Full truncation keeps the internal state signed and truncates
only inside the coefficients and on output. Heston's variance leg
had the same fault.
- levy_stable_sample skewed the opposite way from the convention it
documented: a positive beta stretched the lower tail. Now pinned by
a test in both directions.
- hawkes_process resummed the whole history at every proposal, so a
run with n events cost n squared and an eight-thousand-unit horizon
did not finish. The excitation is now carried forward, which is one
exponential per proposal.
- first_passage_time_sim checked for crossings only at grid points and
so missed the excursions that cross and return within one step,
biasing the times upward. Conditional on the two endpoints the
probability that the bridge between them touched the barrier has a
closed form, so those crossings are now counted.
hurst_dfa now skips windows below sixteen points and documents that it
wants increments rather than an integrated series: removing a straight
line from eight points takes real fluctuation with it and biases the
exponent up by enough to make white noise look persistent.
Co-Authored-By: Claude <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
The Miri job has been hitting its 45-minute cap and being cancelled, so the check has not reported since the workflow was repaired. The build was never the problem -- the cached sysroot and crate finish in under twenty seconds. The scope was. A libtest filter is a substring match rather than a path prefix, so `core::` selected `graph::core::` and `verification::core::` alongside the intended module. Those two arrived after the job was written and quietly doubled it from 27 tests to 54, adding brute-force combinatorial tests that Miri, interpreting at roughly a hundredth of native speed, cannot finish: girth_matches_brute_force alone ran for nine and a half minutes and hamiltonian_path_matches_brute_force for five before the cancel landed. Skipping `::core::` anchors the filter to the top level. A top-level path starts with `core::` and so has nothing preceding it to match, while every nested one does, which keeps the anchor working for any module named `core` added later. That restores the intended 27 tests, which the cancelled run showed completing in nine and a half minutes -- comfortably inside the existing budget, so the timeout stays as it is. Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
queueing.rs covers the birth-death queues and their closed forms -- M/M/1, M/M/c, the finite-capacity and infinite-server variants, Erlang B and C, Pollaczek-Khinchine and Kingman -- along with Jackson networks, an event-driven simulator, and continuous-time Markov chains with transient solutions by uniformization. The tests lean on two theorems that hold across the whole module rather than checking each formula in isolation. Little's law is a statement about areas under a sample path, so L = lambda W has to hold for every model at every admissible parameter, and the simulator measures its time averages by integrating a merged event list rather than by invoking the law, which makes the comparison evidence instead of arithmetic. Summing n p_n against the reported mean ties each distribution to the means derived separately from Erlang C. The models also have to nest: M/M/c at c = 1 is M/M/1, M/M/1/K converges to M/M/1 as the buffer grows, and P-K with an exponential second moment reproduces M/M/1 while P-K with a deterministic one halves it. The simulator is checked against both, so M/D/1 separates the general formula from the exponential special case. Priorities are checked against Cobham's formula and against Kleinrock's conservation law, which fixes sum rho_k Wq_k however the queue is ordered. Three defects surfaced in writing those tests. Forming p_n as a^n / n! overflows both halves independently and returns NaN for a tail the recursion handles without trouble; every such expression is now a running product. Uniformization at a long horizon failed in both directions at once -- past Lt = 745 the leading Poisson weight underflows to zero, and substituting the smallest subnormal made the recurrence climb by roughly e^Lt on its way to the mode and overflow to infinity -- so the weights are now carried as logarithms and materialised only over the window where they are representable. The transient solver's state-space bound followed the free drift even for a stable queue, which is positive-recurrent and stays near its stationary law however long it runs. timeseries.rs covers correlation structure, the two opposed stationarity tests, ARMA/ARIMA/SARIMA, exponential smoothing, GARCH, Granger causality, vector autoregressions, cointegration, seasonal decomposition, changepoint detection, the entropy measures, IAAFT surrogates, and the local level model. The p-values for the Dickey-Fuller and KPSS statistics come from tabulated quantiles of their own non-standard null distributions rather than from a t or a chi-squared, which would be wrong rather than approximate; the Engle-Granger table is kept separate from the plain Dickey-Fuller one because the residual being tested is fitted rather than observed. The two tests take opposite nulls, so the pair is checked on the same data in both directions. Elsewhere the tests assert identities: the spectral density integrates to the variance the impulse-response weights give, the partial autocorrelation cuts off exactly past the autoregressive order, the forecast band converges to the process standard deviation for a stationary model and grows without bound for an integrated one, EWMA is GARCH at zero intercept and unit persistence step for step, and every fit is checked by whether the residuals it leaves are white. Holt's method was seeded half a step ahead of where its recursion expects the state, leaving a transient on data it should reproduce exactly -- an exact straight line. Seeding one step before the data, with the seasonal factors taken against the trend line rather than a flat mean, makes both Holt and Holt-Winters exact from the first prediction on a noiseless trend-plus-season. The Hannan-Rissanen pilot order was clamped to a range that inverts on a short series, panicking where the function documents an error. The shared regression kernel solves the normal equations by Cholesky rather than taking a QR of the design. The crate's Householder QR accumulates an explicit n-by-n orthogonal factor, which is O(n^2 k) and reached a billion operations for a long pilot autoregression; the normal equations are O(n k^2) and every regression here has far more rows than columns. That cut the module's tests from 74 seconds to 1. The trade -- squaring the condition number -- is documented, and the substitution is checked against the QR it replaces. The property suite adds the cross-module theorems none of these modules can check alone: Little's law across every closed form at once, the two independent routes from a birth-death chain to its stationary distribution, the identity between a continuous-time chain and its embedded discrete one reweighted by holding times, and the averaged periodogram of a simulated ARMA against the density the model computes from its coefficients -- an FFT that knows nothing about ARMA models against a formula that knows nothing about the FFT. Random stationary autoregressions are drawn through the Barndorff-Nielsen-Schou map from partial autocorrelations, which lands on the stationary region by construction. 3,390 library tests and 144 property tests pass; clippy is clean under --all-targets -D warnings. Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
rmt.rs covers the classical ensembles -- GOE, GUE, Ginibre, Wishart -- the limiting spectral laws they converge to, the local statistics that distinguish a correlated spectrum from uncorrelated levels, and Marchenko-Pastur denoising of a sample correlation matrix. The tests check universality rather than arithmetic. A GOE spectrum is matched against the semicircle by a Kolmogorov-Smirnov distance and its second moment against the value the scaling fixes exactly; a Wishart spectrum has to sit inside the Marchenko-Pastur band, which at an aspect ratio of a quarter spreads purely noisy eigenvalues over a factor of four even though every population eigenvalue is one. The spacing statistics carry the real content: the ratio of adjacent gaps needs no unfolding, since the local density cancels between numerator and denominator, and it separates 0.5307 for the orthogonal class from 2 ln 2 - 1 for independent points -- and separates the orthogonal class from the unitary one at 0.5996, which is the whole content of the symmetry classification. Spectral rigidity separates the same two cases far more sharply, L/15 against something growing like a logarithm. Every surmise is checked to integrate to one with unit mean, which constrains the prefactors rather than merely the shape. Two of these needed care in the test rather than the code. The spectral densities vanish like a square root at both edges, and Marchenko-Pastur at unit aspect ratio picks up an inverse-square-root singularity where its lower edge reaches zero; midpoint quadrature converges at only h^(1/2) there, so a uniform grid measures the quadrature and not the density. A sine substitution cancels both. Separately, the Jacobi eigen-solver is cubic per sweep, so the module's tests were sized down from matrices whose cost dominated the suite; one test also carried a pooled spectrum left over from an earlier approach that was built, sorted, and never read. extreme.rs covers the generalised extreme value and generalised Pareto families with maximum-likelihood fits, return levels, the Hill estimator, the Ferro-Segers extremal index, rank correlations, and five copula families with sampling, Kendall-tau inversion, tail dependence and the Pickands dependence function. The two routes into a tail are checked against each other: fitting a GEV to block maxima and a generalised Pareto to threshold exceedances of the same data recovers the same shape parameter, which is Pickands-Balkema-de Haan and the reason the threshold route is worth preferring. Return level and return period are checked as exact inverses; the mean excess is checked to be linear in the threshold with slope xi/(1-xi), which is what makes it a threshold diagnostic; the extremal index is checked against a moving maximum over a window of m, whose index is 1/m. The copula tests assert uniform margins for every sampler, recover each family's parameter by inverting Kendall's tau, and check that the Pickands function stays between max(t, 1-t) and 1. Tail dependence needed the tests restated. A coefficient that is zero only asymptotically is not zero at a finite quantile: Gumbel's lower coefficient at q = 0.01 with theta = 2 is q^(2^(1/2) - 1), about 0.15, and the Gaussian's decays only logarithmically. So the tests assert the exact finite-q values where a closed form exists, and otherwise assert the asymptotic statement directly -- that the coefficient falls as the quantile tightens. The Gaussian and t comparison is put the same way: at a loose quantile the two are nearly indistinguishable, and the gap widens monotonically as the quantile tightens, with the t settling on its exact limit while the Gaussian decays away. That is the failure mode a correlation-based risk model cannot see, and stating it as a widening gap rather than a fixed threshold is what makes the test mean it. Kendall's tau compared every pair, which is the definition but costs O(n^2) -- minutes on the sample sizes a copula fit wants, and 41 seconds of the module's own tests. Sorting by x and counting inversions in the resulting y sequence by merge sort gives the same discordance count in O(n log n), with the tie corrections handled by Knight's formula. That took the tests to under a second, and the replacement is pinned against the pair-counting definition over randomised samples carrying ties in one coordinate, the other, and both. The property suite adds the identities that cross modules: the even moments of Wigner's semicircle are the Catalan numbers, computed here by a numerical integral of a density and there by an exact integer recurrence in discrete::combinatorics, with nothing shared between the two routes. Alongside them, a noise covariance staying inside its predicted band across random aspect ratios, denoising preserving the trace while never widening a spectrum, and copula parameters surviving a round trip through a rank statistic that ignores the margins entirely. 3,433 library tests and 153 property tests pass; clippy is clean under --all-targets -D warnings. Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Adds the simplex method, a primal-dual interior point method, LP duality, sensitivity analysis, a small modelling language, and the classical models that reduce to a linear program: diet, production planning, transportation, zero-sum games, the Chebyshev centre, and L1 and minimax regression. Placed at optimization/lp.rs rather than the roadmap's opt/lp.rs. The crate already has an optimization module covering exactly this subject area, and two top-level modules named opt and optimization would be a lasting wart for the sake of matching a path. Every signature the roadmap names is present; only the directory differs. Duality is the organising idea and the module commits to one convention, stated in the module documentation and adhered to throughout: duals[i] is the derivative of the reported objective with respect to b[i]. That is the definition that makes shadow prices mean what people expect and makes sensitivity ranges checkable, and it is what the tests check -- perturbing a right-hand side within its reported range moves the objective by exactly the shadow price times the perturbation, over hundreds of random programs and both signs of perturbation. The tests lean on theorems that are exact rather than asymptotic, so they demand equality. Strong duality is an equation: the objective equals the right-hand side dotted with the shadow prices, and nothing in the solver imposes it -- the duals are read off the optimal basis and the objective off the primal solution. Complementary slackness holds at every optimal basis, in both directions: a variable in use has zero reduced cost, and a slack row has zero shadow price. The two solvers walk the feasible region in completely different ways -- one along the boundary vertex to vertex, the other through the middle, never reaching a vertex -- and share only the standardisation step, so their agreement on every instance is the strongest available check on either. The dual of the dual returns the primal value; weak duality brackets every feasible pair; and the minimax theorem falls out as a corollary, since the two players' programs are duals of one another. Bland's rule is used throughout rather than a faster pivoting rule. Degeneracy is not hypothetical -- Beale's example returns to its starting basis after six pivots under Dantzig's rule and cycles forever -- and Bland's rule cannot cycle because the basis sequence it visits is lexicographically monotone. Beale's example is in the tests, alongside a problem where three constraints meet at a single vertex so that the ratio test ties at every pivot. Two things the tests caught that were wrong in the tests rather than the code. The Chebyshev centre is not unique: in a box four wide and six tall the largest inscribed circle has radius two and slides freely up and down, so only the coordinates the touching faces pin down are determined. The radius is unique and the fit is checked instead by the properties that do hold -- the ball fits inside every face, and touches at least one, so nothing larger fits. Separately, a Matrix cannot be constructed with a zero dimension, which made two emptiness guards unreachable; an untestable branch is worse than none, so they were removed rather than left as decoration. Total unimodularity gets a test of its own: the transportation problem's constraint matrix is totally unimodular, so integral supplies and demands give an integral optimum straight from the simplex method, with no branch and bound anywhere. The regression fits are checked against each other and against ordinary least squares under all three norms -- each must win under the norm it minimises -- and the minimax fit is checked to be pinned by at least three residuals of alternating sign, which is what distinguishes a minimax fit from merely a fit with a large residual. 3,463 library tests and 162 property tests pass; clippy is clean under --all-targets -D warnings. Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
integer.rs adds branch and bound over the linear relaxation, Chvatal-Gomory rounding cuts, the knapsack family, subset sum and partition, bin packing, set cover, facility location, cutting stock by column generation, the dynamic programming classics, and exact cover by Algorithm X with sudoku, n-queens and AC-3 on top of it. network.rs adds transshipment, the critical path method and PERT, linear programming formulations of shortest path and maximum flow, Clarke-Wright vehicle routing, and the sequencing rules. Every method is checked against an independent exact answer rather than against itself. The knapsack table and the branch-and-bound search tree share no code and must agree on every instance, and both are checked against enumeration. Each scheduling rule is checked against every permutation of the jobs, on the objective it provably optimises and on nothing else -- and shortest-processing-time is checked to be genuinely worse than earliest-due-date on maximum lateness, because a rule that looked good on every objective would mean the test was not measuring what it claimed. The greedy methods are checked against their proven ratios, against exact answers: first-fit-decreasing within 11/9 OPT + 6/9, set cover within the harmonic number, longest-processing-time within 4/3 - 1/(3m). A bound nobody tests against an optimum is not a guarantee. The two linear programming formulations exist to cross-check the graph module. Shortest path and maximum flow both have totally unimodular constraint matrices, so a general simplex solver answers them exactly, and it agrees with Dijkstra and with the augmenting-path search over hundreds of random graphs. The two routes share nothing but the graph itself. The first cut generator was wrong. It rounded a variable's bound toward whichever side the objective was not pushing, which is a branch rather than a cut: it removed integer solutions, and the test caught it discarding an optimum of 34 while reporting a bound of 33. Replaced with the Chvatal-Gomory rounding cut, whose validity rests on an argument that can be stated -- scale a row, round the coefficients down, and the left-hand side becomes an integer, so it is bounded by the floor of the right -- and which therefore holds for every non-negative integer point. The test now enumerates every integer point in the box and checks that none is removed. Separately, Miri went red on the previous commit, in core::interval::tests::test_sqrt_exp and not in anything this change touches. Reproduced locally and deterministic, so it is a newer nightly Miri rather than a flake, and it will recur every run until addressed. Instrumenting it showed two calls to exp(1.0) returning values several ulps apart, while Interval::exp widens its bounds by a fixed two. That widening silently assumes the host evaluates the elementary functions to within one ulp -- true of every real platform's libm, not required by the language, and deliberately untrue under Miri. So the enclosures this module returns are rigorous conditional on that assumption, which is now documented on the widening helpers where it belongs. The test itself joins the five siblings in the same file that already carry cfg_attr(miri, ignore) for the same root cause. It still runs on every commit in the ordinary test job; what the attribute records is that one interpreter deliberately violates the test's precondition. Widening every interval for all users to accommodate an interpreter that randomises would make the library worse, not more rigorous. 3,504 library tests and 170 property tests pass; clippy is clean under --all-targets -D warnings. Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Adds spin.rs (spin operators at any spin, coherent states, XXZ chains with a matrix-free Hamiltonian, Lanczos, Krylov time evolution, the transverse-field Ising chain and its Jordan-Wigner solution, Larmor and Rabi dynamics, Bloch equations, free induction decay) and solid_state.rs (tight binding, SSH, graphene, Kronig-Penney, densities of states, occupations, Debye and Einstein and Sommerfeld heat capacities, phonons, Landau levels, Hofstadter, transport, semiconductors, BCS, Josephson, Anderson localisation). Three defects the tests caught. Lanczos returned correct eigenvalues and wrong eigenvectors. eigen_symmetric_tridiagonal hands back each eigenvector as a *row*, and I indexed it as a column. The eigenvalues are unaffected -- they come from the tridiagonal projection, which the indexing does not touch -- so the error was invisible in every energy and showed only in the residual. It had propagated into the Krylov time step as well, where it made the evolution non-unitary. Larmor precession turned the wrong way. dM/dt = gamma M x B puts the angular velocity at -gamma B, so a positive gyromagnetic ratio precesses clockwise seen from +z; the closed form turned anticlockwise and disagreed with the Bloch integrator in the same module. effective_mass_from_band called every band flat. The guard compared the curvature against an absolute 1e-30, and a real band in SI units curves by about 1e-38, so the function refused the case it exists for. The threshold is relative to the band's own scale now. Four test premises of mine were wrong. The magnon band tops out at 4 j s, not 2 j s. The critical Ising chain is *not* the most entangled one at a given size -- deep in the ordered phase the ground state is the symmetry-broken cat and carries a full bit across every cut, more than the critical chain -- so the test now checks the *scaling*, which is what distinguishes criticality: the half-chain entropy climbs at c/6 = 1/12 of a bit per doubling while both phases saturate. The free induction decay test asked for a tail below 1e-3 from a record that was truncated rather than decayed. And I read a 20 per cent gap in the Anderson localisation ratio at the band centre as the Kappus-Wegner anomaly; at 150,000 sites instead of 20,000 both ratios are 4.0 and the gap was sampling noise, so the claim is gone and the test says why the chain has to be long. Adds tests/properties/quantum_matter_props.rs: the angular momentum algebra at every representation up to spin eight, coherent states pointing where they were asked to, Lanczos eigenpairs certified by their own residual and by the variational principle, Krylov evolution unitary at every step size, the Ising chain against its free-fermion energy at every field, SSH edge counts against the bulk winding number on random couplings, and the occupation functions' exact symmetries -- including the point at which the boson-fermion gap falls below what a double can represent, where the test stops demanding a strict inequality. Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Roadmap section 16. `statistical_mechanics.rs` becomes a directory so the new material can sit beside the thermodynamics already there; the roadmap calls the home `statmech/`, but every earlier session has kept new modules under the existing names and this follows that. ising.rs -- Ising2D with Metropolis, heat-bath and Wolff cluster updates, correlation functions and lengths, autocorrelation times, the Onsager magnetisation and energy, the exact one-dimensional chain by transfer matrix, brute-force enumeration for small systems, Potts and XY models with plaquette vorticity, Wang-Landau sampling with canonical reconstruction, parallel tempering, Binder crossings and a fluctuation-dissipation check. lattice_models.rs -- site and bond percolation with a disjoint-set spanning test, cluster size distributions, exact self-avoiding walk enumeration and Rosenbluth sampling, the connective constant, lattice random walks and Polya return probabilities, Flory exponents, Kasteleyn dimer counts, KPZ ballistic deposition with interface widths and growth exponents, the Abelian sandpile, and the Clauset power-law fit. Defects found and fixed while writing the tests: - `sample` drove the chain with a *state-dependent* stopping rule: a "sweep" of Wolff steps ran until the flipped total reached the lattice size, so a measurement was always taken just after a large cluster. That biases the sample toward ordered configurations, and it showed: -3.90 per site against the exact -3.29. Each measurement now follows a fixed number of updates, and `autocorrelation_time` returns the work per update alongside tau so the two updates can still be compared honestly. - `correlation_length_estimate` fitted a length from a single snapshot, which is a fit to noise -- the spread of the correlation function at large separation is comparable to its mean and does not shrink with the lattice. Added `sample_correlations` to average over a run, and the estimate now takes the ensemble average and its background. - `connective_constant_estimate` read consecutive ratios, but the walk counts alternate with parity, so consecutive-ratio Richardson made the answer worse rather than better. Now averaged over a parity pair before extrapolating: 2.63928 against the true 2.638158. Defects in the tests themselves, recorded rather than quietly patched: - The negative control on the connective constant demanded the raw ratio be 0.05 away from the truth; it is 0.0494. It now compares the two errors directly, which is what it was meant to show. - `saw < 4^n` is not strict at n = 1: the first step cannot revisit, so the counts coincide there. - The Potts relation was written `2 * ising_tc_exact() / 2`, which is just `ising_tc_exact()`. The relation is `ising_tc_exact() / 2`. - Kasteleyn's formula is exact only to rounding, so the 1-by-2 count comes back as 0.9999999999999999 and a strict `>= 1.0` fails. - The one-dimensional magnetisation saturates in `beta * h`, not in `h`; at beta = 0.05 a field of 50 reaches 0.988, not 1. - The weighted and unweighted polymer means happened to coincide at 388/6 for the weights I first picked, so the control could not have failed. tests/properties/statmech_props.rs -- 29 property tests. The strong ones are structural rather than statistical: percolation is checked for monotonicity under a common random number, which couples two lattices exactly and makes spanning monotone lattice by lattice rather than on average; the total vorticity of an XY torus is checked to vanish before and after thermalisation, since windings can only be created in pairs; the dimer count is checked against the Fibonacci recurrence for a two-row strip; the Monte Carlo sampler is checked against an exact enumeration of the same sixteen-spin lattice under both updates; and the density of states is checked against that enumeration at five temperatures, with the fluctuation heat capacity checked against a finite difference of the mean energy -- two independent routes that agree only if both are right. 3717 lib tests and 267 property tests pass in debug; clippy is clean under --all-targets -D warnings. Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Roadmap section 17, first half. `statistical_mechanics/md.rs` carries the pair potentials, a cell-list force evaluation with the minimum-image convention, velocity Verlet, Berendsen/Nose-Hoover/Langevin thermostats and a Berendsen barostat, the radial distribution function, structure factor, Lindemann ratio, mean squared displacement, velocity autocorrelation and vibrational spectrum, plus Ewald summation, the second virial coefficient, kinetic-theory lengths, a Green-Kubo transport integral, WHAM umbrella sampling, steered pulling and the Jarzynski estimator. Everything is in reduced Lennard-Jones units with k_B = 1, stated in the module header and in `lj_reduced_units_note`. The equations of motion are consistent under any consistent choice of units and silently wrong under an inconsistent one, so there is no SI path through this module at all. Two deviations from the roadmap's signatures, both recorded in the source. `MdSystem` carries its own cell list rather than the general-purpose `SpatialHash`, which owns a copy of every position and knows nothing about periodic images. And it carries `unwrapped` positions alongside the wrapped ones, because a mean squared displacement taken from wrapped coordinates saturates at the box size and reports no diffusion however freely the particles are moving. Defects found while writing the tests: - The cell list allocated one cell per cutoff per box edge with no bound. A dilute system -- a four-hundred-sigma box with a cutoff of a tenth -- asks for sixty-four billion cells for a few hundred particles, and the process aborted on a 1.5 TB allocation. The counts are now halved until the grid is comparable to the particle count; a cell larger than the cutoff is still correct, only less selective. - `green_kubo_viscosity_lite` integrated the stress autocorrelation to half the record. Past a few correlation times that estimate is noise, and integrating thousands of such lags accumulates a random walk as large as the signal: on an Ornstein-Uhlenbeck series whose integral is 1.01 by construction, it returned 0.14. It now stops at the first non-positive lag, which costs a few per cent of the tail and returns 0.99. - The integrator caches the force between steps -- that is what makes velocity Verlet one evaluation per step rather than two -- but `pos` is public, so writing to it left the cache stale and the next step integrated the previous configuration's forces. The symptom was quiet: energy that almost conserved, and a trajectory that was no longer reversible. Added `refresh_forces` and documented the requirement at the field and at the method. The reversibility property test is what found it, which is the kind of thing it exists for. Defects in the tests themselves, recorded rather than quietly patched: - The Madelung constant is defined by the energy of *one* ion in the field of all the others, while a lattice energy counts each pair once, so the total per ion is half of it. The implementation was right and the expectation was off by that factor of two. - A Lennard-Jones pair truncated at a tenth of sigma is not a free gas: it has a 10^12 core hidden just inside the cutoff, and two diffusing particles eventually find it and are ejected at enormous speed. The ideal-gas fixtures now use a genuinely zero potential. - The Lennard-Jones force changes sign at 2^(1/6) sigma, not at sigma -- the *energy* crosses zero at sigma and the pair is still repelling there. - Particle 0 sits at the lower x in the two-body fixture, so a repulsive radial force points along -x; the sign expectation had the geometry backwards. - The Debye structure factor below 2 pi / L measures the sample's extent rather than its structure and rises toward N, so the search for a Bragg peak found that instead. Documented and the search window moved above it. - Continuity at the cutoff was tested against a fixed tolerance, which only tested the step size I happened to pick. It now measures the gap at two step sizes and checks that halving one halves the other. - Three FCC cells at liquid density give a box only two cutoffs across, so the cell-list comparison was silently running the fallback against itself. Raised to five. - `energy_drift` is a fitted slope, so a partial oscillation genuinely does register as a small trend. The property test asked for more than a linear fit can give; it now checks the closed form on a pure trend and compares a wobble against a trend twenty times its amplitude. The tests lean on closed forms wherever one exists: the harmonic pair's exact solution, the Einstein relation D = T / (m gamma) with its exp(-gamma t) velocity autocorrelation, the hard-sphere B2 = 2 pi d^3 / 3, the Madelung constant, the Fibonacci-like recurrences of the reduced-unit fixtures, and a WHAM inversion of histograms built exactly from a chosen profile so the recovery has no statistical error to hide behind. Where no closed form exists the check is a comparison against an independent route -- the cell list against the direct O(N^2) loop, the fluctuation heat capacity against a finite difference, velocity Verlet against explicit Euler on the same trajectory. tests/properties/md_props.rs adds 19 property tests. The strongest is exact reversibility: run forward, negate the velocities, run the same number of steps, and every particle returns to within 1e-7 of where it started. 3757 lib tests and 286 property tests pass in debug; clippy is clean under --all-targets -D warnings, and the module checks on nightly-2025-11-21. Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Roadmap section 17, second half. `statistical_mechanics/kinetics.rs` carries reaction networks and their stoichiometry, an adaptive implicit integrator for the rate equations, Gillespie's exact stochastic algorithm and explicit tau-leaping, enzyme saturation with fits and the three inhibition mechanisms, equilibrium composition by Newton on the logarithms, the Brusselator, Oregonator and chemical Lotka-Volterra, Eyring and transition-state rate theory with the Kramers correction, nucleation and Avrami transformation, and the acid-base and electrochemical relations. The elementary single-formula relations already in `chemistry` -- the Arrhenius rate, the equilibrium constant from a free energy, the Nernst potential, pH from a proton concentration -- are used rather than duplicated. What is here is the part that needs a solver. Defects found while writing the tests: - `rate_equations` was built on BDF2, which assumes a uniform step. An adaptive controller varies the step every step, so the history was at the wrong spacing, that inconsistency dominated the error estimate, and the controller shrank the step in response until the integration stalled and gave up. Replaced with backward Euler plus Richardson extrapolation: a one-step method has no history to get wrong, and the extrapolated value is second order anyway. - The step-doubling error estimate can be fooled outright on an oscillatory system. An L-stable method damps hard at a step much longer than the period, so the coarse and fine solutions both collapse toward the fixed point, agree closely with each other, and report a small error -- and the controller then grows the step further. A run can step clean over whole oscillations while its error estimate reports success. The step is now also bounded by the solution's own timescale, |c| / |dc/dt|, which looks at the dynamics rather than at the difference between two equally wrong answers. - `mass_action_rates` looked up only the reactant concentrations, so a composition too short to cover the products was silently accepted and the wrong system integrated. It now checks every species the network mentions. - `steady_state_approx_check` skipped a fraction of the *steps* before measuring, and the adaptive integrator front-loads its steps into the induction period -- exactly the region where the approximation is not claimed to hold. It now skips an initial *time*, fifty complex-filling times, and refuses a run that ends before then. Defects in the tests themselves, recorded rather than quietly patched: - Two tests sampled a trace by step index rather than by time and drew the wrong conclusion from the transient: the stiff network's quasi-equilibrium read 0.33 instead of 0.5, and the "has this settled" helper reported a swing where there was none. Every such question here has to be asked of a time window. - The same helper then demanded a dense tail, which fails on precisely the runs that are most obviously converged: a settled system produces one enormous final step, and that single sample is itself the evidence. - The Lineweaver-Burk bias is a statement about *additive* noise. My fixture applied noise proportional to the rate, which survives the transform unchanged -- the two fits erred by 0.2202 and 0.2199, and the test proved nothing until the noise model was corrected. - Henderson-Hasselbalch fails for a *dilute* buffer, not a lopsided one. At 0.101 M acid with 0.1 M base the shortcut and the full balance agree to five decimals; at a micromolar they differ by more than a unit. - `jmak_avrami` reaches exactly 1.0 in double precision once (k t)^n passes about 37, so a strict "less than one" was testing the float format rather than the curve. - 2.302_585 is ln(10) to seven digits, a relative error of 4e-8, which exceeded the 1e-9 tolerances two property tests used. The tests lean on closed forms where they exist -- first-order decay, the quadratic for a weak acid, the logistic ignition time, the Fibonacci-free exact Poisson moments -- and on independent routes where they do not: the Gillespie mean against the rate equations for a linear network, where the two agree exactly in the mean; tau-leaping against Gillespie as the leap shortens; the inhibition mechanisms refitted rather than inspected; and the implicit integrator's step count against the seven million an explicit method would need on the same stiff system. tests/properties/kinetics_props.rs adds 13 property tests, several over randomly generated mass-balanced networks so the conservation law being checked is one the integrator has no way to know about. 3793 lib tests and 299 property tests pass in debug; clippy is clean under --all-targets -D warnings, and the module checks on nightly-2025-11-21. Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Roadmap section 18, first module. `biophysics.rs` becomes a directory so the population-scale models can sit beside the membrane and transport relations already there; the roadmap calls the home `bio/`, but every earlier session has kept new modules under the existing names and this follows that. epidemiology.rs -- SIR, SIS, SIRS, SEIR, SEIRS and MSIR on an adaptive Runge-Kutta with Richardson extrapolation; R0, the herd immunity threshold and the final-size equation by bisection; branching-process extinction; vaccination, demography, two competing strains and an age-structured model with its next-generation R0; the network epidemic threshold from the adjacency spectral radius; an exact stochastic SIR and an SIR on a contact graph; and the Cori and Wallinga-Teunis reproduction-number estimators with serial-interval and SEIR fitting. Defects found while writing the tests: - `final_size_equation` bisected with the sign inverted. The function 1 - z - exp(-R0 z) vanishes at zero and is *positive* just above it -- its slope there is R0 - 1 -- and negative at one, which is the opposite of the usual arrangement. Every final size came back as 1e-12. - `sir_with_vaccination` added the vaccinated fraction to the removed class a second time. `sir` already places 1 - s0 - i0 there, which for s0 = 1 - coverage - i0 is exactly the vaccinated, so the population summed to 1 + coverage. - `seir_fit_to_incidence` ran one Nelder-Mead pass, which contracts onto a direction and stops exploring the others -- and this objective has exactly the valley that punishes: the growth rate constrains beta and sigma only in combination. Added a restart from the best point. Defects in the tests themselves, recorded rather than quietly patched: - I asserted competitive exclusion for the two-strain model, which is false for it. Exclusion is a statement about a system that replenishes its susceptibles; a one-off epidemic is a finite race, and a strain with a thousandfold head start out-infects a rival with nearly twice its reproduction number before the susceptibles run out. The documentation said the same wrong thing and is corrected. The test now checks both halves, and the crossover in head start is checked for monotonicity so it reads as a threshold rather than an accident. - I asserted that concentrating contact within age groups raises R0. It does not: with equal group sizes and equal row sums the next-generation matrix has the same dominant eigenvalue either way, 5.0 against 5.0. What raises R0 is heterogeneous *activity* -- under proportionate mixing the eigenvalue is <k^2>/<k gamma> rather than <k>/gamma -- and the test now checks that closed form across three spreads, along with the core-group case where a tenth of the population sustains an R0 of 25. - The same test then asserted R0 was *below* a crude average of the groups' own reproduction numbers. It is above it, and that is the whole lesson. - At R0 = 1.1 the epidemic grows at 0.025 per unit time and needs some 550 time units merely to climb from a millionth. A fixed horizon of 400 truncated it and the final size came out short -- a run that had not finished, not a defect. The horizon now scales with the growth rate and the test asserts the epidemic actually ended before comparing. The tests lean on closed forms where they exist: the final size against the transcendental equation it solves, the peak against S = 1/R0, the SIS endemic equilibrium at 1 - 1/R0, the demographic equilibrium at S = 1/R0 with a damped oscillation on the way, the complete graph's spectral radius of n - 1 and the star's of sqrt(n - 1), extinction at (1/R0)^i0 measured over three thousand stochastic runs, and the Cori estimator inverted against a renewal process built with a known R and serial interval. Where a claim could only be established numerically -- the heterogeneous-activity eigenvalue -- it was computed independently before being written down. tests/properties/epidemiology_props.rs adds 12 property tests, including a check of the network threshold against a direct power iteration on the adjacency matrix and a comparison of the integrated final size with the implicit solution over random parameters. 3811 lib tests and 311 property tests pass in debug; clippy is clean under --all-targets -D warnings, and the module checks on nightly-2025-11-21. Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Roadmap section 18, second module. The adaptive integrator moves from `epidemiology.rs` up to `biophysics/mod.rs` so both modules share one copy. population.rs -- logistic, Gompertz and Richards growth in closed form, the Allee effect, Lotka-Volterra with its conserved quantity, Rosenzweig- MacArthur with the enrichment threshold, two-species competition with its four outcomes, Leslie matrices with the stable age distribution and Euler-Lotka, Ricker and Beverton-Holt maps with a bifurcation diagram, the Levins metapopulation, Wright-Fisher and Moran drift with the exact fixation probability, heterozygosity decay, Hardy-Weinberg with its chi-square test, one-locus selection with the balanced polymorphism, mutation-selection balance, Hamilton's rule, the Price equation, coalescent times, Watterson's theta, nucleotide diversity, Tajima's D and Fst. Defects found while writing the tests: - The shared integrator reported a numerical breakdown when it had in fact finished. Accumulated time overshoots the end by a rounding residue -- of order 1e-14 on a span of 235 -- leaving the loop condition true and a final "step" smaller than the working precision. An interval too short to bother with and an error controller forced into an impossible step are different conditions and are now distinguished. - `richards` was written with `e^(-r nu t)`, which solves the tidier-looking dN/dt = r N (1 - (N/K)^nu) and destroys the Gompertz limit: the effective rate is r nu, so letting nu fall at fixed r freezes the curve at its initial value instead of approaching anything. It now uses `e^(-r t)`, solving dN/dt = (r/nu) N (1 - (N/K)^nu), where r is the intrinsic rate in both limits and the family is a genuine interpolation rather than two special cases with a gap between them. Defects in the tests themselves, recorded rather than quietly patched: - I asserted a four-cycle in the Ricker map at r = 2.5. It begins at 2.526; 2.5 is still inside the two-cycle window. The windows narrow fast and guessing at their edges gets them wrong, so the cascade was measured before being asserted. - Locating the bifurcations by counting attractor points then failed at the 8-to-16 split, giving a Feigenbaum ratio of 2750. The reason is physical: convergence at a bifurcation is algebraic rather than geometric, so no finite transient settles and a 4-cycle is indistinguishable from an 8-cycle just below the split. They are now found from the multiplier of the cycle, which has no such problem, and the second ratio comes out 4.59 against Feigenbaum's 4.669. - My doc comment claimed a recessive allele at mu = 1e-6, s = 0.1 sits three hundred times commoner than one with h = 0.1. The factor is thirty-one -- sqrt(s/mu) h -- and the quoted frequency was wrong by a decade as well. Both are corrected and the test now checks the closed form rather than a round number picked by eye. The tests lean on closed forms throughout: each growth law is checked against a central difference of its own differential equation rather than against a remembered curve; Beverton-Holt against its exact solution at every step; the Leslie eigenvalue against the Euler-Lotka root, two entirely separate computations; Moran fixation against 4,000 simulations per case; Wright-Fisher against both the martingale property and the closed-form variance p(1-p)(1-(1-1/2N)^t); heterozygosity decay against the same simulation; coalescent intervals against 4N/(k(k-1)) over 20,000 runs; and the Price equation as the exact identity it is. tests/properties/population_props.rs adds 17 property tests over random parameters, including the three growth laws against their own equations, the competition criterion against the integrated outcome, and Fst, Hardy-Weinberg and the Price equation as identities. 3834 lib tests and 328 property tests pass in debug; clippy is clean under --all-targets -D warnings, and the module checks on nightly-2025-11-21. Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Roadmap section 18, third module. seq_align.rs carries Needleman-Wunsch, Smith-Waterman, Gotoh affine gaps, banded alignment and Hirschberg's linear-space traceback; the BLOSUM62 and PAM250 matrices; reverse complement, transcription, the genetic code, translation and ORF finding; Wallace and nearest-neighbour melting temperatures; Hamming, p, Jukes-Cantor and Kimura distances; a k-mer index, minimizers and a suffix-array search; centre-star multiple alignment with profiles, consensus and PSSM scoring; and a de Bruijn assembler. Every algorithm that returns an alignment also returns a score, and the two are checked against each other by rescoring: a dynamic program that reports a maximum it did not reach is the commonest way for one of these to be wrong, and a score-only comparison cannot see it. `alignment_score` and `alignment_score_affine` exist for that purpose and are public, since the same check is worth having outside the tests. Four of the algorithms compute the same optimum by different means -- the quadratic table, Hirschberg's linear-space recursion, Gotoh's three tables with a free opening cost, and a band wide enough to hold the whole table -- and the property tests require all four to agree on random inputs. Any disagreement is a defect in one of them. Defects in the tests themselves, recorded rather than quietly patched: - I asserted that the affine and linear optima coincide whenever the affine alignment has no gaps. They need not: with a costly opening the affine optimum takes mismatches where the linear one buys gaps, so the two optima legitimately differ. What is true, and is now checked, is that the same gapless alignment scores identically under either model. - I asserted that identity scores highest for every letter of the BLOSUM and PAM alphabets. That is a statement about *residues*: B, Z and X are ambiguity codes whose scores are averages, and BLOSUM62 gives X/A zero against X/X of minus one. The property now runs over the twenty standard residues, which is what it was ever about. - A minimizer assertion I wrote ended in `|| true` and could not fail. It is replaced by the real statement -- that each selected k-mer is minimal over some window containing it -- computed from the text directly rather than through the function under test. The tests lean on structure rather than stored numbers: the genetic code is checked for fourfold degeneracy in the third position and for methionine and tryptophan being the only single-codon residues, rather than against a table; the reverse complement for being an involution that preserves GC; the k-mer index and suffix-array search against a naive scan, which is the only thing there that is obviously right; Jukes-Cantor by inverting its own closed form and by reducing to Kimura at the one-to-two transition ratio it implicitly assumes; and the assembler for never inventing a k-mer no read holds while losing none that they do. The two melting-temperature models are checked to *disagree*: Wallace ignores stacking, which is a small error at fourteen bases and thirty degrees at sixty. A test that only checked their agreement would have been asserting something false. 3852 lib tests and 339 property tests pass in debug; clippy is clean under --all-targets -D warnings, and the module checks on nightly-2025-11-21. Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Roadmap section 18, fourth module. `PhyloTree` stores a parent index and
a branch length per node, which makes root-walking and MRCA direct at the
cost of making "children of" a search. On top of it: Newick parse and
emit, leaves, height, total length, patristic distance, splits as rooted
clades and as unrooted bipartitions, and Robinson-Foulds. Then UPGMA and
neighbour joining, Fitch parsimony, Felsenstein pruning under JC69,
column bootstrap support, birth-death simulation with the extinct
lineages pruned away, the gamma statistic and lineage-through-time.
Two design points are worth stating because they decide what the tests
can assert. `bipartitions()` exists separately from `splits()` because a
neighbour-joining root is an artefact: on a four-taxon tree {A,B} and
{C,D} name the same branch, and comparing them as rooted clades would
report two replicates that found the same tree as disagreeing. Bootstrap
support is therefore computed on bipartitions, which is also the
convention support values are reported under. And `birth_death_tree`
stops at the first event *after* the target count is reached rather than
at the branching that reaches it: stopping on the branching leaves the
last internode interval exactly zero, which biases the gamma statistic
by about +sqrt(3/n) -- measurably, +0.50 at 40 tips before the fix and
+0.03 after, against a theoretical mean of zero.
The strongest tests are the ones with an exact answer to check against:
- Felsenstein pruning is compared against enumerating all 4^internal
ancestral assignments on a five-tip tree. Pruning is a rearrangement
of that sum, so the two must agree to rounding, and they do.
- The maximum-likelihood branch length for a pair of sequences is found
by scanning and lands on the closed-form Jukes-Cantor distance.
- Neighbour joining is checked by inverting the patristic map: given
distances that came from a tree it returns that tree's distances
exactly. UPGMA does the same on ultrametric input.
- Patristic distances satisfy the four-point condition -- two of the
three pairings equal, the third no larger -- which is the defining
property of a tree metric.
- Gamma on pure-birth trees has mean 0.03 +- 0.05 and standard deviation
1.00 over 250 replicates, which is its null distribution. That checks
the simulator and the statistic against each other; either being wrong
breaks it.
- Unequal rates fool UPGMA where neighbour joining holds: with two fast
and two slow tips the slow pair is closest in the matrix without being
related, average linkage joins them, and the Q correction does not.
Defects found and fixed while writing the tests, in my own text rather
than in the code:
- The `gamma_statistic` and `birth_death_tree` docs said extinction
pushes gamma negative. It pushes it positive: near the present
lineages have not had time to die, so the reconstructed tree's nodes
crowd toward the tips. Measured +1.5 at mu/lambda = 0.5 against 0.0
for pure birth. The bias runs opposite to the slowdown test, which is
why a significantly negative gamma is read as conservative evidence.
- A first property test asserted that a Newick round trip preserves leaf
*index* order. It does not, and should not: Newick encodes the tree,
not the node numbering. Rewritten to compare distances by label.
- A hand-checked tree's total length was written as 8 where the string
says 1+1+2+3.
- Two test sequences differed at 15 of 20 sites, which is the point
where the Jukes-Cantor correction stops being finite; the module
correctly refused them and the test, not the code, was wrong.
3889 lib tests and 355 property tests pass in debug; clippy is clean
under --all-targets -D warnings; checked on nightly-2025-11-21.
Co-Authored-By: Claude <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Roadmap section 18, fifth module, which completes the section. Conductance models -- Hodgkin-Huxley, Morris-Lecar in both its type I and type II parameterisations, FitzHugh-Nagumo -- alongside the integrate-and-fire family: LIF with an exact F-I curve to check the simulation against, Izhikevich with the five published presets, and AdEx with both adaptation currents. Then spike train statistics, Poisson trains, PSTH and raster, the spike-triggered average, a von Mises tuning fit, exponential and alpha synapses, the STDP window and its all-to-all pairing, Izhikevich's random network, Hopfield storage and recall, Wilson-Cowan, the passive cable, and the drift-diffusion decision process with its gambler's-ruin accuracy. `nernst_potential` and `goldman_potential` already exist in `biophysics` and are not repeated; the module header says so. Three defects the tests found in the code, all of them things that would have returned plausible numbers rather than failing: - The sealed-end boundary condition in `cable_equation_1d` used the one-sided difference V[n-1] = V[n-2]. That imposes a zero gradient only to first order, and the boundary decides the convergence rate of the whole solution: halving the spacing cut the error by 2.03 where a second-order scheme must cut it by 4. Replaced with a ghost node reflected through the end, giving 2 V[n-2] - (2+k) V[n-1] = 0. The ratio is now 4.02. - `hopfield_recall` updated every unit simultaneously against the previous state. That has no Lyapunov property -- the energy can rise and the network can settle into a two-cycle between two states, neither of them stored -- and the doc comment claimed the energy never increases, which is true only of sequential updates. Switched to sweeping the units in index order, which makes the claim true and the recall convergent. - `reaction_time_ddm` capped each trial at a step budget divided by the trial count. A decision time has a long tail, so that threw away exactly the slow trials the distribution is about, and errored on an ordinary run of 2000 trials. The budget is now shared across trials. Two documentation errors, both about which way an effect runs: - `hodgkin_huxley` did not say that a strongly *hyperpolarising* current is what the fixed step cannot follow. beta_m grows exponentially as the membrane hyperpolarises, so below about -25 uA/cm^2 it reaches thousands per millisecond and 0.01 ms is no longer stable; the function reports the breakdown, and there is now a test that it does. Depolarising currents integrate cleanly at 500 uA/cm^2 and simply block. - The F-I curve's doc said the rate jumps at "the rheobase" without distinguishing the two rheobases. A 2.24 uA/cm^2 step makes the model fire once and then sit still; repetitive firing needs 6.3, where the rate is already 50 Hz. Both numbers are now tested, and the fact that they differ is its own test. The tests that carry the most weight are the ones with an exact answer: - The simulated LIF rate matches the closed-form 1/(t_ref + tau ln(...)) to within 2% across five currents and randomised parameters. - Simulated decision accuracy matches 1/(1 + exp(-2 A a / sigma^2)) to within sampling error, and the formula is shown to depend only on the combination A a / sigma^2. - The cable solution matches cosh((L-x)/lambda)/cosh(L/lambda) and converges at second order. - The von Mises fit recovers its own generating parameters to 1e-7 over forty randomised curves, because the log form is linear and solved rather than searched. - Hodgkin-Huxley's all-or-none response: a 0.5 ms pulse at 10 uA/cm^2 fails entirely, and doubling a suprathreshold one from 20 to 40 moves the peak by 1.4 mV rather than doubling it. - Type I and type II excitability begin differently: the saddle-node parameterisation starts at 1.4 Hz where the Hopf one jumps to 7.1 Hz, from silence in both cases. - A Hopfield probe can settle in a spurious state deeper than the pattern it started from -- asserted directly, since descending the energy finds a minimum and not the right one. 3935 lib tests and 376 property tests pass in debug; clippy is clean under --all-targets -D warnings; checked on nightly-2025-11-21. Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Roadmap section 19a, first module, and a new top-level `finance/` since nothing existing is a home for it. Black-Scholes-Merton with a continuous dividend yield, the Greeks, implied volatility, put-call parity, CRR binomial and trinomial lattices with American exercise, Monte Carlo for European, Asian, barrier and lookback payoffs, Longstaff-Schwartz least squares, Merton jump diffusion, Heston by simulation, Crank-Nicolson on the log-price grid, the SVI smile with a fit, and a delta-hedging simulation. Three defects the tests found in the code: - `merton_jump_price` weighted its Poisson sum with intensity `lambda` where the risk-neutral intensity is `lambda (1 + k)`. The series still converged and still looked like a price -- it matched Black-Scholes at `lambda = 0` and produced a plausible smile -- but the discount factors no longer summed to `e^(-rT)`, so the call and put prices violated put-call parity. That is an arbitrage in a model whose whole purpose is to be free of them, and only the parity property caught it: no price-level comparison would have. - `trinomial` indexed the previous layer by an offset from the terminal width instead of the constant shift of one that the recursion actually has, and read one past the end. Every call panicked. - `implied_volatility` stopped on a price tolerance. Where vega is small that is meaningless: a call struck at 70 with the share at 100 and eighteen days to run prices identically at 5% and at 20% volatility to the last bit of a double, and the solver returned 10% -- the midpoint of its first bracket -- as though it had measured something. It now bisects on the volatility bracket, and returns `None` when vega falls below `1e-8` relative to the price, because the price does not determine a volatility there and reporting one is reporting rounding noise. One property of my own writing that turned out to be false. I had asserted that averaging two consecutive binomial step counts beats either, since the error oscillates in sign. It does oscillate -- 100 steps lands 2.0e-2 below the exact price and 101 steps 1.7e-2 above it, which is now a test -- but across randomised parameters the averaging helped on only 11 of 40 draws, so the phase is not predictable and the claim is not a theorem. Replaced with two properties that are: - The CRR lattice is arbitrage-free at *every* step count, because its up-probability is chosen to make the price an exact martingale. Parity holds to rounding even at seven steps, where the price is nowhere near the continuous answer. - The trinomial is arbitrage-free only in the limit. Matching the first two moments of the log price makes the price a martingale to O(dt^2), so its parity residual is real at coarse steps -- 2.2e-3 on a two-and-a-half year option at seven steps -- and falls as one over the square of the step count, reaching 4.1e-8 by sixteen hundred. The test asserts the rate, which says the violation is a discretisation artefact rather than a defect, and the doc now says which lattice to use when an exactly consistent call and put matter more than smooth convergence. The strongest tests are the model-free identities and the degenerate cases, since both have exact targets: - Put-call parity across 600 randomised parameter sets, and the bounds arbitrage would close. - Homogeneity: doubling the spot and strike together doubles the price. - The symmetry C(S,K,r,q) = P(K,S,q,r), exact to 1e-10. - Merton with no jumps reproduces Black-Scholes to 1e-13, and Heston with no volatility of volatility reproduces it within its own standard error. - A knock-in and a knock-out priced on identical paths sum to the barrier-free price to 1e-9, since every path pays into exactly one. - Crank-Nicolson converges at second order: errors of 3.4e-3, 8.6e-4 and 2.1e-4 as the grid doubles twice, ratios of 4.00 and 4.01. - The Greeks match Richardson-extrapolated finite differences of the price they differentiate, over randomised parameters. - Delta hedging's residual risk falls as one over the square root of the rebalance count, and hedging at 20% into a 30% world loses money well outside the sampling error. 3967 lib tests and 396 property tests pass in debug; clippy is clean under --all-targets -D warnings; checked on nightly-2025-11-21. Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Roadmap section 19a, second module. Compounding conventions and conversion between them, discount factors, NPV, IRR and XIRR, bond pricing and yield solving, Macaulay and modified duration, convexity, zero-curve bootstrapping with linear interpolation in the zero rate, forward rates, Nelson-Siegel with a fit, Vasicek and CIR bond prices, and level-payment amortisation. Four defects the tests found in the code: - `ns_fit` fitted the basis `[1, slope, slope - decay]`, which is linear in the three coefficients and so solves exactly, but then converted back with `b1 = beta1 - beta2`. Matching coefficients gives `b1 = beta1` directly -- the `slope - decay` vector already carries the `+ b2 slope` term. The curve fit was off by 280 basis points while every internal residual looked fine, because the *fit* was right and only the reported parameters were wrong. - `cir_bond_price` at zero volatility evaluated `1^infinity`. Its `A(t)` factor is a base tending to one raised to `2 kappa theta / sigma^2`; with `sigma = 0` the exponent is infinite, the base is exactly 1.0, and IEEE resolves `powf(1.0, inf)` to 1.0 -- silently dropping the whole factor and returning `e^(-B r0)`. That gave 0.9464 where the deterministic answer is 0.8339. Now special-cased to the Vasicek limit, with a doc note that the closed form is dominated by rounding below about `sigma = 1e-5`. - `bracket_and_solve` pinned its lower end near -100%. Discounting a hundred-period schedule at -99.99% raises 1e-4 to the hundredth power and overflows, so every bond with a deeply negative yield failed to solve at all. It now walks down toward -1 and keeps the last point where the function is still finite. - `nelson_siegel` computed `(1 - e^-x)/x` directly. That cancels catastrophically for small `t/tau`: at 1e-10 it keeps about six digits, and the short-end limit was wrong in the seventh. Switched to `-exp_m1(-x)/x`, which is accurate to the denormals. One claim of my own that was false. I asserted a short-rate bond price is always below one. It is for CIR, whose square-root diffusion keeps the rate non-negative. It is not for Vasicek: with weak mean reversion the convexity term `sigma^2/(2 kappa^2)` can exceed `theta` outright, the long-run yield goes negative, and the bond is worth more than the unit it pays. That is the Gaussian model's known feature, so it is now asserted where it applies and demonstrated where it does not. The tests lean on identities and on substituting answers back, since almost everything here is defined as the solution to an equation: - A rate converted between five compounding conventions returns exactly, and discounts identically at every horizon, not only at one year. - IRR, XIRR and yield-to-maturity each zero the equation they were solved from, over randomised cashflows. - A bond prices at par exactly when its coupon equals its yield, at every rate and term. - Bootstrapping recovers the curve its bonds were priced from to 1e-12, and reprices them to 1e-10. - A forward rate makes rolling equal to holding to 1e-12, which is the no-arbitrage identity it is defined by. - Duration and convexity match Richardson-extrapolated differences of the price, and a zero-coupon bond's duration is exactly its maturity. - Both short-rate models with no diffusion reproduce the deterministic integral `exp(-int r)` to 1e-13 and each other. - Cashflows changing sign twice get no internal rate of return: both 10% and 20% zero the value of one such series, and reporting either as *the* return would be a mistake the code refuses to make. - A mortgage at 5% over thirty years does not repay more principal than interest until period 195 of 360, and a higher rate pushes that later. 3988 lib tests and 412 property tests pass in debug; clippy is clean under --all-targets -D warnings; checked on nightly-2025-11-21. Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Roadmap section 19a, third and fourth modules, which complete the section. portfolio.rs: simple and log returns, the Markowitz frontier in closed form, minimum-variance and tangency weights, risk parity and risk contributions, Sharpe, Sortino, max drawdown, Calmar, information ratio, CAPM beta, and both Kelly fractions. risk.rs: historical, parametric and Cornish-Fisher value at risk, expected shortfall, a GARCH(1,1) one-step forecast, a moving-average crossover backtest, and Kupiec's coverage test. One real defect, and it was invisible from the outside. `risk_parity_weights` iterated `w_i <- w_i / (C w)_i`. At rest that gives `(C w)_i` equal across assets -- which is the *minimum-variance* condition, not equal risk contribution -- so it returned the minimum-variance weights under another name, and every plausibility check on them passed. The damped update `w_i <- sqrt(w_i / (C w)_i)` has the right fixed point: `w_i (C w)_i` is then the same constant for every asset. The contributions are now 1/3 each on a three-asset problem where they used to be 0.73, 0.16 and 0.11, which is what gave it away. Two things I had to get right about Cornish-Fisher, neither of which is a bug so much as a limit that had to be measured before it could be documented: - Its kurtosis term carries the factor `z^3 - 3z`, which is zero at `z = -sqrt(3)`, an alpha of about 4.2%. So the same sample gets the correction applied one way at 1% and the other way at 5%. On a uniform sample the 1% estimate moves from 0.673 to 0.591 against a true 0.490 -- the right direction -- while the 5% estimate moves from 0.476 to 0.483, the wrong one. - It is asymptotic, not convergent. With a skew of -0.39 and an excess kurtosis of 0.85 it improves on the Gaussian fit, moving a 1% VaR from 0.0197 to 0.0234 against a historical 0.0300. With a skew of -4.6 and an excess kurtosis of 33.8 it returns 0.0729 where the sample's own quantile is 0.0309. Both regimes are now tested, the second as a known failure rather than papered over, and the doc says the moments have to be looked at before the number is trusted. The guard was also replaced. It had been an arbitrary bound on how far the correction may move the quantile; it is now the standard validity condition -- the corrected quantile must be increasing in z, since a quantile function that decreases is not one. The tests lean on perturbation and on coherence, both of which check a solution without trusting the formula it came from: - The minimum-variance portfolio's variance rises in every budget- preserving direction, and the excess scales exactly quadratically with the step, which is what a minimum looks like. - Every frontier point has strictly less variance than any portfolio with the same expected return, tested by perturbing along directions orthogonalised against both constraints. - The tangency portfolio has the highest Sharpe ratio in its neighbourhood. - Expected shortfall is positively homogeneous, translation-equivariant and subadditive on every sample. Value at risk satisfies the first two and fails the third: two independent bonds each defaulting in 4% of scenarios have a *negative* 95% VaR individually and a VaR of 0.5 combined, so diversifying raised the measured risk. That is the reason the regulatory measure changed, and it is now a test. - A beta built into a series comes back out of it to 1e-9, alpha included. - Kelly is the vertex of the growth rate: growth falls on both sides and returns to exactly the risk-free rate at twice the fraction. - The crossover backtest on a monotone series matches buy-and-hold over the period it was invested, exactly -- any lookahead would beat it. - Kupiec's statistic is exactly zero at the expected breach count and grows in both directions from it. A note on CI wall time, which I have raised before and which now has a cause: the `test` job runs the suite twice, once plainly and once instrumented under `cargo llvm-cov --summary-only`. That is why it sits near eighteen minutes while the local debug suite takes three and a half. Splitting coverage into its own job would roughly halve it, but that is a workflow change rather than a module, so I have left it alone. 4017 lib tests and 424 property tests pass in debug; clippy is clean under --all-targets -D warnings; checked on nightly-2025-11-21. Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Roadmap section 19b, first part, under the existing astrophysics/ directory rather than a new astro/. Kepler's equation solved for elliptic and hyperbolic orbits, conversions among the true, eccentric and mean anomalies in both directions, state_from_elements as the inverse of the existing OrbitalElements::from_state_vectors, two-body propagation by Lagrange's f and g functions, the orbital period and vis-viva. The element set, the state-to-elements conversion and the geometric quantities already live in orbital_elements.rs and are reused rather than repeated. Three defects the tests found, and the first is the kind that only a conservation law catches: - The hyperbolic branch of the Lagrange coefficients had two sign errors. Substituting E = i H and sqrt(a) = i sqrt(-a) into the elliptic form cancels the imaginary units in `f` and `g_dot` but not in `g` or `f_dot`, and I had carried the elliptic signs through. The trajectory came out smooth and plausible -- it receded, it curved the right way -- while the specific energy drifted by 4e-4 over a hundred seconds and by half its own value over a thousand. Energy conservation is now 1e-15. - `kepler_solve_elliptic` wrapped its converged root to [0, 2pi). At a mean anomaly of zero Newton lands on zero from either side, and an undershoot of one ulp came back as a full turn. The root lies in the same revolution as the mean anomaly, so it is clamped rather than wrapped. - `kepler_solve_hyperbolic` seeded small mean anomalies with the textbook `M/(e-1)`. That diverges as the orbit approaches parabolic: at e = 1.001 it puts the first guess at four hundred, where cosh overflows and the iteration has no derivative left. Replaced with `asinh(M/e)`, which inverts the leading term and is bounded everywhere. Two corrections to my own documentation: - `vis_viva` rejected an infinite semi-major axis, which is exactly the parabolic case `1/a = 0` where it should return escape speed. The doc claimed the formula covers all three conics while the guard refused one of them. - The same function's error condition said "outside the orbit". The real boundary is `r > 2a`, where the kinetic energy runs out. For a bound orbit that reaches past apoapsis, so between `a(1+e)` and `2a` the formula answers with the speed a body of that energy *would* have, which is not a speed anything reaches. The doc now says so and the property test pins the boundary at 2a rather than at apoapsis. The tests are built on the three kinds of invariant orbital mechanics supplies: - Inverse pairs. Kepler's equation solved and read forward composes to the identity to 1e-11 at eccentricities up to 0.9999; the three anomalies cycle back to 1e-9; elements and state vectors invert each other to 1e-7 in the angles and 1e-8 in the state. - Conserved quantities. Energy, the angular momentum *vector* and the eccentricity vector are all unchanged by propagation, over spans up to three periods forward and back. The eccentricity vector is the one that pins the orbit's orientation within its plane, which the other two do not. - Group structure. Propagating by t1 then t2 equals propagating by t1 + t2, and -t undoes t. Two-body motion is a one-parameter flow and a propagator that is not one is wrong somewhere. Beyond those: propagation is checked against an entirely separate route -- convert to elements, add n dt to the mean anomaly, convert back -- and the two agree to 1e-6 of the radius. Kepler's second law appears as an inequality that holds everywhere on the orbit: between periapsis and apoapsis the true anomaly leads the mean, and past apoapsis it lags. A parabolic orbit has neither an elliptic nor a hyperbolic anomaly and is refused rather than forced into the wrong branch. 4031 lib tests and 437 property tests pass in debug; clippy is clean under --all-targets -D warnings; checked on nightly-2025-11-21. Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Roadmap section 19b, second part, completing the section. lambert.rs solves the boundary-value problem by universal variables and builds porkchop grids of departure characteristic energy from it. maneuvers.rs adds what propulsion.rs, tidal.rs and lagrange.rs do not already cover: combined burns, the sphere of influence, patched-conic escape, gravity assist deflection, the Oberth effect, J2 nodal drift, sun-synchronous inclination and ground tracks. The Hohmann and bi-elliptic transfers, the plane-change delta-v, Tsiolkovsky, staging, the Roche limit and the Hill radius already existed and are referenced rather than repeated. Two numerical defects, both in the same function and both invisible until the geometry pushed on them: - `stumpff_c` evaluated its positive branch as `(1 - cos u)/u^2`. At the single-revolution boundary `z = 4 pi^2` the cosine is within an ulp of one, the subtraction keeps no digits, and the result comes back as zero or negative -- which made the flight time infinite and the whole upper bracket unusable. Rewritten as `2 sin^2(u/2)/z`, identical in exact arithmetic and accurate at both ends: the sine is *small* there rather than large, so squaring it loses nothing. Verified to follow the expected quadratic vanishing to a part in a thousand at 1e-10 from the boundary. - The lower bracket walked z downward looking for a bound that does not exist, doubling to -1e12 and giving up. It now stops at the first z that is either fast enough or has no positive chord, the latter being a valid lower bracket since the bisection treats a missing solution as "too fast". One documentation claim of mine was wrong. I had written that a flight time shorter than the minimum-energy transfer has no solution. It does: within one revolution a transfer exists for every positive flight time, and hurrying simply costs more without limit. The minimum-energy transfer is a particular duration, not a floor on one -- which is now a test, with the departure speed scanned across two orders of magnitude of flight time and its minimum shown to be interior, falling before and rising after. Four of my own test claims were also wrong and are recorded rather than quietly fixed: that the saving from a combined burn grows with the plane change (it depends on which order the separate burns are taken in, so only the inequality is a theorem); that the sphere of influence is *inside* the equal-force radius (it is three and a half times outside it, 924,000 km against 259,000 -- a smaller exponent on a ratio below one gives a larger answer); that a 500 km orbit at 45 degrees drifts 4.9 degrees a day (5.4); and that a mirror symmetry through the equator flips the prograde flag (it does not -- reflecting z leaves the z component of `r_a x r_b` unchanged, so the *same* flag gives the mirrored solution, which is now the property tested). Lambert reproduces Vallado's example 7-5 to six decimal places. The property that carries the most weight is independent of that: given an arc generated by `propagate_kepler`, the solver recovers the very velocities that generated it, and flying its answer lands on the target. Over three hundred randomised geometries the departure velocity agrees to a part in 1e8 or better. That last figure is a conditioning limit, not a tolerance chosen for convenience. The velocities come out as `(r2 - f r1)/g`, and as the transfer angle approaches pi that numerator is a difference of two nearly equal vectors. The worst residual over three thousand draws was 1.2e-8, and its transfer angle was 179.99 degrees. Exactly pi is refused, since the plane is then undefined and infinitely many orbits connect the points; the approach to it is merely imprecise, and the docs now say so. Other properties: a porkchop cell reproduces the Lambert solution it came from exactly; a combined burn never exceeds either sequential ordering and is symmetric in the two speeds and even in the angle; a flyby's turn depends only on the combination `r_p v^2 / mu`, which is checked by scaling two of the three and finding the turn unmoved; the Oberth gain matches `v dv + dv^2/2` exactly and periapsis beats apoapsis by more than a factor of four on a 0.7-eccentricity ellipse; the J2 drift is westward prograde, exactly zero at the pole and eastward retrograde, and the sun-synchronous solver inverts it to a part in 1e9; and a ground track's latitude is bounded by the orbit's inclination and attains it, while successive ascending nodes walk west by exactly one body rotation per orbital period. 4047 lib tests and 448 property tests pass in debug; clippy is clean under --all-targets -D warnings; checked on nightly-2025-11-21. Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Roadmap section 19b, final part. Two modules under the existing astrophysics/ directory. time_systems.rs: julian_date and jd_to_calendar, gmst and local_sidereal, tle_epoch_to_jd, and the J2000/JULIAN_CENTURY constants. coords.rs: equatorial_to_horizontal and its inverse, the ecliptic pair, mean_obliquity and precession_approx, sun_position_approx, moon_position_approx, planet_position_low_precision over Standish's elements, rise_set_times, and tle_parse_lite. The TLE reader parses and checksums only. A TLE's numbers are *defined* by SGP4 -- they are mean elements in Brouwer's theory, not osculating ones -- so converting them to a state vector without SGP4 would give an answer that is wrong by kilometres while looking entirely reasonable. The doc comment says so rather than leaving the omission to be guessed at. Two defects the tests found: - jd_to_calendar used the textbook Julian-calendar branch below JD 2299161 (the 1582 reform) while julian_date is proleptic Gregorian throughout, so the two stopped inverting each other before the reform: 1 January -4712 went out as JD 38 and came back as 8 February. Made proleptic on both sides; the historian's convention is a different function, not this one. - equatorial_to_horizontal took the altitude with asin. Near the zenith the argument is within an ulp of one, where asin has a square root's conditioning -- 1e-16 in becomes 1.5e-8 out. Now atan2(up, hypot(south, east)), which is well conditioned everywhere, and the round trip closes to 1e-15 at the pole. Also corrected in my own tests: precession in right ascension is m = 46.12"/yr, not the 50.29" of general precession in longitude; those are different quantities. 18 unit tests and 14 property tests. Suite is 4,065 lib + 462 property tests, green in debug, clippy clean under --all-targets -D warnings, and checked on nightly-2025-11-21. Still outstanding for you, unchanged from the last two sessions: the CI test job runs the suite twice (cargo test, then cargo llvm-cov rebuilds instrumented and reruns it), which is the 14-18 minute wall time; splitting coverage into its own job is a workflow change I have left alone. And PR #4 now spans sessions 4-37. Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Roadmap section 19c, first part. New fem/ module; fem1d.rs holds fem_1d_poisson, fem_1d_general and fem_1d_quadratic, the Bc enum, a Fem1dSolution wrapper that can evaluate between nodes, the L2 and H1 error norms, and convergence_rate. Assembly is five-point Gauss per element into a symmetric banded matrix, solved by L D L^T without pivoting -- which is unconditionally stable while the problem is coercive, and a reaction term negative enough to lose coercivity is a genuinely singular operator rather than something to pivot around. Flux conditions use the outward normal at both ends, so the same Neumann value means the same physical thing on the left and the right, and Robin is p du/dn + alpha u = g in the same convention, which keeps the matrix symmetric. Dirichlet data is eliminated symmetrically rather than by zeroing a row. A pure flux problem is reported as Singular, detected by the exact criterion rather than a threshold: the constant function is in the kernel exactly when every row of the assembled matrix sums to zero. A reaction term, a Dirichlet end or a nonzero Robin coefficient each independently removes it, and a Robin end with a zero coefficient is a flux condition that pins nothing -- the boundary case the check has to get right instead of treating "Robin" as a keyword. Two things the tests turned up, both now documented rather than papered over: - Nodal exactness for Poisson is exact only up to the quadrature of the *load*. With a transcendental f the residual nodal error falls off as the five-point rule does, around h^11, not as the h^2 of the solution -- so the two are separated by measuring the rate rather than by loosening a tolerance, and an assembly error would still show up as second order. - Linear elements have a constant derivative per element, so the stiffness quadrature reproduces the element *average* of p exactly. The discrete bilinear form therefore still agrees with the true one on the element space, and the discrete solution is the exact a-orthogonal projection rather than an approximation of one. The property test asserts the Pythagoras identity that follows, which holds to 1e-17, instead of the Cea inequality it implies -- an equality cannot be satisfied by accident, and the earlier inequality form was nearly non-strict for exactly this reason. Quadratic elements are exact at element vertices and merely third-order at the midsides, for a reason worth stating: the vertex Green's function is piecewise linear and lies in the space, while the midside one kinks inside an element and does not. Both halves are asserted. 13 unit tests and 18 property tests, the latter covering Galerkin orthogonality, Ritz minimisation and its exact quadratic excess, nodal exactness, the patch test at both degrees, superposition, the discrete maximum principle, nested-refinement energy monotonicity, the exact conservation law from testing against the constant, reflection symmetry, the sharp Poincare constant, and the h^2/h^3 convergence orders that identify the spaces. Suite is 4,078 lib + 480 property tests, green in debug, clippy clean under --all-targets -D warnings, checked on nightly-2025-11-21. Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Roadmap section 19c, second part. fem2d.rs holds FemMesh2 with its three generators (rect, disk, from_delaunay), refine_uniform, quality_min_angle, the assembled stiffness and mass matrices, element_gradient, dirichlet_energy, barycentric interpolate, and fem_2d_poisson / fem_2d_reaction_diffusion solved by Jacobi-preconditioned conjugate gradients. The linear triangle needs no quadrature for the stiffness term at all -- the shape function gradients are constant, so the element matrix is the gradient product times the area, exactly. Coefficients and the source are sampled at the centroid, a one-point rule of the same order as the element itself. FemMesh2::new orients every triangle counterclockwise rather than rejecting a clockwise one: the sign of the area is a labelling convention, while a zero area is not and is refused. The boundary is derived from the edge counts, and an edge in three triangles is reported as NotManifold. The disk generator puts 6k points on ring k so the arc spacing tracks the radial spacing, and merges consecutive rings by angle, which keeps the triangles from going thin at the rim the way a fixed point count per ring would. Two things worth recording: - Uniform refinement leaves quality_min_angle *exactly* unchanged, because the four children of a triangle are all similar to their parent. The property test asserts equality to 1e-13 rather than a bound, since a quality measure that drifts under midpoint refinement is measuring something other than shape. - The off-diagonal stiffness entry for an edge is minus half the cotangent of the opposite angle. That single identity is why the Delaunay condition and the discrete maximum principle are the same statement, so the test checks the cotangent formula directly, one triangle at a time, and asserts that the entry turns positive exactly when the opposite angle turns obtuse. One test tolerance was replaced rather than loosened. The patch test's gradient check failed on a sliver in a Delaunay mesh at 1e-8. The gradient of a linear field amplifies a nodal error by the sum of the shape function gradient magnitudes, which on a sliver is large -- that is what makes slivers bad. The bound asserted is now that amplification itself, read off by differentiating each shape function's own indicator vector, which is both sharper and a statement about the method. 10 unit tests and 13 property tests: Euler's formula on all three generators, conformity and boundary-cycle structure, orientation, exact area and shape preservation under refinement, the cotangent identity, zero stiffness row sums and the mass matrix totalling the area, the patch test, Ritz minimisation with its exact quadratic excess, energy monotonicity under refinement, superposition, the discrete maximum principle on a nonobtuse mesh, rotation invariance, the inverse-square domain scaling, and second-order convergence. Suite is 4,088 lib + 493 property tests, green in debug, clippy clean under --all-targets -D warnings, checked on nightly-2025-11-21. Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Roadmap section 19c, third part. fem_2d_helmholtz solves
-lap u - k^2 u = f; fem_eigenvalues_drum and fem_eigenmodes_drum solve
the generalised problem K phi = lambda M phi over the interior nodes.
Helmholtz needs a different solver from Poisson for a structural
reason rather than a numerical one: once k^2 passes the first Dirichlet
eigenvalue the operator stops being positive definite, and conjugate
gradients is a minimisation method with nothing left to minimise. A
dense LU is used instead, and the O(n^3) cost is documented rather than
hidden. At an eigenvalue the operator is genuinely singular -- that is
resonance, not a numerical accident -- and is reported as such.
The eigenproblem goes through the Cholesky factor of the mass matrix,
which turns it into a standard symmetric one. The consistent mass
matrix is used rather than a lumped one on purpose: lumping shifts the
eigenvalues downwards, and it is precisely their being *upper* bounds
that makes them useful. Every discrete eigenvalue is a Rayleigh
quotient minimised over a subspace of the true admissible space, so it
cannot fall below the true one, and a nested refinement can only lower
it. Both halves are asserted.
Validation is against analytic spectra rather than against itself: the
unit square's pi^2(m^2 + n^2), and the circular membrane's Bessel
zeros taken from the crate's own bessel_j_zeros, which ties the solver
to the Part 3 membrane.
Two things the tests turned up:
- I had the disk's mode order wrong. j_{2,1} = 5.136 comes in below
j_{0,2} = 5.520, so the fourth and fifth modes of a circular drum are
the doubly degenerate two-nodal-diameter pair and the second radially
symmetric mode is only sixth. The test now asserts that ordering
explicitly rather than assuming it.
- The doubled eigenvalue 5 pi^2 of the square splits on the rectangle
mesh, because every cell is cut along the same diagonal and so the
mesh is not symmetric under exchanging x and y. That is a mesh
artefact of the same order as the discretisation error, and the test
now asserts both that it happens and that it stays small.
The two drum tests were 12.5 s of the debug suite. Replacing the
absolute tolerance on a 16x16 mesh with the h^2 error-ratio check
across the 6x6 and 12x12 meshes, and taking the disk to eight rings,
brought that to 2.8 s while making the assertion stronger -- a rate
identifies the space where a tolerance on one mesh does not.
6 unit tests and 7 property tests, the latter covering Betti
reciprocity (which holds for the indefinite operator too, since
symmetry has nothing to do with definiteness), the simple pole of the
resonant response, the upper-bound and refinement-monotonicity of the
spectrum, the inverse-square domain scaling and rotation invariance,
Courant's nodal domain theorem for the first two modes, and
mass-orthonormality with the Rayleigh quotient returning its own
eigenvalue.
Suite is 4,094 lib + 500 property tests, green in debug, clippy clean
under --all-targets -D warnings, checked on nightly-2025-11-21.
Co-Authored-By: Claude <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Roadmap section 19c, fourth part. fem_2d_elasticity_plane_stress on the constant-strain triangle, element_strain / element_stress / strain_energy / von_mises_stress, and fem_2d_heat_transient marching the theta scheme. Two defects in my own earlier code, both found by the elasticity patch test failing at 1e-9 when it should have been at rounding: - apply_dirichlet wrote a bare 1.0 on the diagonal of every pinned row. That is the textbook recipe and it is wrong for any problem whose natural scale is not one: an elasticity matrix has diagonal entries of order Young's modulus, so a row of 1 among rows of 1e10 gave the assembled system a condition number of 1e10 that the physics never had. The pinned rows now carry the mean free diagonal instead, which leaves the solution identical -- the row still says u_i = g -- and removed a factor of 1300 from the patch test error. - The tolerance pcg_jacobi takes is relative to the norm of the right-hand side, and I was multiplying it by the data scale at all three call sites. With elasticity data of order 1e8 that turned a requested 1e-13 into an actual 1e-5. Passing it as the pure number it is brought the patch test to 1e-13 relative, which is where it should have been all along. The Poisson and heat solvers were affected the same way; only the small numbers involved had hidden it. Both fixes are in shared code, so the Poisson, Helmholtz, eigenvalue and heat paths all get the accuracy too. The elasticity solver checks the three rigid body motions explicitly against the constraints and reports Singular if any survives, rather than letting the solver discover it. Stress is per triangle, not per node, because the strain of a linear displacement field is constant on an element -- and the doc says why averaging that to the nodes before showing it to anyone is how a coarse mesh comes to look convincing. A plane-stress subtlety worth the paragraph it gets: equal biaxial tension has a von Mises value equal to the tension, not zero. The three-dimensional intuition that hydrostatic stress cannot yield a material does not survive into plane stress, because a state that is hydrostatic in plane has a free surface out of it and so is not hydrostatic at all. Asserted directly. The theta scheme's doc and tests carry the A-stable / L-stable distinction: for a mode too stiff to resolve, backward Euler's amplification factor tends to zero and Crank-Nicolson's tends to minus one, so the stiff mode dies under one and survives under the other while flipping sign every step. That is why a discontinuous initial condition rings under Crank-Nicolson and why the remedy is to start with backward Euler steps. 11 unit tests and 9 property tests: rigid motions in the kernel, the uniform-strain patch test with a shape-independent stress, the uniaxial, pure-shear and biaxial closed forms across random materials, frame indifference of von Mises, Clapeyron's theorem with the load and modulus scalings, exact conservation of heat for the insulated problem, the theta scheme's amplification factor reproduced to 1e-7 on a discrete eigenmode, the A-stable/L-stable contrast, and the march settling onto the Poisson solution. Suite is 4,105 lib + 509 property tests, green in debug, clippy clean under --all-targets -D warnings, checked on nightly-2025-11-21. Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Roadmap section 19c, fifth part. fdtd.rs holds fdtd_courant_check and
its two-dimensional counterpart, fdtd_1d marching the Yee scheme with a
choice of perfect conductor or first-order Mur absorbing ends, an
Fdtd1d result type carrying both fields and the conserved energy, and
photonic_crystal_bandgap_1d from the Bloch dispersion relation.
Fields are normalised to E and eta_0 H, which removes the free-space
impedance from every line of the update and, more importantly, makes
the two terms of the energy comparable -- in unnormalised units one
would be 1e5 times the other and their sum would be numerical nonsense.
Three things worth recording, two of them my own errors caught by
probing the claims before writing tests around them:
- Fdtd1d::energy had the magnetic half-steps off by one. Snapshot k of
the magnetic history holds H^{k-1/2}, so the pair straddling E^n is
h[n] and h[n+1], not h[n-1] and h[n]. With the wrong pair the
invariant held exactly only for a uniform permittivity at a Courant
number of one -- where an extra symmetry rescues it -- and drifted by
parts in 1e6 otherwise. With the right pair it is conserved to 9e-16
for a graded permittivity at any admissible Courant number, which is
what the algebra says it should be.
- The Mur update was reading the edge cell's neighbour from two steps
back rather than one. The correct first-order form needs only values
from within the step, so the cross-step history is gone entirely.
- The stability limit belongs to the *fastest* wave in the grid, not to
vacuum. A permittivity below one -- a plasma above its cutoff, or an
engineered medium -- has a phase speed above c and tightens the bound
by exactly its index. Checking only the nominal Courant number would
let such a grid through to blow up, so the check is against the
smallest permittivity present, and the tests assert the threshold is
sharp on both sides.
The conserved quantity is the one leapfrog actually has, not the
obvious sum of squares. The tests assert both halves: the leapfrog form
is constant to the last bit, and the naive form is not -- asserting the
wrong one would be asserting a tolerance rather than an invariant.
Validation is against closed forms rather than against itself: the
magic time step translates a pulse bit for bit and a Courant number of
0.6 demonstrably does not; a dielectric interface reproduces the
Fresnel amplitudes to a few percent with the sign inversion off a
denser medium; the Mur end leaks a factor of a thousand less than a
wall; and the quarter-wave stack's gaps sit at exactly the odd
multiples of its design frequency with relative width
(4/(m pi)) arcsin(|na-nb|/(na+nb)) to nine digits, while the even
multiples close.
9 unit tests and 9 property tests. Suite is 4,114 lib + 518 property
tests, green in debug, clippy clean under --all-targets -D warnings,
checked on nightly-2025-11-21. CI confirmed green on all five jobs for
6f74795 before this push.
Co-Authored-By: Claude <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Roadmap section 19c, sixth part. fdtd_2d_tm marches the transverse magnetic Yee scheme with a Berenger split-field perfectly matched layer; waveguide_cutoff_check_fdtd infers a parallel-plate guide's cutoff from the evanescent decay it shows when driven below it, and waveguide_cutoff_numerical gives the value it should find. Two departures from the roadmap's signature, both to make the thing testable rather than merely runnable: - The source is a closure over the step index, as in fdtd_1d, not a frequency. Taking the waveform is the only way to switch the drive *off*, and with a source still running the field near it is the source's own and says nothing about what the boundary reflected. Every absorption measurement here depends on that. - The result carries an envelope alongside the final field. One snapshot of a driven oscillation is whatever phase it landed on; what a steady-state calculation is for is the amplitude, and that cannot be reconstructed from a single frame. The layer depth is per axis. A waveguide needs its ends absorbed and its plates conducting: absorbing the plates would stop it being a waveguide, and leaving the ends conducting lets the switch-on transient rattle around forever and swamp the field being measured. That was not a hypothetical -- it is what the first working version did. Three errors of mine that probing the claims caught before any test was written around them: - The drive advanced its phase by omega per step rather than by omega * dt. That simulates a frequency 1/S times too high while still producing a perfectly clean exponential, so it would have passed any test that only checked the profile was exponential. With the units right the measured decay is 0.20832 against a predicted 0.208360. - The decay fit ran the whole length of the guide. An evanescent field reaches a numerical floor within a few decay lengths, and the flat stretch beyond does not merely add scatter -- it drags the fitted slope towards zero, and on a long enough guide reports no decay at all. The fit now stops while the signal is still fifty times above the floor, and a correlation below 0.999 in the log is reported as a failed measurement rather than returned as a number. - The fit also began a full guide width downstream, to let higher modes die. They are never excited: the source is the mode's own transverse pattern and the discrete sine vectors are exactly orthogonal. Backing off that far merely threw away the dynamic range a fast-decaying high mode needs, and starting three cells out is what let modes two and three work at all. What comes back is the *numerical* cutoff, (2/S) arcsin(S sin(ky/2)), not the textbook m pi / a. The grid has its own dispersion relation and is always the slower of the two; the shortfall is (ky/2)^2 (1-S^2)/6, which the tests check directly. Reporting the continuum figure would be reporting what the answer ought to be rather than what the simulation has. Measured against the numerical value the agreement is a couple of parts in a thousand across widths 12 to 20 and modes one to three. A mode near the grid's resolution limit -- three half waves across ten cells -- decays within a couple of cells and leaves too little profile above the floor to fit. That returns NoConvergence, which is the honest answer. 5 unit tests and 5 property tests added: the matched layer returning a thousandth of what a conductor does and improving with depth, exact mirror and transpose symmetry of the Yee grid bit for bit, linearity, the numerical cutoff zeroing the decay to rounding and closing on the continuum value as the square of the cell size, and the measured decay recovering the grid's own cutoff more closely than the continuum one. Suite is 4,119 lib + 523 property tests, green in debug, clippy clean under --all-targets -D warnings, checked on nightly-2025-11-21. CI confirmed green on all five jobs for 7fff238 before this push. Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
chebyshev_points and cheb_diff_matrix, chebyshev_collocation_bvp for -(p u')' + q u = f with Dirichlet ends, spectral_poisson_periodic and its inverse spectral_second_derivative, and spectral_convergence_demo. The roadmap asked for the periodic solver to wrap the existing fft_poisson_2d. It does not, and the doc says why: that function divides by the eigenvalue of the *five-point* Laplacian, which makes the discrete residual vanish to rounding -- exactly what a pressure projection wants, since there the finite-difference divergence is the thing that must be zero -- but leaves it second-order accurate against the continuum. spectral_poisson_periodic divides by the true symbol -k^2 instead. The two solve different problems and both are right, and a test demonstrates the split directly: the spectral solution satisfies the continuum equation to 1e-12 and the three-point difference equation only to O(h^2), with the difference residual falling by exactly four per refinement. The differentiation matrix takes its diagonal as minus the sum of the rest of its row rather than from the closed form. The two agree analytically and differ in floating point by cancellation that grows with n; the sum makes the matrix annihilate constants by construction, which matters because the constant is the one thing every derivative operator must kill and an error in it pollutes everything else. The convergence claim is stated as what it actually is. "Spectral methods converge exponentially" is false as a property of the method and true as a property of smooth data, so the tests ask which model the errors follow rather than how small they are: for analytic data log(error) is linear in n, for data with k continuous derivatives it is linear in log(n), and comparing the two correlations separates the cases without naming a rate. |x|^3 comes out at n^-2.2 and |x|^5 at n^-4.7, with the smoother one more accurate at every size. One test needed its window chosen per function rather than fixed. Geometric convergence runs into the rounding floor, and past that the recorded errors are cancellation noise; cos(2x) is entire and is there by n = 16 while 1/(2+x) has a pole a unit from the interval and is still converging at n = 24. Fitting a model across the floor is fitting nothing, and that is what the first version of the test did -- it reported that an entire function follows a power law. Cross-validation is against a different discretisation rather than against itself: twenty-four collocation points and four hundred linear elements solve the same variable-coefficient problem and agree to the accuracy of the weaker one. 10 unit tests and 7 property tests: exactness on every polynomial the space holds, D applied twice giving the second derivative, centro-antisymmetry and zero row sums, the Jacobian being the only thing an interval change introduces, the periodic solver exact within the band and mean-free with a nonzero source mean dropped, linearity, the collocation patch test, and the finite-element cross-check. Section 19c is now complete: fem1d, fem2d, fdtd, spectral_pde. Suite is 4,129 lib + 530 property tests, green in debug, clippy clean under --all-targets -D warnings, checked on nightly-2025-11-21. CI confirmed green on all five jobs for 062a3a9 before this push. Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Roadmap section 19d, first part. New learn/ module. nn.rs holds Mlp with Act and Loss, forward, backward, numerical_grad_check, train_sgd and train_adam, predict, preactivations, conv2d_forward, and linear_regression_gd_check. The module is built around the observation that a learning algorithm is unusually easy to test badly. A falling training curve is not evidence of anything: gradient descent reduces the loss under a wrong gradient too, just more slowly and towards somewhere else. So no test here uses a loss curve as its main assertion. What settles backpropagation is the central difference, and that is asserted across random architectures, activations, losses and inputs, to eight digits. Softmax and cross-entropy are fused rather than composed. Taken separately the activation has a Jacobian and the loss has a gradient; taken together the product collapses to exactly p - y at the logits. That cancellation is worth having for accuracy as much as speed, since computing the two separately loses precision precisely where the network is confident. Cross-entropy therefore requires a softmax output, and a softmax output with squared error is refused rather than silently computing something else -- it would need the full Jacobian, which is not implemented. One test was rewritten rather than tuned. The rectifier gradient check disagrees with a central difference whenever a pre-activation lands within the difference step of zero, because the derivative genuinely does not exist at the kink. My first version asserted a pass rate, and at 90 of 120 random architectures that rate was neither meeting the threshold nor meaning anything -- it depends on the widths and depths drawn. The test now asserts the *cause*: every disagreement is required to have a pre-activation within a thousand difference steps of zero, which a genuinely wrong gradient would fail while sitting nowhere near one. Making that checkable is why preactivations is public. Other exact properties asserted: softmax invariant under a shift of its input, including at magnitudes where the naive computation overflows and where the order of the outputs must still match the order of the logits; a bias-free rectifier network positively homogeneous at any depth; permuting a hidden layer's units together with the next layer's columns leaving the computed function untouched, which is why two networks cannot be compared parameter by parameter; convolution linear and, away from the padding it cannot be shift invariant in, commuting with a shift; a uniform kernel giving the window mean and a delta kernel the identity, both exactly. XOR gets its classical treatment: a single layer cannot get below the 0.125 that predicting the mean costs, and one hidden layer solves it. And descent on linear least squares is checked against the closed form through the normal equations rather than against itself, with the step size taken as the reciprocal of the largest eigenvalue of X^T X, which is the largest step for which descent on a quadratic converges at all. 11 unit tests and 9 property tests. Suite is 4,140 lib + 539 property tests, green in debug, clippy clean under --all-targets -D warnings, checked on nightly-2025-11-21. CI confirmed green on all five jobs for 9d6e269 before this push. Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Roadmap section 19d, second part. gp.rs holds KernelFn (Rbf, Matern32, Matern52, Periodic, Linear, and closure under Sum and Product), Gp with fit by Cholesky, predict returning mean and variance, log_marginal_likelihood, optimize_hyperparams by Nelder-Mead over the logarithms, condition_estimate, sample_prior and sample_posterior. Regression here is conditioning rather than fitting: there is no optimisation in fit, only a factorisation, and the answer is exact given the kernel. Two consequences are asserted with `==` rather than a tolerance, because they are identities: - The posterior variance does not depend on the observations at all. Two processes fitted to the same inputs with entirely different targets return variances agreeing bit for bit. Uncertainty in a Gaussian process is a statement about where the data is, not about what it said. - A periodic kernel repeats exactly. The separation enters through a sine of half the lag over the period, so k(x, x+p) equals k(x,x) to the last bit rather than to a tolerance. Softmax-style fusion has an analogue here: log|K| is read off the Cholesky diagonal rather than formed as a determinant, which for any sizeable n underflows. The value is cross-checked in the tests against an independent LU determinant and solve that share no code with it. Hyperparameters are optimised in log space, which keeps every one positive without a constraint and makes the search scale-free -- a length scale of 0.01 and one of 100 are the same distance from 1, which is how they should be treated when nothing is known about the scale. condition_estimate is public because of what the tests found. Fitting a squared exponential to points spaced well inside its length scale gives a covariance matrix that is singular to working precision, and the jitter that makes the Cholesky succeed is then what limits the interpolation accuracy: the error is the jitter times the condition number, which reached 1e-4 at a condition number of 2e5 in the property tests. My first version asserted a fixed 1e-8 and was measuring the conditioning rather than the method. The tolerance is now that product, and the doc says the remedy is a shorter length scale, a rougher kernel or a nonzero noise -- statements about the model, not about the arithmetic. The estimate is documented as a lower bound, since the factor's diagonal says nothing about how the off-diagonal mass is arranged, which is why the tolerance carries a safety factor over it. One other test was rewritten rather than tuned: I had asserted that unit noise leaves a residual above 0.5, which depends on the data's amplitude and spacing and is not a property of anything. It now asserts the limit that is -- overwhelming noise collapses the posterior mean onto the prior's. 10 unit tests and 8 property tests: exact interpolation and vanishing variance at noiseless data, the variance ignoring the targets while the mean is exactly linear in them, conditioning never raising uncertainty anywhere, the prior recovered far from data, every kernel giving a positive semi-definite Gram matrix including the compound ones, the marginal likelihood against an independent determinant, tuning never lowering it, and both samplers reproducing the distributions they came from within their own standard errors. Suite is 4,150 lib + 547 property tests, green in debug, clippy clean under --all-targets -D warnings, checked on nightly-2025-11-21. CI confirmed green on all five jobs for 32def5b before this push. Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
…t trap Roadmap section 19d, third part. cluster.rs holds kmeans and kmeans_once with kmeans_pp_init and elbow_data, dbscan, hierarchical_agglomerative with four linkages and dendrogram_cut, gaussian_mixture_em, silhouette_score, adjusted_rand_index, davies_bouldin, knn_classify and knn_regress. THE RNG FINDING, which is the important part of this commit. A property test asserting that two independent random partitions score near zero on the adjusted Rand index failed with a score of exactly one. The cause is in monte_carlo: Rng is a plain linear congruential generator that returns its raw state, and for such a generator bit k has period at most 2^(k+1). Taking `next_u64() % m` for a power of two m reads exactly those bits. Measured: `% 2` gives 0,1,0,1 for ever, `% 4` gives 0,3,2,1, `% 8` has period eight. Two "random" label sequences drawn one after another are therefore perfectly correlated -- not a subtle statistical weakness but no randomness at all. A modulus with an odd factor mixes in higher bits and is fine, which is why this went unnoticed. Thirty-one sites across eleven files were drawing small integers this way, most of them in tests written over many earlier sessions. They were not producing wrong answers -- the invariants they assert hold for any input -- but their coverage was a repeating cycle of length two, four or eight rather than the random spread the code reads as. Fixed by adding Rng::below, which takes its answer from the top of the word, documenting the hazard on next_u64, and converting all thirty-one sites. The full suite passes with the widened randomisation, so no latent defect was hiding behind the narrow coverage -- but that was worth finding out rather than assuming. Two defects in this session's own code, both found the same way: - gaussian_mixture_em initialised every covariance at the *global* spread of the data. A component wide enough to cover the whole dataset claims every point almost equally, so the first maximisation dragged all the means back to the global mean and threw away the k-means initialisation that had just been computed. It converged to a local optimum eighty log-units worse than the right one. Seeding each covariance from its own cluster's scatter fixes it. - kmeans ran Lloyd's algorithm once. On three well-separated blobs about one run in two hundred lands on a stable configuration with two centres inside one blob and one spanning the other two -- inertia 800 against the best 21, and no number of iterations escapes it because no single point wants to move. It now restarts ten times and keeps the lowest inertia, with kmeans_once left public so a single monotone trajectory can still be observed. One assertion was wrong rather than the code: I had claimed the restarted run always beats a single one. Best-of-ten is a minimum over its own draws and says nothing about an independent eleventh, so the test now compares the two distributions, which is the claim restarts actually support. Centroid linkage is documented and tested as *inverting* rather than quietly producing dendrograms that cannot be drawn: merging two clusters puts their centre between them, which can be nearer a third than either original was. The property test requires an inversion to actually occur across random point sets, so the caveat cannot rot. DBSCAN's core points are asserted invariant to input order while border points are explicitly not, which is the algorithm as defined. 11 unit tests and 11 property tests. Suite is 4,161 lib + 558 property tests, green in debug, clippy clean under --all-targets -D warnings, checked on nightly-2025-11-21. CI confirmed green on all five jobs for e67b208 before this push. Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Roadmap section 19d, fourth part, completing learn/. tree.rs holds gini and entropy, decision_tree_fit and regression_tree_fit sharing one split search, tree_predict and tree_predict_value, feature_importance, random_forest_fit with bagging and per-split feature subsampling, and gradient_boosting_lite under squared loss. The tests are built around what a tree can *express* rather than how often it happens to be right, because the second depends on the data and the first does not: - A tree of depth d has at most 2^d leaves and can therefore name at most 2^d distinct classes, whatever it is given. My first attempt asserted an accuracy bound instead -- that one split cannot get more than half of four quadrants right -- and it was wrong: unequal quadrant counts let a leaf's majority exceed a quarter, and the stump reached fifty-five per cent. The bound on leaves holds always. - Splits are decided by the order of a column's values, not their magnitudes, so any increasing affine rescaling of any feature leaves the tree computing the same function, node for node. That is asserted across random columns and factors spanning six orders of magnitude, and it is what distinguishes trees from every distance-based method in this crate -- k-means, k-nearest-neighbours and a Gaussian process all give different answers under the same rescaling. - Feature importances are nonnegative and sum to exactly the total weighted impurity the tree removed, so dividing the credit among the columns neither creates nor loses any. - A regression tree's every leaf is a mean of training targets, so no prediction can leave their range -- checked a million units outside the training data. That is the same statement as "a tree never extrapolates", which is what makes it safe against runaway outputs and useless for trends. - Boosting's loss falls at every round, and a tree of depth zero cannot split, so a round of them must change nothing at all. The recorded first loss is the variance of the targets exactly, since the model starts at their mean, and gbm_predict is checked against the fit's own bookkeeping rather than trusted. Gini and entropy are asserted at their exact values -- zero for a pure node, exactly 1 - 1/k and exactly ln k for k equal classes -- and shown blind to the order of the counts and to scaling them all together, since they see proportions. The module documentation says why forests and boosting are opposite strategies rather than variants: a forest averages deep overfitted trees whose errors are decorrelated, while boosting adds shallow underfitted ones each fitted to the previous residual. Which is also why a forest's round count is harmless and boosting's has to be stopped early. 7 unit tests and 7 property tests. Suite is 4,168 lib + 565 property tests, green in debug, clippy clean under --all-targets -D warnings, checked on nightly-2025-11-21. Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Roadmap section 19d, final part, completing Part 4. src/units.rs became src/units/mod.rs unchanged, with two submodules alongside it. quantity.rs: Dim as seven i8 SI exponents with exact mul/div/pow/sqrt, Quantity with dimension-checked arithmetic and about thirty constructors, parse_unit, parse_quantity, unit_convert, si_prefixes_format, and the 2022 CODATA table. dimensional.rs: buckingham_pi, is_dimensionless_group, dimensionless_groups_named, natural_units_power and natural_units_convert, planck_units. Buckingham's theorem is a rank computation, so it is done over the crate's exact Rational rather than in floating point. A group is *exactly* in the null space or it is not, and one whose dimensions cancel to 1e-16 rather than to zero is a rounding error about to be reported as physics. The returned exponents stay rational for the same reason: Reynolds happens to have integer exponents, a general null space basis does not, and rounding it would silently change the group. is_dimensionless_group compares each row against zero, not against a tolerance. Two parser decisions the tests forced, both now documented rather than implicit: - Juxtaposition is multiplication. The CODATA table's own units are written "J s" and "1/mol", and the first version rejected both -- which was caught by a test that parses every unit in the table rather than trusting it. Whitespace now separates factors and a bare "1" is a valid term. - Parentheses are not supported, so "J/(mol K)" is refused as an unknown unit rather than quietly parsed as something else. A "/" applies to the single term after it, so the table writes "J/mol/K". Refusing is the safer of the two ways not to support them. The gram carries the SI prefixes rather than the kilogram, so kg comes out at exactly one and mg at 1e-6 -- the kilogram being the only base unit whose name already contains a prefix. And unit names are resolved whole before any prefix is split off, which is what makes m a metre, mm a millimetre, min a minute and T a tesla. That is a rule rather than a deduction and the doc says so, because any other rule gives different answers for the same strings. The CODATA table is checked for internal consistency rather than transcribed and trusted: the seven constants that are exact by definition since the 2019 SI revision are asserted exactly, the gas constant is the product of two of them, epsilon_0 mu_0 c^2 is one, and the fine-structure and Rydberg constants are recomputed from the others. The Planck units are derived from hbar, c and G rather than copied, and the tests check their defining relations -- l_P = c t_P, E_P = m_P c^2, and the Schwarzschild radius of the Planck mass being twice the Planck length, which is the statement that gravity and quantum mechanics meet there. Natural units refuse a dimension involving amperes, kelvin, moles or candela rather than guessing at a convention to absorb them, and the property test checks the bookkeeping is consistent by requiring the converted magnitudes to multiply. 6 unit tests in quantity, 5 in dimensional, and 8 property tests. Suite is 4,180 lib + 573 property tests, green in debug, clippy clean under --all-targets -D warnings, checked on nightly-2025-11-21. CI confirmed green on all five jobs for 566d8d6 before this push. Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
…able
Completes the last two items on the Part 4 roadmap: the exact/symbolic
-> units/ cross-reference, and the constant consolidation.
dimensional_check_formula walks an Expr and returns its dimension,
enforcing the two rules a hand derivation drops: every term of a sum
has to have the same dimension, and a transcendental's argument has to
be dimensionless. Neither can be checked by evaluating the formula --
both sides of `x + v` are perfectly good floats -- so this is a check
numerical testing cannot do.
Two rules the tests forced, both of which came out of running the
checker on real `diff` output rather than out of theory:
* Zero is the additive identity of every dimension at once, so a
literally-zero term joins any sum. `diff` does not simplify, so the
product rule leaves `0 * t` sitting beside `v * 1`, and a checker
that refused that sum would be useless on anything differentiated.
The waiver does not reach inside the zero term: its subexpressions
are still checked.
* A Const exponent is read as the dyadic rational it exactly is, via
Rational::from_f64_exact. `0.5` is one half, so Pow(x, 0.5) is a
square root -- which matters because that is exactly how `diff`
writes the derivative of one. `0.1` is not one tenth but the
power-of-two fraction the float holds, and no dimension is
divisible by that denominator, so `l^0.1` is reported as a root
that does not exist rather than rounded into one that does.
Constants: math::constants was already the one table, but four modules
carried their own copies, two of them at different values.
* chemistry::FARADAY was 96485.0, which differs from N_A e in the
sixth digit. math::constants now computes FARADAY from its two
exact factors and chemistry re-exports it.
* particle_physics::FINE_STRUCTURE was 7.297e-3 against ALPHA's
7.2973525693e-3, a disagreement at 5e-5 relative.
* habitable_zone and magnetosphere each had their own SOLAR_TEMPERATURE
and SOLAR_LUMINOSITY; atmosphere had its own STANDARD_PRESSURE.
kinetics_props built the Nernst slope from its own transcribed 96485.0
and 8.314462618, at 1e-9 tolerance, so it was asserting the
transcription rather than the formula. It now uses the crate's R and
FARADAY, and caught the Faraday change -- which is what it should have
been doing all along.
Tests. The two constant tables are pinned to each other: the eight
constants fixed by the 2019 SI redefinition must agree bit-for-bit,
since nothing but a transcription error could move them, and the eleven
measured ones to 1e-8, which separates the 2018-to-2022 CODATA revision
(1e-13 to 1.5e-9 observed) from a mistyped digit. Derived constants are
checked against their definitions rather than against themselves, and
every unit string in the CODATA table has to parse.
For the checker, the cross-checks matter more than the direct ones: it
agrees with Quantity arithmetic, which implements the same algebra with
no code in common; it certifies the groups buckingham_pi finds by an
exact null space over the rationals as dimensionless; and the
derivative of anything has the dimension of the thing over the
variable's, which makes it a check on the differentiator too. The
property generator builds an expression alongside the dimension its
construction guarantees, so agreeing is evidence rather than tautology.
4193 lib tests, 577 property tests, clippy clean, nightly-2025-11-21
clean.
Co-Authored-By: Claude <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
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
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.
Completes Phase A and opens Phase B of
docs/ROADMAP_PART4.md.exact/symbolicdiscrete/primesdiscrete/number_theoryA rare branch that no test reached
After the primes tests passed, I replaced the entire strong Lucas test — the second half of BPSW — with
return true. The whole suite stayed green.is_prime_bigintshort-circuits below 2⁶² and otherwise runs Miller-Rabin on random bases first, so essentially every composite is rejected before Lucas is reached.It is now tested directly: over the odd numbers below 20,000 it must accept every prime and exactly the strong Lucas pseudoprimes for Selfridge's parameters. The implementation reproduces that set — 5459, 5777, 10877, 16109, 18971 — computed independently rather than read off the code. The test also pins the property BPSW rests on: none of the six base-2 strong pseudoprimes below 20,000 (2047, 3277, 4033, 4681, 8321, 15841) is a Lucas pseudoprime, so each test catches what the other misses. The mutation now fails.
This is the third rare-branch gap found this way, after Knuth's add-back correction and the Lentz convergence threshold.
Other findings
Two missing simplifier rules.
split_coeffdidn't recognise negation, so its −1 stayed an opaque factor andx − xhashed under two keys instead of cancelling. And products of exponentials didn't combine, soexp(x)·exp(−x)couldn't reach 1.Dualwas missingatan,sinh,cosh. The roadmap asks that symbolic derivatives be cross-checked against forward-mode automatic differentiation; those functions made it impossible. Backfilled withtanh, tested against closed forms andcosh² − sinh² = 1.Notes
π(10⁹) = 50,847,534 in 87 ms in a debug build, via the Lucy_Hedgehog recurrence over the distinct values of n/i — O(√n) state, so the roadmap's headline value is an ordinary test rather than something skipped for cost.
simplify(diff(sin²+cos²)) == 0falls out of ordinary term collection, not a hard-coded trig rewrite: the product rule yields+2·cos·sinand−2·cos·sin, and they collect.Three specification readings are documented in the commits:
quadratic_diophantine_solveas a definite form,frobenius_numberwith a unit coin, andstern_brocot_nthindexed breadth-first.Verification
2,981 lib tests, 107 property tests,
cargo clippy --all-targets -- -D warningsclean. Each commit was verified by extracting its staged tree into a clean checkout and confirming the committed tree hash matched what was tested there.🤖 Generated with Claude Code
https://claude.ai/code/session_01VUkMGnfUbaopkUJLWeYCYi
Generated by Claude Code