{"messages":[{"content":"You are a constitutional council ranking individual git commits for ownership allocation.\n\nCompare these two commits. Decide which contributed more lasting value to the project.\n\nJudge substance, not spectacle:\n- Prefer correct, lasting design and real bugfixes over churn, formatting, renames, or generated noise.\n- Prefer clarity and necessity over sheer line count. A small precise change can beat a large diffuse one.\n- Do not favor a side merely because its patch is longer or noisier.\n- Weight what the change does for the project, not the contributor's name.\n\nReturn ONLY a JSON object: {\"winner\": \"A\" or \"B\", \"ratio\": \"N:M\", \"explanation\": \"...\"}\nThe explanation must cite concrete differences in the patches (1-3 sentences).\n\nSide A — contributor: tommy-mor\nSide A — commit message:\n[07cc0e1c] Bind HTTP before replaying the ledger so boots stay reachable.\n\nImport-time JsonlStore replay of a 450MB evidence log blocked uvicorn from\nlistening; defer load to a startup task and keep /api/health always live.\n\nCo-authored-by: Cursor <cursoragent@cursor.com>\n\nSide A — unified diff (full patch):\ndiff --git a/constitution.py b/constitution.py\nindex eb2779c3c6775ed26bb932320635ffdfc6edd439..17e322a8884a09d7579676ae221b705f46039f6e 100644\n--- a/constitution.py\n+++ b/constitution.py\n@@ -27,8 +27,9 @@ Run: uv run constitution.py\n from decimal import Decimal, getcontext, DefaultContext\n from datetime import datetime, timezone\n from fastapi import FastAPI, Request, Response\n-from fastapi.responses import PlainTextResponse, HTMLResponse\n+from fastapi.responses import PlainTextResponse, HTMLResponse, JSONResponse\n from starlette.middleware.sessions import SessionMiddleware\n+from starlette.requests import Request\n import json, time, os, asyncio, httpx, pathlib, subprocess, hashlib, re, fcntl, base64\n import sympy as sp  # type: ignore[reportMissingImports]\n from tenacity import retry, retry_if_exception, stop_after_attempt, wait_exponential\n@@ -269,7 +270,42 @@ class Evidence:\n     payload: dict\n \n \n-store = JsonlStore(JSONL_PATH)\n+class DeferredJsonlStore:\n+    \"\"\"\n+    Delay ledger replay until after the HTTP server binds.\n+\n+    A 450MB+ evidence ledger can take minutes to parse; doing that at import\n+    time makes Fly health checks and the public site time out on every boot.\n+    \"\"\"\n+\n+    def __init__(self, path: pathlib.Path | str):\n+        self.path = pathlib.Path(path)\n+        self._inner: JsonlStore | None = None\n+        self.ready = False\n+\n+    def load_sync(self) -> None:\n+        self._inner = JsonlStore(self.path)\n+        self.ready = True\n+\n+    def _require(self) -> JsonlStore:\n+        if self._inner is None:\n+            self.load_sync()\n+        assert self._inner is not None\n+        return self._inner\n+\n+    def read(self) -> list:\n+        if self._inner is None:\n+            return []\n+        return self._inner.read()\n+\n+    async def append(self, e) -> None:\n+        await self._require().append(e)\n+\n+    async def atomic(self, fn):\n+        return await self._require().atomic(fn)\n+\n+\n+store = DeferredJsonlStore(JSONL_PATH)\n _LEDGER_LOCK = asyncio.Lock()\n \n \n@@ -2211,10 +2247,29 @@ async def get_ledger(offset: int = 0, limit: int = 100, full: int = 0):\n     return rows\n \n \n+@app.middleware(\"http\")\n+async def ledger_ready_gate(request: Request, call_next):\n+    if request.url.path == \"/api/health\":\n+        return await call_next(request)\n+    if isinstance(store, DeferredJsonlStore) and not store.ready:\n+        return JSONResponse(\n+            {\n+                \"ok\": False,\n+                \"loading\": True,\n+                \"message\": \"replaying ledger into memory\",\n+            },\n+            status_code=503,\n+        )\n+    return await call_next(request)\n+\n+\n @app.get(\"/api/health\")\n async def get_health():\n-    \"\"\"Liveness only — must not touch the ledger (health checks during ranking).\"\"\"\n-    return {\"ok\": True}\n+    \"\"\"Liveness only — must not touch the ledger (boot + ranking).\"\"\"\n+    return {\n+        \"ok\": True,\n+        \"ledger_ready\": bool(getattr(store, \"ready\", True)),\n+    }\n \n \n @app.get(\"/api/epoch\")\n@@ -3703,8 +3758,17 @@ async def callback(request: Request, code: str):\n \n @app.on_event(\"startup\")\n async def startup():\n-    if os.environ.get(\"DISABLE_EPOCH_LOOP\") != \"1\":\n-        asyncio.create_task(epoch_loop())\n+    async def boot():\n+        if isinstance(store, DeferredJsonlStore) and not store.ready:\n+            print(\"replaying ledger…\", flush=True)\n+            await asyncio.to_thread(store.load_sync)\n+            n = len(store.read())\n+            print(f\"ledger ready: {n} events\", flush=True)\n+        if os.environ.get(\"DISABLE_EPOCH_LOOP\") != \"1\":\n+            asyncio.create_task(epoch_loop())\n+\n+    # Do not block bind/health on ledger replay.\n+    asyncio.create_task(boot())\n \n \n if __name__ == \"__main__\":\n\n\nSide B — contributor: tommy-mor\nSide B — commit message:\n[0e81b5e1] Solve Rank Centrality exactly instead of capping power iteration.\n\n`compute_scores_from_edges` ran power iteration for at most 10 000 sweeps and\nreturned whatever vector it had reached, with no signal when the cap bound. On\na preference chain the spectral gap is Theta(1/n^2), so the cap binds past\nn ~= 1500 -- and worse, the L1 *step* tolerance is an absolute test, so on a\n1024-node chain with mixed ratios it declared success after 4325 sweeps while\nleaving 486 of 1024 items at the wrong rank. A 1024-node star fails outright:\nit needs 14 154 sweeps.\n\nNew `stationary` module owns every candidate solver -- power, Aitken-accelerated\npower, Gauss-Seidel/SOR, preconditioned BiCGSTAB, dense LU, dense GTH, and a\nsparse GTH state reduction with minimum-degree elimination -- behind one `solve`\nthat splits disconnected components, tries the sparse direct path first (trees\nand chains reduce in O(n) with zero fill), and falls back to Gauss-Seidel then\npower iteration on the dense graphs elimination cannot help with. Every answer\ncarries its own recomputed residual and a `converged` flag; nothing comes back\nclaiming success it did not earn.\n\nScores are also produced in log space. A chain of n items each preferred 2:1\nspans 2^(n-1), which leaves the f64 range at n ~= 1075, so `pi` genuinely cannot\norder a long chain no matter how it is computed. `Solution::log_pi` stays exact\nthere and is what `ranked_items_subset` now sorts on.\n\nEdges are sorted before aggregation, so the result no longer depends on\n`HashMap` iteration order. The legacy builder produced a different bit pattern\non 7 of 7 rebuilds of the same graph; the sort costs ~10% of chain construction.\n\nAdds `benches/stationary_solvers.rs` (~80s) and `examples/solver_probe.rs`\n(~30s), which reports time, residual, and rank displacement against closed-form\nanswers per topology. No new dependencies.\n\nCo-authored-by: Cursor <cursoragent@cursor.com>\n\nSide B — unified diff (full patch):\ndiff --git a/server/Cargo.toml b/server/Cargo.toml\nindex 7de0c944b4a8c7e3123b818432fce82f09e6aeca..cdc2a78be0ffbd6ac27a02a332093daf5eb8b990 100644\n--- a/server/Cargo.toml\n+++ b/server/Cargo.toml\n@@ -49,4 +49,8 @@ harness = false\n name = \"ingest_replay\"\n harness = false\n \n+[[bench]]\n+name = \"stationary_solvers\"\n+harness = false\n+\n \ndiff --git a/server/benches/stationary_solvers.rs b/server/benches/stationary_solvers.rs\nnew file mode 100644\nindex 0000000000000000000000000000000000000000..ff442b86e61861633ce2be1a016eae53aec73e6e\n--- /dev/null\n+++ b/server/benches/stationary_solvers.rs\n@@ -0,0 +1,171 @@\n+//! Head-to-head cost of the stationary-distribution solvers behind Rank Centrality.\n+//!\n+//! Deliberately small: a couple of representative sizes per topology, ten\n+//! samples each, ~1 s of measurement per point. The point is the *ordering*\n+//! between solvers, which is stable at this resolution; for a finer sweep plus\n+//! accuracy numbers use `cargo run --release -p slugsocial-server --example\n+//! solver_probe`, which finishes in seconds.\n+//!\n+//! What each topology is here to show:\n+//!\n+//! - **chain** — spectral gap `Θ(1/n²)`. Power iteration needs `Θ(n²)` sweeps\n+//!   and blows the 10 000 cap past n ≈ 1500; sparse GTH is `O(n)` because a\n+//!   tree eliminates with zero fill.\n+//! - **star** — uniformization by `d_max = n-1` makes the chain almost purely\n+//!   lazy, so power iteration crawls (it fails the cap at n = 1024) while\n+//!   Gauss–Seidel finishes in two sweeps.\n+//! - **clique / sparse** — well-conditioned. Iterative wins outright and the\n+//!   direct methods are the ones paying.\n+\n+use std::hint::black_box;\n+use std::time::Duration;\n+\n+use criterion::{criterion_group, criterion_main, BenchmarkId, Criterion};\n+\n+use slugsocial_server::ranking::chain_from_edges;\n+use slugsocial_server::stationary::{\n+    bicgstab, dense_gth, dense_lu, power, power_aitken, solve, sor, sparse_gth, RankChain,\n+    SolveOptions, SparseGthOutcome,\n+};\n+\n+mod common;\n+use common::{build_group, Topology};\n+\n+const TOL: f64 = 1e-8;\n+const CAP: usize = 10_000;\n+\n+fn opts() -> SolveOptions {\n+    SolveOptions {\n+        tol: TOL,\n+        max_iters: CAP,\n+        ..SolveOptions::default()\n+    }\n+}\n+\n+fn chain_of(topo: Topology, n: usize) -> RankChain {\n+    let g = build_group(topo, n);\n+    chain_from_edges(g.idx_to_item.len(), g.edges.iter().map(|(&k, &w)| (k, w)))\n+}\n+\n+fn quick<'a>(\n+    c: &'a mut Criterion,\n+    name: &str,\n+) -> criterion::BenchmarkGroup<'a, criterion::measurement::WallTime> {\n+    let mut g = c.benchmark_group(name);\n+    g.sample_size(10)\n+        .warm_up_time(Duration::from_millis(300))\n+        .measurement_time(Duration::from_secs(1));\n+    g\n+}\n+\n+/// Every solver on the same input, per topology.\n+fn solver_shootout(c: &mut Criterion) {\n+    let cases: Vec<(Topology, usize, bool)> = vec![\n+        // (topology, n, include the O(n³) dense arms)\n+        (Topology::Chain, 256, true),\n+        (Topology::Chain, 1024, true),\n+        (Topology::Star, 1024, true),\n+        (Topology::Clique, 256, true),\n+        (Topology::RandomSparse { degree: 6 }, 1024, false),\n+    ];\n+\n+    for (topo, n, dense) in cases {\n+        let chain = chain_of(topo, n);\n+        let mut g = quick(c, &format!(\"solvers/{}\", topo.label()));\n+\n+        g.bench_with_input(BenchmarkId::new(\"power\", n), &chain, |b, ch| {\n+            b.iter(|| black_box(power(ch, opts()).residual))\n+        });\n+        g.bench_with_input(BenchmarkId::new(\"power+aitken\", n), &chain, |b, ch| {\n+            b.iter(|| black_box(power_aitken(ch, opts()).residual))\n+        });\n+        g.bench_with_input(BenchmarkId::new(\"sor\", n), &chain, |b, ch| {\n+            b.iter(|| black_box(sor(ch, opts(), 1.0).residual))\n+        });\n+        g.bench_with_input(BenchmarkId::new(\"bicgstab\", n), &chain, |b, ch| {\n+            b.iter(|| black_box(bicgstab(ch, opts()).residual))\n+        });\n+        g.bench_with_input(BenchmarkId::new(\"sparse-gth\", n), &chain, |b, ch| {\n+            let unlimited = SolveOptions {\n+                direct_work_budget: u64::MAX,\n+                dense_core_max: usize::MAX,\n+                ..opts()\n+            };\n+            b.iter(|| match sparse_gth(ch, unlimited, TOL) {\n+                SparseGthOutcome::Solved(s) => black_box(s.residual),\n+                SparseGthOutcome::TooDense { .. } => unreachable!(),\n+            })\n+        });\n+        if dense {\n+            g.bench_with_input(BenchmarkId::new(\"dense-gth\", n), &chain, |b, ch| {\n+                b.iter(|| black_box(dense_gth(ch, TOL).residual))\n+            });\n+            g.bench_with_input(BenchmarkId::new(\"dense-lu\", n), &chain, |b, ch| {\n+                b.iter(|| black_box(dense_lu(ch, TOL).residual))\n+            });\n+        }\n+        g.finish();\n+    }\n+}\n+\n+/// The shipped path, including the cost of deciding which solver to use.\n+fn hybrid_dispatch(c: &mut Criterion) {\n+    let cases: Vec<(String, Topology, usize)> = vec![\n+        (\"chain\".into(), Topology::Chain, 4000),\n+        (\"star\".into(), Topology::Star, 4000),\n+        (\"clique\".into(), Topology::Clique, 400),\n+        (\n+            \"sparse-d6\".into(),\n+            Topology::RandomSparse { degree: 6 },\n+            4000,\n+        ),\n+        (\n+            \"components\".into(),\n+            Topology::ManyComponents {\n+                components: 64,\n+                size: 64,\n+            },\n+            0,\n+        ),\n+    ];\n+    let mut g = quick(c, \"hybrid\");\n+    for (label, topo, n) in cases {\n+        let chain = chain_of(topo, n);\n+        g.bench_with_input(BenchmarkId::new(label, chain.n), &chain, |b, ch| {\n+            b.iter(|| black_box(solve(ch, opts()).residual))\n+        });\n+    }\n+    g.finish();\n+}\n+\n+/// Building the chain (aggregate, pairwise-normalize, sort) versus solving it.\n+/// Sorting is what buys run-to-run determinism; this is where to see its price.\n+fn chain_construction(c: &mut Criterion) {\n+    let mut g = quick(c, \"chain_from_edges\");\n+    for (topo, n) in [\n+        (Topology::RandomSparse { degree: 6 }, 4000usize),\n+        (Topology::Clique, 400),\n+    ] {\n+        let group = build_group(topo, n);\n+        let total = group.idx_to_item.len();\n+        g.bench_with_input(\n+            BenchmarkId::new(topo.label(), group.edges.len()),\n+            &group,\n+            |b, grp| {\n+                b.iter(|| {\n+                    let ch = chain_from_edges(total, grp.edges.iter().map(|(&k, &w)| (k, w)));\n+                    black_box(ch.nnz())\n+                })\n+            },\n+        );\n+    }\n+    g.finish();\n+}\n+\n+criterion_group!(\n+    benches,\n+    solver_shootout,\n+    hybrid_dispatch,\n+    chain_construction\n+);\n+criterion_main!(benches);\ndiff --git a/server/examples/solver_probe.rs b/server/examples/solver_probe.rs\nnew file mode 100644\nindex 0000000000000000000000000000000000000000..e3933530ac9efab379dc10ec2a467eb30ad203f8\n--- /dev/null\n+++ b/server/examples/solver_probe.rs\n@@ -0,0 +1,738 @@\n+//! Head-to-head comparison of stationary-distribution solvers for Rank Centrality.\n+//!\n+//! Run: `cargo run --release -p slugsocial-server --example solver_probe`\n+//!\n+//! Reports, per (topology, solver): wall time, backward error `‖πP − π‖₁`, and\n+//! how badly the produced *order* disagrees with ground truth. On a chain the\n+//! ground truth order is known exactly by construction, and on any tree the\n+//! exact scores are known in closed form from detailed balance, so \"how wrong\n+//! is the current implementation\" is answerable without trusting any solver.\n+\n+use std::time::Instant;\n+\n+use slugsocial_server::path_types::ItemId;\n+use slugsocial_server::ranking::chain_from_edges;\n+use slugsocial_server::reducer::{GroupState, VoteData};\n+use slugsocial_server::stationary::{\n+    bicgstab, dense_gth, dense_lu, power, power_aitken, solve, sor, sparse_gth, Method, RankChain,\n+    Solution, SolveOptions, SparseGthOutcome,\n+};\n+\n+const TOL: f64 = 1e-8;\n+const CAP: usize = 10_000;\n+\n+// ---------------------------------------------------------------------------\n+// Graph generators (mirrors server/benches/common/mod.rs so numbers line up)\n+// ---------------------------------------------------------------------------\n+\n+struct Rng(u64);\n+impl Rng {\n+    fn new(s: u64) -> Self {\n+        Rng(s | 1)\n+    }\n+    fn next_u64(&mut self) -> u64 {\n+        let mut x = self.0;\n+        x ^= x >> 12;\n+        x ^= x << 25;\n+        x ^= x >> 27;\n+        self.0 = x;\n+        x.wrapping_mul(0x2545_F491_4F6C_DD1D)\n+    }\n+    fn below(&mut self, n: usize) -> usize {\n+        (self.next_u64() % n as u64) as usize\n+    }\n+}\n+\n+fn item_name(i: usize) -> String {\n+    format!(\"~/bench/i{i:06}\")\n+}\n+\n+fn vote(a: usize, b: usize, l: i32, r: i32, ts: i64) -> VoteData {\n+    VoteData {\n+        ts,\n+        a: ItemId::parse(&item_name(a)).unwrap(),\n+        b: ItemId::parse(&item_name(b)).unwrap(),\n+        ratio_left: l,\n+        ratio_right: r,\n+        body: \"synthetic\".into(),\n+        principal: \"bench\".into(),\n+        delegate: None,\n+        thread_tag: \"bench\".into(),\n+    }\n+}\n+\n+#[derive(Clone, Copy)]\n+enum Topo {\n+    /// `0 > 1 > 2 > …`, every vote at a fixed ratio. Ground-truth order known.\n+    Chain {\n+        ratio: i32,\n+    },\n+    /// Chain with per-edge varying ratios (what `benches/common` builds).\n+    ChainVaried,\n+    Star,\n+    Clique,\n+    Sparse {\n+        degree: usize,\n+    },\n+    Components {\n+        count: usize,\n+        size: usize,\n+    },\n+}\n+\n+fn pairs(t: Topo, n: usize) -> Vec<(usize, usize, i32)> {\n+    let mut out = Vec::new();\n+    match t {\n+        Topo::Chain { ratio } => {\n+            for i in 0..n.saturating_sub(1) {\n+                out.push((i, i + 1, ratio));\n+            }\n+        }\n+        Topo::ChainVaried => {\n+            for i in 0..n.saturating_sub(1) {\n+                out.push((i, i + 1, 2 + (i % 5) as i32));\n+            }\n+        }\n+        Topo::Star => {\n+            for i in 1..n {\n+                out.push((0, i, 2 + (i % 5) as i32));\n+            }\n+        }\n+        Topo::Clique => {\n+            let mut k = 0;\n+            for i in 0..n {\n+                for j in (i + 1)..n {\n+                    out.push((i, j, 2 + (k % 5)));\n+                    k += 1;\n+                }\n+            }\n+        }\n+        Topo::Sparse { degree } => {\n+            let mut rng = Rng::new(0xC0FFEE ^ n as u64);\n+            for i in 0..n.saturating_sub(1) {\n+                out.push((i, i + 1, 2));\n+            }\n+            for k in 0..(degree.saturating_sub(2) * n / 2) {\n+                let a = rng.below(n);\n+                let b = rng.below(n);\n+                if a != b {\n+                    out.push((a, b, 2 + (k % 5) as i32));\n+                }\n+            }\n+        }\n+        Topo::Components { count, size } => {\n+            for c in 0..count {\n+                let base = c * size;\n+                for i in 0..size.saturating_sub(1) {\n+                    out.push((base + i, base + i + 1, 2));\n+                }\n+            }\n+        }\n+    }\n+    out\n+}\n+\n+fn node_count(t: Topo, n: usize) -> usize {\n+    match t {\n+        Topo::Components { count, size } => count * size,\n+        _ => n,\n+    }\n+}\n+\n+fn build_chain(t: Topo, n: usize) -> RankChain {\n+    let total = node_count(t, n);\n+    let mut g = GroupState::new();\n+    for i in 0..total {\n+        g.ensure_item_pub(&item_name(i));\n+    }\n+    for (k, (i, j, ratio)) in pairs(t, total).into_iter().enumerate() {\n+        g.apply_vote(vote(i, j, ratio, 1, k as i64));\n+    }\n+    chain_from_edges(total, g.edges.iter().map(|(&k, &w)| (k, w)))\n+}\n+\n+// ---------------------------------------------------------------------------\n+// Accuracy metrics\n+// ---------------------------------------------------------------------------\n+\n+/// Kendall-tau distance between the order induced by `a` and by `b`, as a\n+/// fraction of all pairs. 0 = identical ranking, 0.5 = as good as random.\n+fn kendall_distance(a: &[f64], b: &[f64]) -> f64 {\n+    let n = a.len();\n+    if n < 2 {\n+        return 0.0;\n+    }\n+    // Sort indices by `a` descending, then count inversions with respect to `b`.\n+    let mut idx: Vec<usize> = (0..n).collect();\n+    idx.sort_by(|&x, &y| b[y].partial_cmp(&b[x]).unwrap_or(std::cmp::Ordering::Equal));\n+    let mut seq: Vec<f64> = idx.iter().map(|&i| a[i]).collect();\n+    let mut buf = seq.clone();\n+    let inv = count_inversions(&mut seq, &mut buf);\n+    let total = (n as f64) * (n as f64 - 1.0) / 2.0;\n+    inv as f64 / total\n+}\n+\n+/// Inversions = pairs out of descending order (strict; ties are concordant).\n+fn count_inversions(v: &mut [f64], buf: &mut [f64]) -> u64 {\n+    let n = v.len();\n+    if n < 2 {\n+        return 0;\n+    }\n+    let mid = n / 2;\n+    let (l, r) = v.split_at_mut(mid);\n+    let (bl, br) = buf.split_at_mut(mid);\n+    let mut inv = count_inversions(l, bl) + count_inversions(r, br);\n+    let (mut i, mut j, mut k) = (0usize, 0usize, 0usize);\n+    while i < l.len() && j < r.len() {\n+        if l[i] >= r[j] {\n+            buf[k] = l[i];\n+            i += 1;\n+        } else {\n+            // every remaining element of `l` is smaller than r[j]: inversion\n+            inv += (l.len() - i) as u64;\n+            buf[k] = r[j];\n+            j += 1;\n+        }\n+        k += 1;\n+    }\n+    while i < l.len() {\n+        buf[k] = l[i];\n+        i += 1;\n+        k += 1;\n+    }\n+    while j < r.len() {\n+        buf[k] = r[j];\n+        j += 1;\n+        k += 1;\n+    }\n+    v.copy_from_slice(&buf[..n]);\n+    inv\n+}\n+\n+/// The pre-change chain builder, kept verbatim so determinism can be measured\n+/// rather than argued about: adjacency rows and their weight sums are produced\n+/// by iterating a `HashMap`, whose order varies per instance.\n+fn legacy_chain(n: usize, edges: impl Iterator<Item = ((usize, usize), f64)>) -> RankChain {\n+    use std::collections::{HashMap, HashSet};\n+    let mut raw: HashMap<(usize, usize), f64> = HashMap::new();\n+    for ((src, dst), w) in edges {\n+        if src >= n || dst >= n || w <= 0.0 {\n+            continue;\n+        }\n+        *raw.entry((src, dst)).or_insert(0.0) += w;\n+    }\n+    let keys: Vec<(usize, usize)> = raw.keys().copied().collect();\n+    let mut normalized: HashMap<(usize, usize), f64> = HashMap::new();\n+    for (i, j) in keys {\n+        if normalized.contains_key(&(i, j)) {\n+            continue;\n+        }\n+        let w_ij = *raw.get(&(i, j)).unwrap_or(&0.0);\n+        let w_ji = *raw.get(&(j, i)).unwrap_or(&0.0);\n+        let total = w_ij + w_ji;\n+        if total <= 0.0 {\n+            continue;\n+        }\n+        normalized.insert((i, j), w_ij / total);\n+        if w_ji > 0.0 {\n+            normalized.insert((j, i), w_ji / total);\n+        }\n+    }\n+    let mut rows: Vec<Vec<(usize, f64)>> = vec![Vec::new(); n];\n+    let mut neighbors: Vec<HashSet<usize>> = vec![HashSet::new(); n];\n+    for ((src, dst), w) in &normalized {\n+        rows[*src].push((*dst, *w));\n+        neighbors[*src].insert(*dst);\n+        neighbors[*dst].insert(*src);\n+    }\n+    let row_sum: Vec<f64> = rows\n+        .iter()\n+        .map(|r| r.iter().map(|&(_, w)| w).sum())\n+        .collect();\n+    let d_max = neighbors.iter().map(|s| s.len()).max().unwrap_or(0) as f64;\n+    RankChain {\n+        n,\n+        rows,\n+        row_sum,\n+        d_max,\n+    }\n+}\n+\n+/// Pairs whose relative order differs between two score vectors.\n+fn order_flips(a: &[f64], b: &[f64]) -> usize {\n+    let n = a.len();\n+    let mut idx: Vec<usize> = (0..n).collect();\n+    idx.sort_by(|&x, &y| b[y].partial_cmp(&b[x]).unwrap_or(std::cmp::Ordering::Equal));\n+    let mut seq: Vec<f64> = idx.iter().map(|&i| a[i]).collect();\n+    let mut buf = seq.clone();\n+    count_inversions(&mut seq, &mut buf) as usize\n+}\n+\n+/// Adjacent positions where a chain's known-correct descending order breaks,\n+/// split into strict inversions (genuinely the wrong way round) and ties (the\n+/// two items are indistinguishable, so their displayed order is arbitrary).\n+fn chain_order_violations(pi: &[f64]) -> (usize, usize) {\n+    let mut inv = 0;\n+    let mut tie = 0;\n+    for i in 0..pi.len().saturating_sub(1) {\n+        if pi[i] < pi[i + 1] {\n+            inv += 1;\n+        } else if pi[i] == pi[i + 1] {\n+            tie += 1;\n+        }\n+    }\n+    (inv, tie)\n+}\n+\n+/// Exact π for a chain via detailed balance, in log space (no overflow).\n+///\n+/// A chain is a tree, so the Markov chain is reversible and\n+/// `π_i / π_{i+1} = a_{i+1,i} / a_{i,i+1}` holds exactly edge by edge.\n+fn chain_exact_log(t: Topo, n: usize) -> Vec<f64> {\n+    let mut logs = vec![0.0f64; n];\n+    for (i, (_, _, ratio)) in pairs(t, n).into_iter().enumerate() {\n+        // node i preferred over node i+1 at `ratio`:1, so π_i / π_{i+1} = ratio\n+        logs[i + 1] = logs[i] - (ratio as f64).ln();\n+    }\n+    logs\n+}\n+\n+/// Largest absolute error in `ln π_i` against an exact log-space reference,\n+/// restricted to entries the reference says are representable in f64.\n+fn max_log_error(pi: &[f64], exact_log: &[f64]) -> (f64, usize) {\n+    let ref_max = exact_log.iter().cloned().fold(f64::NEG_INFINITY, f64::max);\n+    let pi_max = pi.iter().cloned().fold(0.0f64, f64::max);\n+    if pi_max <= 0.0 {\n+        return (f64::INFINITY, pi.len());\n+    }\n+    let mut worst = 0.0f64;\n+    let mut lost = 0usize;\n+    for i in 0..pi.len() {\n+        let want = exact_log[i] - ref_max;\n+        if want < -700.0 {\n+            continue; // genuinely below f64 range; not the solver's fault\n+        }\n+        if pi[i] <= 0.0 {\n+            lost += 1;\n+            continue;\n+        }\n+        let got = (pi[i] / pi_max).ln();\n+        worst = worst.max((got - want).abs());\n+    }\n+    (worst, lost)\n+}\n+\n+// ---------------------------------------------------------------------------\n+// Runner\n+// ---------------------------------------------------------------------------\n+\n+struct Run {\n+    name: &'static str,\n+    ms: f64,\n+    sol: Solution,\n+}\n+\n+fn time_it(name: &'static str, mut f: impl FnMut() -> Solution) -> Run {\n+    // One warm-up, then the best of three: solver cost here is deterministic,\n+    // so the minimum is the cleanest estimate of the real cost.\n+    let sol = f();\n+    let mut best = f64::INFINITY;\n+    for _ in 0..3 {\n+        let t = Instant::now();\n+        let _ = f();\n+        best = best.min(t.elapsed().as_secs_f64() * 1000.0);\n+    }\n+    Run {\n+        name,\n+        ms: best,\n+        sol,\n+    }\n+}\n+\n+fn run_all(chain: &RankChain, include_dense: bool) -> Vec<Run> {\n+    let opts = SolveOptions {\n+        tol: TOL,\n+        max_iters: CAP,\n+        ..SolveOptions::default()\n+    };\n+    let big_opts = SolveOptions {\n+        tol: TOL,\n+        max_iters: 20_000_000,\n+        ..SolveOptions::default()\n+    };\n+\n+    let mut runs = vec![\n+        time_it(\"power(10k) [current]\", || power(chain, opts)),\n+        time_it(\"power+aitken\", || power_aitken(chain, opts)),\n+        time_it(\"sor(1.0)\", || sor(chain, opts, 1.0)),\n+        time_it(\"bicgstab\", || bicgstab(chain, opts)),\n+        time_it(\"solve() [hybrid]\", || solve(chain, opts)),\n+    ];\n+    if include_dense {\n+        runs.push(time_it(\"sparse-gth(uncapped)\", || force_sparse_gth(chain)));\n+        runs.push(time_it(\"dense-gth\", || dense_gth(chain, TOL)));\n+        runs.push(time_it(\"dense-lu\", || dense_lu(chain, TOL)));\n+    }\n+    // Power iteration with an effectively unlimited cap: what the current\n+    // algorithm would produce if it were allowed to finish.\n+    if chain.n <= 1200 && chain.nnz() <= 8000 {\n+        runs.push(time_it(\"power(unbounded)\", || power(chain, big_opts)));\n+    }\n+    runs\n+}\n+\n+/// Sparse GTH with the budget disabled, for measuring what the direct path\n+/// would cost even on graphs the hybrid would refuse.\n+fn force_sparse_gth(chain: &RankChain) -> Solution {\n+    let unlimited = SolveOptions {\n+        direct_work_budget: u64::MAX,\n+        dense_core_max: usize::MAX,\n+        ..SolveOptions::default()\n+    };\n+    match sparse_gth(chain, unlimited, TOL) {\n+        SparseGthOutcome::Solved(s) => s,\n+        SparseGthOutcome::TooDense { .. } => unreachable!(),\n+    }\n+}\n+\n+fn print_block(label: &str, chain: &RankChain, reference: &[f64], exact_log: Option<&[f64]>) {\n+    println!(\n+        \"\\n--- {label}   n={} nnz={} d_max={}\",\n+        chain.n,\n+        chain.nnz(),\n+        chain.d_max as usize\n+    );\n+    println!(\n+        \"{:<22} {:>9} {:>11} {:>5} {:>7} {:>10} {:>9} {:>10} {:>10} {:>10}\",\n+        \"method\",\n+        \"time(ms)\",\n+        \"resid L1\",\n+        \"conv\",\n+        \"iters\",\n+        \"kendall-d\",\n+        \"log-err\",\n+        \"inv/displ\",\n+        \"tie/maxsh\",\n+        \"kd(log pi)\"\n+    );\n+    for r in run_all(chain, chain.n <= 1500) {\n+        let kd = kendall_distance(&r.sol.pi, reference);\n+        let kd_log = kendall_distance(&r.sol.log_pi, reference);\n+        let (logerr, _) = match exact_log {\n+            Some(e) => max_log_error(&r.sol.pi, e),\n+            None => (f64::NAN, 0usize),\n+        };\n+        let (inv, tie) = chain_order_violations(&r.sol.pi);\n+        let (inv, tie) = if exact_log.is_some() {\n+            // How the ranking would actually be displayed: sort by score\n+            // descending (ties broken by index, as a stable sort does) and see\n+            // how many items land somewhere other than their true rank.\n+            let n = r.sol.pi.len();\n+            let mut order: Vec<usize> = (0..n).collect();\n+            order.sort_by(|&a, &b| {\n+                r.sol.pi[b]\n+                    .partial_cmp(&r.sol.pi[a])\n+                    .unwrap_or(std::cmp::Ordering::Equal)\n+            });\n+            let misplaced = order.iter().enumerate().filter(|&(p, &i)| p != i).count();\n+            let shift = order\n+                .iter()\n+                .enumerate()\n+                .map(|(p, &i)| p.abs_diff(i))\n+                .max()\n+                .unwrap_or(0);\n+            (format!(\"{inv}/{misplaced}\"), format!(\"{tie}/{shift}\"))\n+        } else {\n+            (\"-\".to_string(), \"-\".to_string())\n+        };\n+        println!(\n+            \"{:<22} {:>9.3} {:>11.2e} {:>5} {:>7} {:>10.2e} {:>9} {:>10} {:>10} {:>10.2e}\",\n+            r.name,\n+            r.ms,\n+            r.sol.residual,\n+            if r.sol.converged { \"yes\" } else { \"NO\" },\n+            r.sol.iterations,\n+            kd,\n+            if logerr.is_nan() {\n+                \"-\".to_string()\n+            } else {\n+                format!(\"{logerr:.2e}\")\n+            },\n+            inv,\n+            tie,\n+            kd_log\n+        );\n+    }\n+}\n+\n+fn main() {\n+    let want: Vec<String> = std::env::args().skip(1).collect();\n+    let on = |s: &str| want.is_empty() || want.iter().any(|w| w == s);\n+    println!(\"Rank Centrality stationary solvers — tol={TOL:e}, power cap={CAP}\");\n+\n+    // ---- Chains: ground truth known exactly, and where the cap binds. ----\n+    if on(\"chain\") {\n+        println!(\"\\n================ CHAIN (ratio 2:1, exact order known) ================\");\n+        for n in [100usize, 256, 512, 1024, 2048, 4000] {\n+            let t = Topo::Chain { ratio: 2 };\n+            let chain = build_chain(t, n);\n+            let exact = chain_exact_log(t, n);\n+            // Reference order is the construction order: node 0 best.\n+            let reference: Vec<f64> = (0..n).map(|i| -(i as f64)).collect();\n+            print_block(&format!(\"chain n={n}\"), &chain, &reference, Some(&exact));\n+            println!(\n+                \"   ground truth: π_0/π_{} = 2^{} = 10^{:.0}\",\n+                n - 1,\n+                n - 1,\n+                (n - 1) as f64 * 2f64.log10()\n+            );\n+        }\n+    }\n+    if on(\"varied\") {\n+        println!(\n+            \"\\n================ CHAIN (varied ratios 2..6, benches/common shape) ================\"\n+        );\n+        for n in [256usize, 1024] {\n+            let t = Topo::ChainVaried;\n+            let chain = build_chain(t, n);\n+            let exact = chain_exact_log(t, n);\n+            let reference: Vec<f64> = (0..n).map(|i| -(i as f64)).collect();\n+            print_block(\n+                &format!(\"chain-varied n={n}\"),\n+                &chain,\n+                &reference,\n+                Some(&exact),\n+            );\n+        }\n+\n+        // ---- Other topologies: reference is sparse GTH (subtraction-free). ----\n+    }\n+    if on(\"topo\") {\n+        println!(\"\\n================ OTHER TOPOLOGIES ================\");\n+        let cases: Vec<(String, Topo, usize)> = vec![\n+            (\"star\".into(), Topo::Star, 1024),\n+            (\"clique\".into(), Topo::Clique, 128),\n+            (\"clique\".into(), Topo::Clique, 400),\n+            (\"sparse-d6\".into(), Topo::Sparse { degree: 6 }, 1024),\n+            (\"sparse-d6\".into(), Topo::Sparse { degree: 6 }, 4000),\n+            (\"sparse-d20\".into(), Topo::Sparse { degree: 20 }, 1024),\n+            (\n+                \"components 64x64\".into(),\n+                Topo::Components {\n+                    count: 64,\n+                    size: 64,\n+                },\n+                0,\n+            ),\n+        ];\n+        for (label, t, n) in cases {\n+            let chain = build_chain(t, n);\n+            let reference = solve(\n+                &chain,\n+                SolveOptions {\n+                    tol: TOL,\n+                    max_iters: CAP,\n+                    ..SolveOptions::default()\n+                },\n+            )\n+            .pi;\n+            print_block(\n+                &format!(\"{label} n={}\", node_count(t, n)),\n+                &chain,\n+                &reference,\n+                None,\n+            );\n+        }\n+\n+        // ---- Which method does the hybrid pick, and what does elimination cost? ----\n+    }\n+    if on(\"dispatch\") {\n+        println!(\"\\n================ HYBRID DISPATCH + FILL-IN ================\");\n+        println!(\n+            \"{:<24} {:>7} {:>9} {:>14} {:>10} {:>12}\",\n+            \"graph\", \"n\", \"nnz\", \"picked\", \"iters\", \"resid\"\n+        );\n+        let dispatch: Vec<(String, Topo, usize)> = vec![\n+            (\"chain\".into(), Topo::Chain { ratio: 2 }, 4000),\n+            (\"star\".into(), Topo::Star, 4000),\n+            (\"clique\".into(), Topo::Clique, 200),\n+            (\"clique\".into(), Topo::Clique, 600),\n+            (\"clique\".into(), Topo::Clique, 1200),\n+            (\"sparse-d6\".into(), Topo::Sparse { degree: 6 }, 4000),\n+            (\"sparse-d20\".into(), Topo::Sparse { degree: 20 }, 4000),\n+            (\"sparse-d6\".into(), Topo::Sparse { degree: 6 }, 20000),\n+        ];\n+        for (label, t, n) in dispatch {\n+            let chain = build_chain(t, n);\n+            let opts = SolveOptions {\n+                tol: TOL,\n+                max_iters: CAP,\n+                ..SolveOptions::default()\n+            };\n+            let t0 = Instant::now();\n+            let s = solve(&chain, opts);\n+            let ms = t0.elapsed().as_secs_f64() * 1000.0;\n+            println!(\n+                \"{:<24} {:>7} {:>9} {:>14} {:>10} {:>12.2e}   {:.2} ms {}\",\n+                label,\n+                chain.n,\n+                chain.nnz(),\n+                s.method.label(),\n+                s.iterations,\n+                s.residual,\n+                ms,\n+                if s.converged { \"\" } else { \"NOT CONVERGED\" }\n+            );\n+        }\n+\n+        // ---- Determinism: does HashMap iteration order move the scores? ----\n+    }\n+    if on(\"determinism\") {\n+        println!(\"\\n================ DETERMINISM ================\");\n+        for (label, t, n) in [\n+            (\"sparse-d6 n=1024\", Topo::Sparse { degree: 6 }, 1024usize),\n+            (\"clique n=128\", Topo::Clique, 128),\n+        ] {\n+            let mut first: Option<Vec<f64>> = None;\n+            let mut worst = 0.0f64;\n+            for _ in 0..8 {\n+                let chain = build_chain(t, n);\n+                let pi = force_sparse_gth(&chain).pi;\n+                match &first {\n+                    None => first = Some(pi),\n+                    Some(f) => {\n+                        for i in 0..pi.len() {\n+                            worst = worst.max((pi[i] - f[i]).abs());\n+                        }\n+                    }\n+                }\n+            }\n+            println!(\"{label:<24} max |Δπ| across 8 rebuilds: {worst:.3e}\");\n+        }\n+\n+        // Control: the pre-change builder, rebuilt from scratch each round so\n+        // every rebuild gets a fresh `HashMap` with a different iteration order.\n+        println!(\"\\ncontrol — legacy HashMap-ordered builder, 8 rebuilds of the same graph:\");\n+        for (label, t, n) in [\n+            (\"sparse-d6 n=1024\", Topo::Sparse { degree: 6 }, 1024usize),\n+            (\"clique n=128\", Topo::Clique, 128),\n+            (\"chain n=512\", Topo::Chain { ratio: 2 }, 512),\n+        ] {\n+            let total = node_count(t, n);\n+            let mut g = GroupState::new();\n+            for i in 0..total {\n+                g.ensure_item_pub(&item_name(i));\n+            }\n+            for (k, (i, j, ratio)) in pairs(t, total).into_iter().enumerate() {\n+                g.apply_vote(vote(i, j, ratio, 1, k as i64));\n+            }\n+            let opts = SolveOptions {\n+                tol: TOL,\n+                max_iters: CAP,\n+                ..SolveOptions::default()\n+            };\n+            let mut reference: Option<Vec<f64>> = None;\n+            let mut worst_abs = 0.0f64;\n+            let mut worst_rel = 0.0f64;\n+            let mut flips = 0usize;\n+            let mut distinct = 0usize;\n+            for _ in 0..8 {\n+                let c = legacy_chain(total, g.edges.iter().map(|(&k, &w)| (k, w)));\n+                let pi = power(&c, opts).pi;\n+                match &reference {\n+                    None => reference = Some(pi),\n+                    Some(r) => {\n+                        if pi != *r {\n+                            distinct += 1;\n+                        }\n+                        for i in 0..pi.len() {\n+                            worst_abs = worst_abs.max((pi[i] - r[i]).abs());\n+                            if r[i] > 0.0 {\n+                                worst_rel = worst_rel.max((pi[i] - r[i]).abs() / r[i]);\n+                            }\n+                        }\n+                        flips += order_flips(&pi, r);\n+                    }\n+                }\n+            }\n+            println!(\n+                \"{label:<20} differing rebuilds {distinct}/7  max |Δπ|={worst_abs:.2e}  \\\n+                 max rel={worst_rel:.2e}  rank flips {flips}\"\n+            );\n+        }\n+\n+        // Cost of the sort that buys determinism, isolated from the `HashMap`\n+        // aggregation `chain_from_edges` has always done.\n+        println!(\"\\nsort cost (edge ordering inside chain_from_edges):\");\n+        for (label, t, n) in [\n+            (\"sparse-d6 n=4000\", Topo::Sparse { degree: 6 }, 4000usize),\n+            (\"clique n=400\", Topo::Clique, 400),\n+        ] {\n+            let total = node_count(t, n);\n+            let mut g = GroupState::new();\n+            for i in 0..total {\n+                g.ensure_item_pub(&item_name(i));\n+            }\n+            for (k, (i, j, ratio)) in pairs(t, total).into_iter().enumerate() {\n+                g.apply_vote(vote(i, j, ratio, 1, k as i64));\n+            }\n+            let t0 = Instant::now();\n+            for _ in 0..20 {\n+                let c = chain_from_edges(total, g.edges.iter().map(|(&k, &w)| (k, w)));\n+                std::hint::black_box(c.nnz());\n+            }\n+            let with_sort = t0.elapsed().as_secs_f64() * 1000.0 / 20.0;\n+\n+            let edges: Vec<((usize, usize), f64)> = g.edges.iter().map(|(&k, &w)| (k, w)).collect();\n+            let t0 = Instant::now();\n+            for _ in 0..20 {\n+                let mut keys: Vec<(usize, usize)> = edges.iter().map(|&(k, _)| k).collect();\n+                keys.sort_unstable();\n+                std::hint::black_box(keys.len());\n+            }\n+            let sort_only = t0.elapsed().as_secs_f64() * 1000.0 / 20.0;\n+            println!(\n+                \"{label:<20} chain_from_edges {with_sort:7.3} ms, of which sorting is \\\n+                 {sort_only:6.3} ms  (E={})\",\n+                g.edges.len()\n+            );\n+        }\n+\n+        // ---- Cross-method equivalence on small graphs, at machine precision. ----\n+    }\n+    if on(\"equiv\") {\n+        println!(\"\\n================ EQUIVALENCE (small n, all methods) ================\");\n+        for (label, t, n) in [\n+            (\"chain n=64\", Topo::Chain { ratio: 2 }, 64usize),\n+            (\"star n=64\", Topo::Star, 64),\n+            (\"clique n=64\", Topo::Clique, 64),\n+            (\"sparse-d6 n=200\", Topo::Sparse { degree: 6 }, 200),\n+        ] {\n+            let chain = build_chain(t, n);\n+            let gth = dense_gth(&chain, TOL).pi;\n+            let mut worst: Vec<(String, f64)> = Vec::new();\n+            let opts = SolveOptions {\n+                tol: 1e-14,\n+                max_iters: 200_000,\n+                ..SolveOptions::default()\n+            };\n+            let cands: Vec<(&str, Vec<f64>)> = vec![\n+                (\"sparse-gth\", force_sparse_gth(&chain).pi),\n+                (\"dense-lu\", dense_lu(&chain, TOL).pi),\n+                (\"power(200k)\", power(&chain, opts).pi),\n+                (\"sor\", sor(&chain, opts, 1.0).pi),\n+                (\"bicgstab\", bicgstab(&chain, opts).pi),\n+            ];\n+            for (name, pi) in cands {\n+                let rel = (0..chain.n)\n+                    .filter(|&i| gth[i] > 1e-300)\n+                    .map(|i| (pi[i] - gth[i]).abs() / gth[i])\n+                    .fold(0.0f64, f64::max);\n+                worst.push((name.to_string(), rel));\n+            }\n+            let cells: Vec<String> = worst.iter().map(|(n, v)| format!(\"{n}={v:.1e}\")).collect();\n+            println!(\"{label:<18} max rel err vs dense-GTH: {}\", cells.join(\"  \"));\n+        }\n+    }\n+    let _ = Method::Power;\n+}\ndiff --git a/server/src/lib.rs b/server/src/lib.rs\nindex e102e73ec167fa35dcc22d4ec8c6f40202982429..7845de0b8512738fe20ca9dffd84e9b8c8134d50 100644\n--- a/server/src/lib.rs\n+++ b/server/src/lib.rs\n@@ -16,6 +16,7 @@ pub mod ranking;\n pub mod reducer;\n pub mod scope_rank;\n pub mod state;\n+pub mod stationary;\n pub mod timeago;\n pub mod views;\n pub mod write_cmd;\ndiff --git a/server/src/ranking.rs b/server/src/ranking.rs\nindex 38e6d09b4370e5f8cbae09c0e5760b4e7f1ef7db..2d391bb6fdf4d22247e5524d459f2d1609bd6d70 100644\n--- a/server/src/ranking.rs\n+++ b/server/src/ranking.rs\n@@ -2,6 +2,7 @@ use std::collections::{HashMap, HashSet};\n \n use crate::path_types::ItemId;\n use crate::reducer::GroupState;\n+use crate::stationary::{self, RankChain, SolveOptions, Solution};\n \n #[derive(Debug, Clone)]\n pub struct RankedItem {\n@@ -105,14 +106,21 @@ pub fn ranked_items(group: &mut GroupState, max_iters: usize, tol: f64) -> Vec<R\n     items\n }\n \n-fn compute_scores_from_edges(n: usize, edges: impl Iterator<Item = ((usize, usize), f64)>, max_iters: usize, tol: f64) -> Vec<f64> {\n-    if n == 0 {\n-        return vec![];\n-    }\n-    if n == 1 {\n-        return vec![1.0];\n-    }\n-\n+/// Build the Rank Centrality Markov chain from raw directed vote weights.\n+///\n+/// Rank Centrality (Negahban, Oh, Shah 2012, §3.1):\n+///   P_ij = (1/d_max) * a_ij           for i ≠ j compared\n+///   P_ii = 1 - (1/d_max) * Σ_k a_ik\n+/// where d_i is the *degree* (number of distinct neighbors compared) and\n+/// d_max = max_i d_i. Using the unweighted degree — not the sum of\n+/// pairwise-normalized weights — is what guarantees aperiodicity: it\n+/// forces P_ii > 0 for every non-maximum-degree node, and for max-degree\n+/// nodes whenever any neighbor weight is below 1 (i.e. not a unanimous\n+/// loss). Without this, regular comparison graphs (e.g. a pure star at\n+/// ratio 2:1) produce a bipartite chain that oscillates instead of\n+/// converging — see issue #146. (Direct solvers are immune either way:\n+/// π depends only on the off-diagonals, which d_max scales uniformly.)\n+pub fn chain_from_edges(n: usize, edges: impl Iterator<Item = ((usize, usize), f64)>) -> RankChain {\n     // Collect raw edges into a map for pairwise normalization.\n     let mut raw: HashMap<(usize, usize), f64> = HashMap::new();\n     for ((src, dst), w) in edges {\n@@ -125,90 +133,80 @@ fn compute_scores_from_edges(n: usize, edges: impl Iterator<Item = ((usize, usiz\n     // Pairwise normalization: a_ij = A_ij / (A_ij + A_ji).\n     // This ensures repeated votes on the same pair don't inflate influence\n     // beyond what the ratio implies.\n-    let keys: Vec<(usize, usize)> = raw.keys().copied().collect();\n-    let mut normalized: HashMap<(usize, usize), f64> = HashMap::new();\n+    //\n+    // Sorted, not read straight off the `HashMap`: the summation order of every\n+    // downstream reduction has to be a function of the graph alone, or scores\n+    // wobble in their low bits between runs and near-ties can flip.\n+    let mut keys: Vec<(usize, usize)> = raw.keys().copied().collect();\n+    keys.sort_unstable();\n+\n+    let mut normalized: Vec<((usize, usize), f64)> = Vec::with_capacity(keys.len());\n+    let mut neighbors: Vec<HashSet<usize>> = vec![HashSet::new(); n];\n     for (i, j) in keys {\n-        if normalized.contains_key(&(i, j)) {\n-            continue;\n-        }\n         let w_ij = *raw.get(&(i, j)).unwrap_or(&0.0);\n         let w_ji = *raw.get(&(j, i)).unwrap_or(&0.0);\n         let total = w_ij + w_ji;\n         if total <= 0.0 {\n             continue;\n         }\n-        normalized.insert((i, j), w_ij / total);\n-        if w_ji > 0.0 {\n-            normalized.insert((j, i), w_ji / total);\n-        }\n+        normalized.push(((i, j), w_ij / total));\n+        neighbors[i].insert(j);\n+        neighbors[j].insert(i);\n     }\n \n-    // Rank Centrality (Negahban, Oh, Shah 2012, §3.1):\n-    //   P_ij = (1/d_max) * A_ij           for i ≠ j compared\n-    //   P_ii = 1 - (1/d_max) * Σ_k A_ik\n-    // where d_i is the *degree* (number of distinct neighbors compared) and\n-    // d_max = max_i d_i. Using the unweighted degree — not the sum of\n-    // pairwise-normalized weights — is what guarantees aperiodicity: it\n-    // forces P_ii > 0 for every non-maximum-degree node, and for max-degree\n-    // nodes whenever any neighbor weight is below 1 (i.e. not a unanimous\n-    // loss). Without this, regular comparison graphs (e.g. a pure star at\n-    // ratio 2:1) produce a bipartite chain that oscillates instead of\n-    // converging — see issue #146.\n-    let mut out_edges: Vec<Vec<(usize, f64)>> = vec![Vec::new(); n];\n-    let mut neighbors: Vec<HashSet<usize>> = vec![HashSet::new(); n];\n+    let d_max = neighbors.iter().map(|s| s.len()).max().unwrap_or(0);\n+    RankChain::from_normalized(n, normalized, d_max)\n+}\n \n-    for ((src, dst), w) in &normalized {\n-        out_edges[*src].push((*dst, *w));\n-        neighbors[*src].insert(*dst);\n-        neighbors[*dst].insert(*src);\n+/// Stationary distribution for the chain induced by `edges`, with convergence\n+/// diagnostics attached. See [`stationary::Solution`].\n+///\n+/// This is the API that replaces \"return whatever the iteration reached\":\n+/// callers that care can inspect `converged`, `residual` and `method`, and\n+/// non-convergence is logged rather than swallowed.\n+pub fn solve_scores_from_edges(\n+    n: usize,\n+    edges: impl Iterator<Item = ((usize, usize), f64)>,\n+    max_iters: usize,\n+    tol: f64,\n+) -> Solution {\n+    if n == 0 {\n+        return stationary::trivial(vec![]);\n     }\n-\n-    let weight_sum: Vec<f64> = out_edges\n-        .iter()\n-        .map(|es| es.iter().map(|(_, w)| *w).sum())\n-        .collect();\n-    let d_max = neighbors.iter().map(|s| s.len()).max().unwrap_or(0);\n-    if d_max == 0 {\n-        return vec![1.0 / n as f64; n];\n+    if n == 1 {\n+        return stationary::trivial(vec![1.0]);\n     }\n-    let d_max_f = d_max as f64;\n-\n-    let mut scores = vec![1.0 / n as f64; n];\n-    let mut next = vec![0.0f64; n];\n-\n-    for _ in 0..max_iters {\n-        next.fill(0.0);\n-        for i in 0..n {\n-            let stay_prob = (d_max_f - weight_sum[i]) / d_max_f;\n-            next[i] += scores[i] * stay_prob;\n-\n-            if out_edges[i].is_empty() {\n-                continue;\n-            }\n-            for &(dst, w) in &out_edges[i] {\n-                next[dst] += scores[i] * (w / d_max_f);\n-            }\n-        }\n \n-        let diff: f64 = scores\n-            .iter()\n-            .zip(next.iter())\n-            .map(|(a, b)| (a - b).abs())\n-            .sum();\n+    let chain = chain_from_edges(n, edges);\n+    let solution = stationary::solve(\n+        &chain,\n+        SolveOptions {\n+            tol,\n+            max_iters,\n+            ..SolveOptions::default()\n+        },\n+    );\n \n-        scores.clone_from_slice(&next);\n-        if diff < tol {\n-            break;\n-        }\n+    if !solution.converged {\n+        tracing::warn!(\n+            n,\n+            method = solution.method.label(),\n+            iterations = solution.iterations,\n+            residual = solution.residual,\n+            \"rank centrality did not reach tolerance; this ranking may be misordered\"\n+        );\n+    } else if solution.underflowed {\n+        tracing::debug!(\n+            n,\n+            method = solution.method.label(),\n+            \"rank centrality scores span more than f64 holds; use log scores to order the tail\"\n+        );\n     }\n+    solution\n+}\n \n-    let sum: f64 = scores.iter().sum();\n-    if sum.is_finite() && sum > 0.0 {\n-        for s in &mut scores {\n-            *s /= sum;\n-        }\n-    }\n-    scores\n+fn compute_scores_from_edges(n: usize, edges: impl Iterator<Item = ((usize, usize), f64)>, max_iters: usize, tol: f64) -> Vec<f64> {\n+    solve_scores_from_edges(n, edges, max_iters, tol).pi\n }\n \n /// Rank-centrality within a subset of items (an induced subgraph), using the group's aggregated edges.\n@@ -231,20 +229,25 @@ pub fn ranked_items_subset(group: &GroupState, idxs: &[usize], max_iters: usize,\n         Some(((s, d), w))\n     });\n \n-    let scores = compute_scores_from_edges(idxs.len(), edges_iter, max_iters, tol);\n+    let solved = solve_scores_from_edges(idxs.len(), edges_iter, max_iters, tol);\n \n     // Filter out entries where idx_to_item doesn't have the slot (shouldn't happen, but be safe).\n-    let mut items: Vec<RankedItem> = idxs\n+    // Sorting on the *log* score, not `score`: on a long preference chain the\n+    // true distribution spans more decades than f64 holds, so `score` ties off\n+    // at the bottom while the log scores still order it correctly.\n+    let mut items: Vec<(RankedItem, f64)> = idxs\n         .iter()\n         .enumerate()\n         .filter_map(|(j, &orig)| {\n             let item = group.idx_to_item.get(orig)?.clone();\n-            Some(RankedItem { item, score: *scores.get(j).unwrap_or(&0.0) })\n+            let score = *solved.pi.get(j).unwrap_or(&0.0);\n+            let log_score = *solved.log_pi.get(j).unwrap_or(&f64::NEG_INFINITY);\n+            Some((RankedItem { item, score }, log_score))\n         })\n         .collect();\n \n-    items.sort_by(|a, b| b.score.partial_cmp(&a.score).unwrap_or(std::cmp::Ordering::Equal));\n-    items\n+    items.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap_or(std::cmp::Ordering::Equal));\n+    items.into_iter().map(|(item, _)| item).collect()\n }\n \n pub fn group_summary_scores(\n@@ -314,6 +317,82 @@ mod tests {\n         );\n     }\n \n+    /// A long preference chain — the natural shape of a deep ontology or a\n+    /// hand-ordered list — is where power iteration got the answer wrong. With\n+    /// varied ratios and 1200 items the old solver reported success (its L1\n+    /// step fell below 1e-8 after 4325 sweeps, well inside the cap) while\n+    /// leaving 47% of items at the wrong rank, because an absolute step\n+    /// tolerance says nothing about entries that are themselves 1e-300.\n+    #[test]\n+    fn deep_chain_is_ranked_in_exactly_the_right_order() {\n+        let n = 1200usize;\n+        let mut g = mk_group();\n+        let name = |i: usize| format!(\"i{i:05}\");\n+        for i in 0..n - 1 {\n+            let left = 2 + (i % 5) as i32;\n+            g.apply_vote(vote(i as i64, &name(i), &name(i + 1), left, 1));\n+        }\n+\n+        let idxs: Vec<usize> = (0..g.idx_to_item.len()).collect();\n+        let ranked = ranked_items_subset(&g, &idxs, 10000, 1e-8);\n+        assert_eq!(ranked.len(), n);\n+        for (i, r) in ranked.iter().enumerate() {\n+            let want = format!(\"https://slug.social/{}\", name(i));\n+            assert_eq!(r.item.as_str(), want, \"position {i} of the chain ranking\");\n+        }\n+    }\n+\n+    /// Same graph, edges presented in a different order: the scores must come\n+    /// back bit-identical, not merely close.\n+    #[test]\n+    fn ranking_does_not_depend_on_edge_iteration_order() {\n+        let n = 200usize;\n+        let mut edges: Vec<((usize, usize), f64)> = Vec::new();\n+        for i in 0..n - 1 {\n+            edges.push(((i + 1, i), 3.0));\n+            edges.push(((i, i + 1), 1.0));\n+        }\n+        for i in 0..n / 4 {\n+            let j = (5 * i + 7) % n;\n+            if i != j {\n+                edges.push(((j, i), 2.0));\n+                edges.push(((i, j), 1.0));\n+            }\n+        }\n+\n+        let reference = compute_scores_from_edges(n, edges.iter().copied(), 10000, 1e-8);\n+        let mut state = 0x5EEDu64;\n+        for _ in 0..5 {\n+            let mut shuffled = edges.clone();\n+            for k in (1..shuffled.len()).rev() {\n+                state ^= state >> 12;\n+                state ^= state << 25;\n+                state ^= state >> 27;\n+                let m = (state.wrapping_mul(0x2545_F491_4F6C_DD1D) % (k as u64 + 1)) as usize;\n+                shuffled.swap(k, m);\n+            }\n+            let got = compute_scores_from_edges(n, shuffled.into_iter(), 10000, 1e-8);\n+            assert_eq!(got, reference);\n+        }\n+    }\n+\n+    /// Convergence is reported, never assumed.\n+    #[test]\n+    fn solver_reports_convergence_for_every_shape() {\n+        for n in [2usize, 3, 64, 900] {\n+            let mut edges: Vec<((usize, usize), f64)> = Vec::new();\n+            for i in 0..n - 1 {\n+                edges.push(((i + 1, i), 2.0));\n+                edges.push(((i, i + 1), 1.0));\n+            }\n+            let solved = solve_scores_from_edges(n, edges.into_iter(), 10000, 1e-8);\n+            assert!(solved.converged, \"chain n={n} reported non-convergence\");\n+            assert!(solved.residual < 1e-10, \"chain n={n}: {:e}\", solved.residual);\n+            assert_eq!(solved.pi.len(), n);\n+            assert_eq!(solved.log_pi.len(), n);\n+        }\n+    }\n+\n     #[test]\n     fn group_ranking_cache_dirty_flow() {\n         let mut g = mk_group();\ndiff --git a/server/src/stationary.rs b/server/src/stationary.rs\nnew file mode 100644\nindex 0000000000000000000000000000000000000000..f13af1526e7abdf53fcc02a6f38dee8c6c2257fe\n--- /dev/null\n+++ b/server/src/stationary.rs\n@@ -0,0 +1,1623 @@\n+//! Stationary-distribution solvers for the Rank Centrality Markov chain.\n+//!\n+//! `ranking::compute_scores_from_edges` builds a chain over compared items and\n+//! needs π with `πP = π`, `Σπ = 1`. This module owns every candidate solver so\n+//! they can be compared against each other on identical input, and so the\n+//! ranking code has exactly one place to ask for \"the stationary distribution,\n+//! and tell me whether you actually got there\".\n+//!\n+//! # The chain\n+//!\n+//! Rank Centrality (Negahban, Oh, Shah 2012, §3.1) uses\n+//!\n+//! ```text\n+//!   P_ij = a_ij / d_max        (i ≠ j, compared)\n+//!   P_ii = 1 - (Σ_k a_ik) / d_max\n+//! ```\n+//!\n+//! where `a_ij = A_ij / (A_ij + A_ji)` and `d_max` is the maximum *unweighted*\n+//! degree. `d_max` is a uniformization constant: it exists only to keep the\n+//! chain aperiodic so power iteration cannot oscillate (issue #146). Because\n+//!\n+//! ```text\n+//!   πP = π   ⟺   πQ = 0,   Q = P - I,   Q_ij = a_ij / d_max  (i ≠ j)\n+//! ```\n+//!\n+//! and `Q` is only defined by its off-diagonal entries (the diagonal is minus\n+//! the row sum), scaling every off-diagonal by the same `1/d_max` leaves π\n+//! unchanged. **Direct solvers therefore do not need `d_max` at all** and are\n+//! structurally immune to the periodicity bug. Only the iterative solvers,\n+//! which literally walk `P`, care.\n+//!\n+//! # Dynamic range\n+//!\n+//! On a tree (a chain is a tree) the chain is reversible and detailed balance\n+//! pins the answer exactly: `π_i / π_j = a_ji / a_ij` for every edge. A chain\n+//! of `n` items each preferred 2:1 over the next therefore has\n+//! `π_max / π_min = 2^(n-1)`, which leaves the f64 range at n ≈ 1075. This is a\n+//! property of the model, not of any solver: past that width the tail of the\n+//! ranking is not representable in double precision and underflows to zero.\n+//! [`Solution::underflowed`] reports it, and [`Solution::log_pi`] stays exact\n+//! there, which is why it — not `pi` — is the sort key for a ranking.\n+//!\n+//! # Strategy ([`solve`])\n+//!\n+//! 1. Split disconnected components; each is solved on its own and keeps the\n+//!    share of the mass its node count started with (what power iteration from\n+//!    a uniform start converges to).\n+//! 2. Sparse GTH state reduction, minimum-degree order, with a degree guard and\n+//!    a work budget. Trees and chains reduce completely in `O(n)`; whatever\n+//!    survives is an irreducible core solved by dense GTH when it fits.\n+//! 3. Otherwise Gauss–Seidel, then power iteration as a second opinion. Graphs\n+//!    that defeat step 2 are exactly the well-connected ones these finish in\n+//!    tens of sweeps.\n+//!\n+//! Nothing returns without a checked residual: see [`Solution::converged`].\n+\n+/// Which solver produced a [`Solution`].\n+#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]\n+pub enum Method {\n+    /// Baseline: repeated `π ← πP` until the L1 step falls under `tol`.\n+    Power,\n+    /// Power iteration with periodic Aitken Δ² extrapolation.\n+    PowerAitken,\n+    /// Gauss–Seidel / SOR sweeps on the balance equations.\n+    Sor,\n+    /// Dense LU with partial pivoting on `Qᵀπ = e`, one row replaced by `Σπ = 1`.\n+    DenseLu,\n+    /// Dense Grassmann–Taksar–Heyman state reduction, `O(n³)`, subtraction-free.\n+    DenseGth,\n+    /// GTH state reduction over a sparse graph with minimum-degree elimination.\n+    SparseGth,\n+    /// Jacobi-preconditioned BiCGSTAB on the singular balance system.\n+    BiCgStab,\n+}\n+\n+impl Method {\n+    pub fn label(self) -> &'static str {\n+        match self {\n+            Method::Power => \"power\",\n+            Method::PowerAitken => \"power+aitken\",\n+            Method::Sor => \"sor\",\n+            Method::DenseLu => \"dense-lu\",\n+            Method::DenseGth => \"dense-gth\",\n+            Method::SparseGth => \"sparse-gth\",\n+            Method::BiCgStab => \"bicgstab\",\n+        }\n+    }\n+\n+    /// Direct methods answer exactly (up to rounding); iterative ones can stall.\n+    pub fn is_direct(self) -> bool {\n+        matches!(self, Method::DenseLu | Method::DenseGth | Method::SparseGth)\n+    }\n+}\n+\n+/// A stationary distribution plus everything needed to judge whether to trust it.\n+#[derive(Debug, Clone)]\n+pub struct Solution {\n+    /// Non-negative, sums to 1 (unless every entry underflowed).\n+    pub pi: Vec<f64>,\n+    /// `ln π_i`, shifted so the largest entry is 0. Direct solvers fill this in\n+    /// exactly even where `pi` has flushed to zero, so it is the only faithful\n+    /// sort key on graphs whose score spread exceeds the f64 range.\n+    /// `f64::NEG_INFINITY` for genuinely unreachable states.\n+    pub log_pi: Vec<f64>,\n+    pub method: Method,\n+    /// Sweeps performed. Zero for direct methods.\n+    pub iterations: usize,\n+    /// `‖πP − π‖₁`, the honest backward error of the answer that was returned.\n+    pub residual: f64,\n+    /// False means the number that came back is *not* the stationary\n+    /// distribution to the requested tolerance. Never silently true.\n+    pub converged: bool,\n+    /// The exact answer spans more than f64 can hold, so the tail of the\n+    /// ranking has been flushed to zero and its internal order is lost.\n+    pub underflowed: bool,\n+}\n+\n+/// Tuning for [`solve`]. Defaults encode the recommended hybrid strategy.\n+#[derive(Debug, Clone, Copy)]\n+pub struct SolveOptions {\n+    pub tol: f64,\n+    pub max_iters: usize,\n+    /// Row-merge budget for the sparse elimination in [`sparse_gth`]. Bounds\n+    /// the work wasted before giving up on the direct path.\n+    pub direct_work_budget: u64,\n+    /// Largest irreducible core still worth an `O(n³)` dense GTH.\n+    pub dense_core_max: usize,\n+}\n+\n+impl Default for SolveOptions {\n+    fn default() -> Self {\n+        SolveOptions {\n+            tol: 1e-8,\n+            max_iters: 10_000,\n+            direct_work_budget: DEFAULT_DIRECT_WORK_BUDGET,\n+            dense_core_max: DEFAULT_DENSE_CORE_MAX,\n+        }\n+    }\n+}\n+\n+/// Measured: sparse elimination runs at roughly 20–40 M row-merge steps per\n+/// second, so 2 M caps the abandoned work at a few tens of milliseconds. Trees\n+/// and chains finish two to three orders of magnitude under it.\n+pub const DEFAULT_DIRECT_WORK_BUDGET: u64 = 2_000_000;\n+\n+/// Dense GTH is `O(n³)`; measured at ~5 ms for n=400 and ~40 ms for n=800, so\n+/// 512 keeps the direct path inside a garden render's budget. Above it, the\n+/// iterative solvers are both faster and (on well-connected graphs) accurate.\n+pub const DEFAULT_DENSE_CORE_MAX: usize = 512;\n+\n+/// The comparison chain in the one form every solver consumes.\n+///\n+/// `rows[i]` holds `(j, a_ij)` for the pairwise-normalized weights, sorted by\n+/// `j` so that every floating-point summation in this module runs in a fixed\n+/// order regardless of how the caller's `HashMap` happened to iterate.\n+#[derive(Debug, Clone)]\n+pub struct RankChain {\n+    pub n: usize,\n+    pub rows: Vec<Vec<(usize, f64)>>,\n+    pub row_sum: Vec<f64>,\n+    /// Maximum unweighted degree; the uniformization constant for `P`.\n+    pub d_max: f64,\n+}\n+\n+impl RankChain {\n+    /// Build from pairwise-normalized weights. `edges` may arrive in any order;\n+    /// they are sorted here so results are bit-for-bit reproducible, and\n+    /// repeated `(i, j)` pairs are coalesced so every row holds each column\n+    /// exactly once (elimination relies on that).\n+    pub fn from_normalized(n: usize, mut edges: Vec<((usize, usize), f64)>, d_max: usize) -> Self {\n+        edges.sort_unstable_by(|a, b| a.0.cmp(&b.0));\n+        let mut rows: Vec<Vec<(usize, f64)>> = vec![Vec::new(); n];\n+        for ((i, j), w) in edges {\n+            // Written positively so a NaN weight is dropped rather than kept.\n+            let usable = i < n && j < n && i != j && w > 0.0;\n+            if !usable {\n+                continue;\n+            }\n+            match rows[i].last_mut() {\n+                Some(last) if last.0 == j => last.1 += w,\n+                _ => rows[i].push((j, w)),\n+            }\n+        }\n+        let row_sum: Vec<f64> = rows\n+            .iter()\n+            .map(|r| r.iter().map(|&(_, w)| w).sum())\n+            .collect();\n+        RankChain {\n+            n,\n+            rows,\n+            row_sum,\n+            d_max: d_max as f64,\n+        }\n+    }\n+\n+    pub fn nnz(&self) -> usize {\n+        self.rows.iter().map(|r| r.len()).sum()\n+    }\n+\n+    /// Incoming edges, `cols[j] = [(i, a_ij)]`, sorted by `i`.\n+    fn columns(&self) -> Vec<Vec<(usize, f64)>> {\n+        let mut cols: Vec<Vec<(usize, f64)>> = vec![Vec::new(); self.n];\n+        for i in 0..self.n {\n+            for &(j, w) in &self.rows[i] {\n+                cols[j].push((i, w));\n+            }\n+        }\n+        cols\n+    }\n+\n+    /// `‖πP − π‖₁` for an arbitrary vector, using the uniformized `P`.\n+    ///\n+    /// This is the yardstick every method is judged by, so it must not depend\n+    /// on how the candidate was produced.\n+    pub fn residual(&self, pi: &[f64]) -> f64 {\n+        let n = self.n;\n+        let mut next = vec![0.0f64; n];\n+        for i in 0..n {\n+            let stay = (self.d_max - self.row_sum[i]) / self.d_max;\n+            next[i] += pi[i] * stay;\n+            for &(j, w) in &self.rows[i] {\n+                next[j] += pi[i] * (w / self.d_max);\n+            }\n+        }\n+        (0..n).map(|i| (next[i] - pi[i]).abs()).sum()\n+    }\n+}\n+\n+/// A `Solution` for an input with nothing to solve (no items, or one).\n+pub fn trivial(pi: Vec<f64>) -> Solution {\n+    Solution {\n+        log_pi: vec![0.0; pi.len()],\n+        pi,\n+        method: Method::SparseGth,\n+        iterations: 0,\n+        residual: 0.0,\n+        converged: true,\n+        underflowed: false,\n+    }\n+}\n+\n+/// Normalize in place to sum 1. Returns whether the tail underflowed to zero.\n+fn normalize(pi: &mut [f64]) -> bool {\n+    let sum: f64 = pi.iter().sum();\n+    if sum.is_finite() && sum > 0.0 {\n+        for p in pi.iter_mut() {\n+            *p /= sum;\n+        }\n+    } else {\n+        let uniform = 1.0 / pi.len() as f64;\n+        for p in pi.iter_mut() {\n+            *p = uniform;\n+        }\n+        return true;\n+    }\n+    pi.iter().any(|&p| p <= 0.0)\n+}\n+\n+fn finish(\n+    pi: Vec<f64>,\n+    chain: &RankChain,\n+    method: Method,\n+    iterations: usize,\n+    tol: f64,\n+) -> Solution {\n+    let log_pi = log_of(&pi);\n+    finish_with_logs(pi, log_pi, chain, method, iterations, tol)\n+}\n+\n+fn log_of(pi: &[f64]) -> Vec<f64> {\n+    let max = pi.iter().cloned().fold(0.0f64, f64::max);\n+    if max <= 0.0 {\n+        return vec![f64::NEG_INFINITY; pi.len()];\n+    }\n+    pi.iter().map(|&p| (p / max).ln()).collect()\n+}\n+\n+fn finish_with_logs(\n+    pi: Vec<f64>,\n+    log_pi: Vec<f64>,\n+    chain: &RankChain,\n+    method: Method,\n+    iterations: usize,\n+    tol: f64,\n+) -> Solution {\n+    let mut pi = pi;\n+    let underflowed = normalize(&mut pi);\n+    let residual = chain.residual(&pi);\n+    Solution {\n+        converged: residual.is_finite() && residual <= tol.max(f64::EPSILON * chain.n as f64),\n+        pi,\n+        log_pi,\n+        method,\n+        iterations,\n+        residual,\n+        underflowed,\n+    }\n+}\n+\n+/// `ln(Σ exp(t))` computed by shifting out the largest term, so a sum of\n+/// astronomically different magnitudes neither overflows nor loses the small ones.\n+fn log_sum_exp(terms: &[f64]) -> f64 {\n+    let max = terms.iter().cloned().fold(f64::NEG_INFINITY, f64::max);\n+    if !max.is_finite() {\n+        return max;\n+    }\n+    let sum: f64 = terms.iter().map(|&t| (t - max).exp()).sum();\n+    max + sum.ln()\n+}\n+\n+/// Rebuild `pi` from exact log-scores, shifted so the maximum is 1.\n+///\n+/// `exp` degrades into subnormals rather than jumping to zero, so this keeps\n+/// roughly 745 decades of usable spread where the linear back-substitution\n+/// keeps 308 — the difference between an exact ranking of a 2400-long chain and\n+/// an exact ranking of a 1000-long one.\n+fn pi_from_logs(log_pi: &[f64]) -> Vec<f64> {\n+    let max = log_pi.iter().cloned().fold(f64::NEG_INFINITY, f64::max);\n+    if !max.is_finite() {\n+        return vec![1.0 / log_pi.len() as f64; log_pi.len()];\n+    }\n+    log_pi.iter().map(|&l| (l - max).exp()).collect()\n+}\n+\n+// ---------------------------------------------------------------------------\n+// Iterative solvers\n+// ---------------------------------------------------------------------------\n+\n+/// One `π ← πP` sweep.\n+fn step(chain: &RankChain, pi: &[f64], next: &mut [f64]) {\n+    next.fill(0.0);\n+    for i in 0..chain.n {\n+        let stay = (chain.d_max - chain.row_sum[i]) / chain.d_max;\n+        next[i] += pi[i] * stay;\n+        for &(j, w) in &chain.rows[i] {\n+            next[j] += pi[i] * (w / chain.d_max);\n+        }\n+    }\n+}\n+\n+/// Baseline power iteration: exactly what `ranking.rs` does today.\n+pub fn power(chain: &RankChain, opts: SolveOptions) -> Solution {\n+    let n = chain.n;\n+    let mut pi = vec![1.0 / n as f64; n];\n+    let mut next = vec![0.0f64; n];\n+    let mut iters = 0usize;\n+\n+    for _ in 0..opts.max_iters {\n+        step(chain, &pi, &mut next);\n+        iters += 1;\n+        let diff: f64 = (0..n).map(|i| (pi[i] - next[i]).abs()).sum();\n+        pi.copy_from_slice(&next);\n+        if diff < opts.tol {\n+            break;\n+        }\n+    }\n+    finish(pi, chain, Method::Power, iters, opts.tol)\n+}\n+\n+/// Power iteration with periodic componentwise Aitken Δ² extrapolation.\n+///\n+/// Aitken assumes the error is dominated by a single geometric mode, which is\n+/// exactly the regime that makes plain power iteration slow. The extrapolated\n+/// vector is only accepted when it lowers the residual, so a bad extrapolation\n+/// costs one wasted sweep instead of divergence.\n+pub fn power_aitken(chain: &RankChain, opts: SolveOptions) -> Solution {\n+    let n = chain.n;\n+    let mut x0 = vec![1.0 / n as f64; n];\n+    let mut x1 = vec![0.0f64; n];\n+    let mut x2 = vec![0.0f64; n];\n+    let mut cand = vec![0.0f64; n];\n+    let mut iters = 0usize;\n+    let mut best_resid = f64::INFINITY;\n+\n+    while iters < opts.max_iters {\n+        step(chain, &x0, &mut x1);\n+        step(chain, &x1, &mut x2);\n+        iters += 2;\n+\n+        let diff: f64 = (0..n).map(|i| (x2[i] - x1[i]).abs()).sum();\n+        if diff < opts.tol {\n+            return finish(x2, chain, Method::PowerAitken, iters, opts.tol);\n+        }\n+\n+        // Aitken: x* ≈ x2 - (Δx1)² / Δ²x0, componentwise.\n+        let mut usable = true;\n+        for i in 0..n {\n+            let d1 = x1[i] - x0[i];\n+            let d2 = x2[i] - x1[i];\n+            let denom = d2 - d1;\n+            if denom.abs() <= f64::MIN_POSITIVE {\n+                cand[i] = x2[i];\n+                continue;\n+            }\n+            let v = x2[i] - d2 * d2 / denom;\n+            if !v.is_finite() || v < 0.0 {\n+                usable = false;\n+                break;\n+            }\n+            cand[i] = v;\n+        }\n+\n+        std::mem::swap(&mut x0, &mut x2);\n+        if !usable {\n+            continue;\n+        }\n+        let sum: f64 = cand.iter().sum();\n+        if !(sum.is_finite() && sum > 0.0) {\n+            continue;\n+        }\n+        for c in cand.iter_mut() {\n+            *c /= sum;\n+        }\n+        let r_cand = chain.residual(&cand);\n+        let r_cur = chain.residual(&x0);\n+        if r_cand < r_cur && r_cand < best_resid {\n+            best_resid = r_cand;\n+            x0.copy_from_slice(&cand);\n+            if r_cand <= opts.tol {\n+                return finish(x0, chain, Method::PowerAitken, iters, opts.tol);\n+            }\n+        }\n+    }\n+    finish(x0, chain, Method::PowerAitken, iters, opts.tol)\n+}\n+\n+/// Gauss–Seidel / SOR on the balance equations.\n+///\n+/// Balance says `π_j · Σ_k a_jk = Σ_{i≠j} π_i a_ij`, i.e. every node's outflow\n+/// equals its inflow. Sweeping that in place (using already-updated components\n+/// within the sweep) propagates information across the whole graph in one pass\n+/// instead of one edge per pass, which is precisely what power iteration cannot\n+/// do on a long chain.\n+pub fn sor(chain: &RankChain, opts: SolveOptions, omega: f64) -> Solution {\n+    let n = chain.n;\n+    let cols = chain.columns();\n+    let mut pi = vec![1.0 / n as f64; n];\n+    let mut iters = 0usize;\n+\n+    for _ in 0..opts.max_iters {\n+        iters += 1;\n+        let mut delta = 0.0f64;\n+        for j in 0..n {\n+            let out = chain.row_sum[j];\n+            if out <= 0.0 {\n+                continue;\n+            }\n+            let inflow: f64 = cols[j].iter().map(|&(i, w)| pi[i] * w).sum();\n+            let gs = inflow / out;\n+            let updated = (1.0 - omega) * pi[j] + omega * gs;\n+            let updated = if updated.is_finite() && updated > 0.0 {\n+                updated\n+            } else {\n+                pi[j]\n+            };\n+            delta += (updated - pi[j]).abs();\n+            pi[j] = updated;\n+        }\n+        // The iteration is only defined up to scale; renormalize so `tol` means\n+        // the same thing here as it does for power iteration.\n+        let sum: f64 = pi.iter().sum();\n+        if sum.is_finite() && sum > 0.0 {\n+            for p in pi.iter_mut() {\n+                *p /= sum;\n+            }\n+        }\n+        if delta < opts.tol {\n+            break;\n+        }\n+    }\n+    finish(pi, chain, Method::Sor, iters, opts.tol)\n+}\n+\n+/// Jacobi-preconditioned BiCGSTAB on `Qᵀπ = 0` with row 0 replaced by `Σπ = 1`.\n+///\n+/// Included as the Krylov representative. GMRES/Arnoldi would need restart\n+/// bookkeeping and a stored basis for the same job; BiCGSTAB gives the same\n+/// short-recurrence answer for a nonsymmetric operator at fixed memory.\n+pub fn bicgstab(chain: &RankChain, opts: SolveOptions) -> Solution {\n+    let n = chain.n;\n+    let cols = chain.columns();\n+\n+    // A x = b where row 0 is Σx = 1 and row j>0 is (inflow - outflow) at j.\n+    let apply = |x: &[f64], out: &mut [f64]| {\n+        out[0] = x.iter().sum();\n+        for j in 1..n {\n+            let inflow: f64 = cols[j].iter().map(|&(i, w)| x[i] * w).sum();\n+            out[j] = inflow - x[j] * chain.row_sum[j];\n+        }\n+    };\n+    // Jacobi preconditioner: the diagonal of that operator.\n+    let diag: Vec<f64> = (0..n)\n+        .map(|j| {\n+            let d = if j == 0 { 1.0 } else { -chain.row_sum[j] };\n+            if d.abs() < 1e-300 {\n+                1.0\n+            } else {\n+                d\n+            }\n+        })\n+        .collect();\n+    let precond = |v: &[f64], out: &mut [f64]| {\n+        for i in 0..n {\n+            out[i] = v[i] / diag[i];\n+        }\n+    };\n+\n+    let mut b = vec![0.0f64; n];\n+    b[0] = 1.0;\n+\n+    let mut x = vec![1.0 / n as f64; n];\n+    let mut ax = vec![0.0f64; n];\n+    apply(&x, &mut ax);\n+    let mut r: Vec<f64> = (0..n).map(|i| b[i] - ax[i]).collect();\n+    let r_hat = r.clone();\n+\n+    let mut rho = 1.0f64;\n+    let mut alpha = 1.0f64;\n+    let mut omega = 1.0f64;\n+    let mut v = vec![0.0f64; n];\n+    let mut p = vec![0.0f64; n];\n+    let mut y = vec![0.0f64; n];\n+    let mut z = vec![0.0f64; n];\n+    let mut s = vec![0.0f64; n];\n+    let mut t = vec![0.0f64; n];\n+    let dot = |a: &[f64], b: &[f64]| -> f64 { (0..a.len()).map(|i| a[i] * b[i]).sum() };\n+    let bnorm = 1.0f64;\n+    let mut iters = 0usize;\n+\n+    for _ in 0..opts.max_iters {\n+        iters += 1;\n+        let rho_new = dot(&r_hat, &r);\n+        if rho_new.abs() < 1e-300 {\n+            break;\n+        }\n+        let beta = (rho_new / rho) * (alpha / omega);\n+        rho = rho_new;\n+        for i in 0..n {\n+            p[i] = r[i] + beta * (p[i] - omega * v[i]);\n+        }\n+        precond(&p, &mut y);\n+        apply(&y, &mut v);\n+        let denom = dot(&r_hat, &v);\n+        if denom.abs() < 1e-300 {\n+            break;\n+        }\n+        alpha = rho / denom;\n+        for i in 0..n {\n+            s[i] = r[i] - alpha * v[i];\n+        }\n+        if dot(&s, &s).sqrt() / bnorm < opts.tol * 1e-2 {\n+            for i in 0..n {\n+                x[i] += alpha * y[i];\n+            }\n+            break;\n+        }\n+        precond(&s, &mut z);\n+        apply(&z, &mut t);\n+        let tt = dot(&t, &t);\n+        if tt.abs() < 1e-300 {\n+            break;\n+        }\n+        omega = dot(&t, &s) / tt;\n+        for i in 0..n {\n+            x[i] += alpha * y[i] + omega * z[i];\n+            r[i] = s[i] - omega * t[i];\n+        }\n+        if dot(&r, &r).sqrt() / bnorm < opts.tol * 1e-2 {\n+            break;\n+        }\n+        if omega.abs() < 1e-300 {\n+            break;\n+        }\n+    }\n+\n+    // BiCGSTAB has no sign constraint; clamp before normalizing.\n+    for xi in x.iter_mut() {\n+        if !xi.is_finite() || *xi < 0.0 {\n+            *xi = 0.0;\n+        }\n+    }\n+    finish(x, chain, Method::BiCgStab, iters, opts.tol)\n+}\n+\n+// ---------------------------------------------------------------------------\n+// Direct solvers\n+// ---------------------------------------------------------------------------\n+\n+/// Dense LU with partial pivoting on the balance system.\n+///\n+/// `πQ = 0` is rank `n-1`, so equation 0 is replaced by the normalization\n+/// `Σπ = 1`, giving a nonsingular `n×n` system solved in `O(n³)`.\n+pub fn dense_lu(chain: &RankChain, tol: f64) -> Solution {\n+    let n = chain.n;\n+    // Row-major A, where A[j][i] is the coefficient of π_i in equation j.\n+    let mut a = vec![0.0f64; n * n];\n+    for v in a.iter_mut().take(n) {\n+        *v = 1.0; // equation 0: Σ π_i = 1\n+    }\n+    for i in 0..n {\n+        for &(j, w) in &chain.rows[i] {\n+            if j != 0 {\n+                a[j * n + i] += w; // inflow to j from i\n+            }\n+        }\n+        if i != 0 {\n+            a[i * n + i] -= chain.row_sum[i]; // outflow from i\n+        }\n+    }\n+    let mut b = vec![0.0f64; n];\n+    b[0] = 1.0;\n+\n+    // Gaussian elimination, partial pivoting.\n+    let mut perm: Vec<usize> = (0..n).collect();\n+    for k in 0..n {\n+        let mut piv = k;\n+        let mut best = a[perm[k] * n + k].abs();\n+        for r in (k + 1)..n {\n+            let v = a[perm[r] * n + k].abs();\n+            if v > best {\n+                best = v;\n+                piv = r;\n+            }\n+        }\n+        if best == 0.0 {\n+            continue;\n+        }\n+        perm.swap(k, piv);\n+        let pk = perm[k];\n+        let pivot = a[pk * n + k];\n+        for r in (k + 1)..n {\n+            let pr = perm[r];\n+            let f = a[pr * n + k] / pivot;\n+            if f == 0.0 {\n+                continue;\n+            }\n+            a[pr * n + k] = 0.0;\n+            for c in (k + 1)..n {\n+                a[pr * n + c] -= f * a[pk * n + c];\n+            }\n+            b[pr] -= f * b[pk];\n+        }\n+    }\n+    let mut x = vec![0.0f64; n];\n+    for k in (0..n).rev() {\n+        let pk = perm[k];\n+        let mut acc = b[pk];\n+        for c in (k + 1)..n {\n+            acc -= a[pk * n + c] * x[c];\n+        }\n+        let d = a[pk * n + k];\n+        x[k] = if d.abs() > 0.0 { acc / d } else { 0.0 };\n+    }\n+    for xi in x.iter_mut() {\n+        if !xi.is_finite() || *xi < 0.0 {\n+            *xi = 0.0;\n+        }\n+    }\n+    finish(x, chain, Method::DenseLu, 0, tol)\n+}\n+\n+/// Rescale threshold for GTH back-substitution.\n+///\n+/// Back-substitution accumulates `π` in *unnormalized* form, and on a strongly\n+/// ordered chain the values grow geometrically. Rescaling the partial vector\n+/// whenever it nears the top of the f64 range keeps every ratio exact (they are\n+/// all defined up to one global scalar) and avoids returning `inf`.\n+const GTH_RESCALE_ABOVE: f64 = 1e250;\n+\n+/// Dense Grassmann–Taksar–Heyman state reduction on a compact rate matrix.\n+///\n+/// Every operation is an addition of non-negative numbers or a division by a\n+/// positive sum: there is no subtraction anywhere, so there is no cancellation\n+/// and no pivoting is needed. That is why GTH is the reference implementation\n+/// here even though dense LU costs the same `O(n³)`.\n+///\n+/// Returns `ln π` shifted so the largest entry is 0. Finishing in logs matters\n+/// because the answer itself can span more than f64's range.\n+fn gth_dense_logs(n: usize, rows: &[Vec<(usize, f64)>]) -> Vec<f64> {\n+    if n == 0 {\n+        return vec![];\n+    }\n+    if n == 1 {\n+        return vec![0.0];\n+    }\n+    let mut a = vec![0.0f64; n * n];\n+    for (i, row) in rows.iter().enumerate().take(n) {\n+        for &(j, w) in row {\n+            if j < n && j != i {\n+                a[i * n + j] += w;\n+            }\n+        }\n+    }\n+\n+    for k in (1..n).rev() {\n+        let s: f64 = (0..k).map(|j| a[k * n + j]).sum();\n+        let reaches_survivors = s > 0.0;\n+        if !reaches_survivors {\n+            // State k cannot reach the surviving block, so it feeds no mass\n+            // back into it.\n+            for i in 0..k {\n+                a[i * n + k] = 0.0;\n+            }\n+            continue;\n+        }\n+        for i in 0..k {\n+            let f = a[i * n + k] / s;\n+            a[i * n + k] = f;\n+            if f == 0.0 {\n+                continue;\n+            }\n+            for j in 0..k {\n+                a[i * n + j] += f * a[k * n + j];\n+            }\n+        }\n+    }\n+\n+    // Cheap path: accumulate linearly, take logs at the end.\n+    let mut x = vec![0.0f64; n];\n+    x[0] = 1.0;\n+    let mut rescaled = false;\n+    for k in 1..n {\n+        let mut acc = 0.0f64;\n+        for i in 0..k {\n+            acc += x[i] * a[i * n + k];\n+        }\n+        x[k] = acc;\n+        if acc > GTH_RESCALE_ABOVE {\n+            let inv = 1.0 / acc;\n+            for v in x.iter_mut().take(k + 1) {\n+                *v *= inv;\n+            }\n+            rescaled = true;\n+        }\n+    }\n+    if !rescaled {\n+        return shift_logs(x.iter().map(|&v| v.ln()).collect());\n+    }\n+\n+    // The span overflowed a single f64 vector, so redo the accumulation in log\n+    // space. The elimination coefficients are unchanged; only the sum differs.\n+    let mut logx = vec![f64::NEG_INFINITY; n];\n+    logx[0] = 0.0;\n+    let mut terms: Vec<f64> = Vec::with_capacity(n);\n+    for k in 1..n {\n+        terms.clear();\n+        for i in 0..k {\n+            let c = a[i * n + k];\n+            if c > 0.0 && logx[i].is_finite() {\n+                terms.push(logx[i] + c.ln());\n+            }\n+        }\n+        logx[k] = log_sum_exp(&terms);\n+    }\n+    shift_logs(logx)\n+}\n+\n+/// Dense GTH over the whole chain. `O(n³)` time, `O(n²)` memory.\n+pub fn dense_gth(chain: &RankChain, tol: f64) -> Solution {\n+    let logs = gth_dense_logs(chain.n, &chain.rows);\n+    let pi = pi_from_logs(&logs);\n+    finish_with_logs(pi, logs, chain, Method::DenseGth, 0, tol)\n+}\n+\n+/// Shift log-scores so the maximum is exactly 0.\n+fn shift_logs(mut logs: Vec<f64>) -> Vec<f64> {\n+    let max = logs.iter().cloned().fold(f64::NEG_INFINITY, f64::max);\n+    if max.is_finite() {\n+        for l in logs.iter_mut() {\n+            *l -= max;\n+        }\n+    }\n+    logs\n+}\n+\n+/// Outcome of attempting the sparse direct path.\n+pub enum SparseGthOutcome {\n+    Solved(Solution),\n+    /// Elimination could not finish cheaply: states were left with degree above\n+    /// the guard, or the work budget ran out. Fall back to an iterative solver.\n+    TooDense {\n+        core: usize,\n+        work: u64,\n+    },\n+}\n+\n+/// `μ = m - n + c` over the symmetrized comparison graph: the number of\n+/// independent cycles, i.e. how far the graph is from being a forest. `O(n + m)`.\n+fn cyclomatic_number(chain: &RankChain) -> usize {\n+    let n = chain.n;\n+    let mut parent: Vec<usize> = (0..n).collect();\n+    fn find(parent: &mut [usize], mut x: usize) -> usize {\n+        while parent[x] != x {\n+            parent[x] = parent[parent[x]];\n+            x = parent[x];\n+        }\n+        x\n+    }\n+\n+    let mut edges = 0usize;\n+    let mut merges = 0usize;\n+    for i in 0..n {\n+        for &(j, _) in &chain.rows[i] {\n+            // Count each compared pair once: the reverse arc exists whenever\n+            // both sides got votes, and the pair is an edge either way.\n+            if j < i && chain.rows[j].binary_search_by(|p| p.0.cmp(&i)).is_ok() {\n+                continue;\n+            }\n+            edges += 1;\n+            let (ri, rj) = (find(&mut parent, i), find(&mut parent, j));\n+            if ri != rj {\n+                parent[ri] = rj;\n+                merges += 1;\n+            }\n+        }\n+    }\n+    // c = n - merges, so μ = m - n + c = m - merges.\n+    edges.saturating_sub(merges)\n+}\n+\n+/// Above this live degree, eliminating a state creates more fill than it\n+/// removes. Trees, chains, series-parallel graphs and most hand-built\n+/// ontologies reduce away completely under this guard; random graphs stall\n+/// almost immediately, which is exactly the signal to go iterative.\n+const MAX_ELIM_DEGREE: usize = 12;\n+\n+/// GTH state reduction over the sparse graph, eliminating minimum-degree states first.\n+///\n+/// Elimination on a graph is fill-bounded: removing a state connects its\n+/// surviving in-neighbours to its surviving out-neighbours. A chain — or any\n+/// tree — has an elimination order with *zero* fill, so this is `O(n)` on\n+/// exactly the topology where power iteration needs `Θ(n²)` sweeps. The two\n+/// failure modes are complementary: what stalls here (dense, well-connected\n+/// graphs) is what power iteration finishes in a few dozen sweeps.\n+///\n+/// States that cannot be eliminated cheaply form an irreducible *core*, handed\n+/// to dense GTH when it is small enough and reported as `TooDense` otherwise.\n+pub fn sparse_gth(chain: &RankChain, opts: SolveOptions, tol: f64) -> SparseGthOutcome {\n+    let n = chain.n;\n+    if n == 1 {\n+        return SparseGthOutcome::Solved(finish(vec![1.0], chain, Method::SparseGth, 0, tol));\n+    }\n+\n+    // Predict the core before doing any work. A graph whose every vertex has\n+    // degree ≥ 3 has at most `2μ - 2` vertices, where `μ = m - n + c` is the\n+    // cyclomatic number, so `2μ` bounds what elimination can leave behind. A\n+    // forest has `μ = 0` and reduces to nothing; a random graph has `μ ≈ m` and\n+    // reduces to almost nothing *but itself*, which is the case worth skipping\n+    // before paying for it.\n+    if n > opts.dense_core_max && 2 * cyclomatic_number(chain) > opts.dense_core_max {\n+        return SparseGthOutcome::TooDense { core: n, work: 0 };\n+    }\n+\n+    // Live out-rates, plus in-neighbour lists that may carry stale or duplicate\n+    // entries — they are compacted against `alive` when a state is popped.\n+    let mut out: Vec<Vec<(usize, f64)>> = chain.rows.clone();\n+    let mut inn: Vec<Vec<usize>> = vec![Vec::new(); n];\n+    for i in 0..n {\n+        for &(j, _) in &chain.rows[i] {\n+            inn[j].push(i);\n+        }\n+    }\n+\n+    let mut alive = vec![true; n];\n+    // Position of each column in the row currently being updated, -1 when\n+    // absent. Turns a row merge into O(len) instead of a map lookup per entry.\n+    let mut slot: Vec<i64> = vec![-1; n];\n+    let mut recovery: Vec<(usize, Vec<(usize, f64)>)> = Vec::with_capacity(n);\n+    let mut work: u64 = 0;\n+\n+    let mut heap: std::collections::BinaryHeap<std::cmp::Reverse<(usize, usize)>> =\n+        std::collections::BinaryHeap::new();\n+    for k in 0..n {\n+        heap.push(std::cmp::Reverse((out[k].len() + inn[k].len(), k)));\n+    }\n+\n+    let mut remaining = n;\n+    while remaining > 1 {\n+        let Some(std::cmp::Reverse((d, k))) = heap.pop() else {\n+            break;\n+        };\n+        if !alive[k] {\n+            continue;\n+        }\n+\n+        // Compact k's neighbour lists, then re-check the heap key against the\n+        // true degree: a stale key means another elimination has changed it.\n+        out[k].retain(|&(j, w)| alive[j] && w > 0.0);\n+        inn[k].sort_unstable();\n+        inn[k].dedup();\n+        inn[k].retain(|&i| alive[i] && i != k);\n+        let cur = out[k].len() + inn[k].len();\n+        if cur != d {\n+            heap.push(std::cmp::Reverse((cur, k)));\n+            continue;\n+        }\n+        if out[k].len().max(inn[k].len()) > MAX_ELIM_DEGREE {\n+            // The cheapest remaining state is already expensive, so every other\n+            // one is too. What is left is the core.\n+            break;\n+        }\n+\n+        let s: f64 = out[k].iter().map(|&(_, w)| w).sum();\n+        let in_list = std::mem::take(&mut inn[k]);\n+        let out_list = std::mem::take(&mut out[k]);\n+\n+        let mut coeffs: Vec<(usize, f64)> = Vec::with_capacity(in_list.len());\n+        for &i in &in_list {\n+            for (p, &(j, _)) in out[i].iter().enumerate() {\n+                slot[j] = p as i64;\n+            }\n+            work = work.saturating_add(out[i].len() as u64 + out_list.len() as u64);\n+\n+            let w_ik = if slot[k] >= 0 {\n+                let p = slot[k] as usize;\n+                let v = out[i][p].1;\n+                out[i][p].1 = 0.0;\n+                v\n+            } else {\n+                0.0\n+            };\n+            let f = if s > 0.0 { w_ik / s } else { 0.0 };\n+            if f > 0.0 {\n+                coeffs.push((i, f));\n+                for &(j, r) in &out_list {\n+                    if j == i {\n+                        continue;\n+                    }\n+                    if slot[j] >= 0 {\n+                        out[i][slot[j] as usize].1 += f * r;\n+                    } else {\n+                        slot[j] = out[i].len() as i64;\n+                        out[i].push((j, f * r));\n+                        inn[j].push(i);\n+                    }\n+                }\n+            }\n+            for &(j, _) in out[i].iter() {\n+                slot[j] = -1;\n+            }\n+            slot[k] = -1;\n+            out[i].retain(|&(j, w)| j != k && w > 0.0);\n+        }\n+\n+        alive[k] = false;\n+        remaining -= 1;\n+        recovery.push((k, coeffs));\n+\n+        for &i in &in_list {\n+            if alive[i] {\n+                heap.push(std::cmp::Reverse((out[i].len() + inn[i].len(), i)));\n+            }\n+        }\n+        for &(j, _) in &out_list {\n+            if alive[j] {\n+                heap.push(std::cmp::Reverse((out[j].len() + inn[j].len(), j)));\n+            }\n+        }\n+\n+        if work > opts.direct_work_budget {\n+            return SparseGthOutcome::TooDense {\n+                core: (0..n).filter(|&i| alive[i]).count(),\n+                work,\n+            };\n+        }\n+    }\n+\n+    let core: Vec<usize> = (0..n).filter(|&i| alive[i]).collect();\n+    if core.len() > opts.dense_core_max {\n+        return SparseGthOutcome::TooDense {\n+            core: core.len(),\n+            work,\n+        };\n+    }\n+\n+    // Solve the irreducible core densely, then unwind the eliminated states.\n+    let mut compact = vec![usize::MAX; n];\n+    for (c, &g) in core.iter().enumerate() {\n+        compact[g] = c;\n+    }\n+    let core_rows: Vec<Vec<(usize, f64)>> = core\n+        .iter()\n+        .map(|&g| {\n+            let mut row: Vec<(usize, f64)> = out[g]\n+                .iter()\n+                .filter(|&&(j, w)| alive[j] && w > 0.0)\n+                .map(|&(j, w)| (compact[j], w))\n+                .collect();\n+            row.sort_unstable_by(|a, b| a.0.cmp(&b.0));\n+            row\n+        })\n+        .collect();\n+\n+    let core_logs = gth_dense_logs(core.len(), &core_rows);\n+    let mut logx = vec![f64::NEG_INFINITY; n];\n+    for (c, &g) in core.iter().enumerate() {\n+        logx[g] = core_logs[c];\n+    }\n+    let mut terms: Vec<f64> = Vec::new();\n+    for (k, coeffs) in recovery.iter().rev() {\n+        terms.clear();\n+        for &(i, f) in coeffs {\n+            if f > 0.0 && logx[i].is_finite() {\n+                terms.push(logx[i] + f.ln());\n+            }\n+        }\n+        logx[*k] = log_sum_exp(&terms);\n+    }\n+\n+    let shifted = shift_logs(logx);\n+    let pi = pi_from_logs(&shifted);\n+    SparseGthOutcome::Solved(finish_with_logs(\n+        pi,\n+        shifted,\n+        chain,\n+        Method::SparseGth,\n+        0,\n+        tol,\n+    ))\n+}\n+\n+// ---------------------------------------------------------------------------\n+// Recommended strategy\n+// ---------------------------------------------------------------------------\n+\n+/// Solve for π using the hybrid strategy: sparse direct first, iterative only\n+/// when the direct path is genuinely too expensive.\n+///\n+/// This never returns an unconverged answer without saying so — check\n+/// [`Solution::converged`].\n+pub fn solve(chain: &RankChain, opts: SolveOptions) -> Solution {\n+    if chain.n == 0 {\n+        return Solution {\n+            pi: vec![],\n+            log_pi: vec![],\n+            method: Method::SparseGth,\n+            iterations: 0,\n+            residual: 0.0,\n+            converged: true,\n+            underflowed: false,\n+        };\n+    }\n+    if chain.n == 1 {\n+        return Solution {\n+            pi: vec![1.0],\n+            log_pi: vec![0.0],\n+            method: Method::SparseGth,\n+            iterations: 0,\n+            residual: 0.0,\n+            converged: true,\n+            underflowed: false,\n+        };\n+    }\n+    if chain.nnz() == 0 {\n+        let n = chain.n;\n+        return Solution {\n+            pi: vec![1.0 / n as f64; n],\n+            log_pi: vec![0.0; n],\n+            method: Method::SparseGth,\n+            iterations: 0,\n+            residual: 0.0,\n+            converged: true,\n+            underflowed: false,\n+        };\n+    }\n+\n+    // A chain over several disconnected components has no unique stationary\n+    // distribution: any split of mass between them is stationary. Power\n+    // iteration from a uniform start picks the one where each component holds\n+    // its share of the nodes, so reproduce that explicitly instead of letting a\n+    // direct solver pick one component and zero the rest.\n+    let comps = weak_components(chain);\n+    if comps.len() > 1 {\n+        return solve_by_component(chain, &comps, opts);\n+    }\n+    solve_connected(chain, opts)\n+}\n+\n+/// Weakly connected components of the comparison graph, each sorted ascending.\n+fn weak_components(chain: &RankChain) -> Vec<Vec<usize>> {\n+    let n = chain.n;\n+    let mut adj: Vec<Vec<usize>> = vec![Vec::new(); n];\n+    for i in 0..n {\n+        for &(j, _) in &chain.rows[i] {\n+            adj[i].push(j);\n+            adj[j].push(i);\n+        }\n+    }\n+    let mut seen = vec![false; n];\n+    let mut comps = Vec::new();\n+    let mut stack = Vec::new();\n+    for start in 0..n {\n+        if seen[start] {\n+            continue;\n+        }\n+        seen[start] = true;\n+        stack.push(start);\n+        let mut comp = Vec::new();\n+        while let Some(x) = stack.pop() {\n+            comp.push(x);\n+            for &y in &adj[x] {\n+                if !seen[y] {\n+                    seen[y] = true;\n+                    stack.push(y);\n+                }\n+            }\n+        }\n+        comp.sort_unstable();\n+        comps.push(comp);\n+    }\n+    comps\n+}\n+\n+fn solve_by_component(chain: &RankChain, comps: &[Vec<usize>], opts: SolveOptions) -> Solution {\n+    let n = chain.n;\n+    let mut log_pi = vec![f64::NEG_INFINITY; n];\n+    let mut converged = true;\n+    let mut iterations = 0usize;\n+    let mut method = Method::SparseGth;\n+\n+    for comp in comps {\n+        let mut compact = vec![usize::MAX; n];\n+        for (c, &g) in comp.iter().enumerate() {\n+            compact[g] = c;\n+        }\n+        let rows: Vec<Vec<(usize, f64)>> = comp\n+            .iter()\n+            .map(|&g| {\n+                chain.rows[g]\n+                    .iter()\n+                    .filter(|&&(j, _)| compact[j] != usize::MAX)\n+                    .map(|&(j, w)| (compact[j], w))\n+                    .collect()\n+            })\n+            .collect();\n+        let row_sum = rows\n+            .iter()\n+            .map(|r| r.iter().map(|&(_, w)| w).sum())\n+            .collect();\n+        let sub = RankChain {\n+            n: comp.len(),\n+            rows,\n+            row_sum,\n+            d_max: chain.d_max,\n+        };\n+        let sol = solve_connected(&sub, opts);\n+        converged &= sol.converged;\n+        iterations = iterations.max(sol.iterations);\n+        if !sol.method.is_direct() {\n+            method = sol.method;\n+        }\n+        // Each component keeps the share of the mass its node count started\n+        // with, spread internally by its own stationary distribution.\n+        let share = (comp.len() as f64 / n as f64).ln() - log_sum_exp(&sol.log_pi);\n+        for (c, &g) in comp.iter().enumerate() {\n+            log_pi[g] = sol.log_pi[c] + share;\n+        }\n+    }\n+\n+    let shifted = shift_logs(log_pi);\n+    let pi = pi_from_logs(&shifted);\n+    let mut out = finish_with_logs(pi, shifted, chain, method, iterations, opts.tol);\n+    // One component missing tolerance condemns the whole vector, even if the\n+    // combined residual happens to look small next to the dominant component.\n+    out.converged &= converged;\n+    out\n+}\n+\n+fn solve_connected(chain: &RankChain, opts: SolveOptions) -> Solution {\n+    match sparse_gth(chain, opts, opts.tol) {\n+        SparseGthOutcome::Solved(sol) if sol.converged => return sol,\n+        // A direct solve that misses tolerance means the conditioning is worse\n+        // than the elimination could handle; let the iterative path try rather\n+        // than hand back an answer nothing has checked.\n+        SparseGthOutcome::Solved(_) | SparseGthOutcome::TooDense { .. } => {}\n+    }\n+\n+    // Gauss–Seidel first: on every graph dense enough to reach this branch it\n+    // converged in tens of sweeps and beat power iteration on both time and\n+    // residual. Power iteration stays as the second opinion.\n+    let sweeps = sor(chain, opts, 1.0);\n+    if sweeps.converged {\n+        return sweeps;\n+    }\n+    let pow = power_aitken(chain, opts);\n+    if pow.converged || pow.residual < sweeps.residual {\n+        pow\n+    } else {\n+        sweeps\n+    }\n+}\n+\n+#[cfg(test)]\n+mod tests {\n+    use super::*;\n+\n+    const TOL: f64 = 1e-8;\n+\n+    fn opts() -> SolveOptions {\n+        SolveOptions::default()\n+    }\n+\n+    /// Build a chain the way production does: raw directed vote weights through\n+    /// `ranking::chain_from_edges`. `votes` are `(winner, loser, w_win, w_lose)`;\n+    /// Rank Centrality walks *toward* the winner, so the loser's row carries the\n+    /// weight (see `GroupState::apply_vote`).\n+    fn chain_from_votes(n: usize, votes: &[(usize, usize, f64, f64)]) -> RankChain {\n+        let mut raw: std::collections::HashMap<(usize, usize), f64> = Default::default();\n+        for &(w, l, ww, wl) in votes {\n+            *raw.entry((l, w)).or_insert(0.0) += ww;\n+            *raw.entry((w, l)).or_insert(0.0) += wl;\n+        }\n+        crate::ranking::chain_from_edges(n, raw.into_iter())\n+    }\n+\n+    fn path_graph(n: usize, ratio: f64) -> RankChain {\n+        let votes: Vec<(usize, usize, f64, f64)> =\n+            (0..n - 1).map(|i| (i, i + 1, ratio, 1.0)).collect();\n+        chain_from_votes(n, &votes)\n+    }\n+\n+    fn star_graph(n: usize) -> RankChain {\n+        let votes: Vec<(usize, usize, f64, f64)> =\n+            (1..n).map(|i| (0, i, 2.0 + (i % 5) as f64, 1.0)).collect();\n+        chain_from_votes(n, &votes)\n+    }\n+\n+    fn clique_graph(n: usize) -> RankChain {\n+        let mut votes = Vec::new();\n+        let mut k = 0;\n+        for i in 0..n {\n+            for j in (i + 1)..n {\n+                votes.push((i, j, 2.0 + (k % 5) as f64, 1.0));\n+                k += 1;\n+            }\n+        }\n+        chain_from_votes(n, &votes)\n+    }\n+\n+    fn sparse_graph(n: usize, degree: usize) -> RankChain {\n+        let mut state = 0x00C0_FFEEu64;\n+        let mut next = || {\n+            state ^= state >> 12;\n+            state ^= state << 25;\n+            state ^= state >> 27;\n+            state.wrapping_mul(0x2545_F491_4F6C_DD1D)\n+        };\n+        let mut votes: Vec<(usize, usize, f64, f64)> =\n+            (0..n - 1).map(|i| (i, i + 1, 2.0, 1.0)).collect();\n+        for k in 0..(degree.saturating_sub(2) * n / 2) {\n+            let a = (next() % n as u64) as usize;\n+            let b = (next() % n as u64) as usize;\n+            if a != b {\n+                votes.push((a, b, 2.0 + (k % 5) as f64, 1.0));\n+            }\n+        }\n+        chain_from_votes(n, &votes)\n+    }\n+\n+    fn max_rel_error(got: &[f64], want: &[f64]) -> f64 {\n+        (0..want.len())\n+            .filter(|&i| want[i] > 1e-280)\n+            .map(|i| (got[i] - want[i]).abs() / want[i])\n+            .fold(0.0f64, f64::max)\n+    }\n+\n+    fn force_sparse(chain: &RankChain) -> Solution {\n+        let unlimited = SolveOptions {\n+            direct_work_budget: u64::MAX,\n+            dense_core_max: usize::MAX,\n+            ..opts()\n+        };\n+        match sparse_gth(chain, unlimited, TOL) {\n+            SparseGthOutcome::Solved(s) => s,\n+            SparseGthOutcome::TooDense { .. } => panic!(\"uncapped elimination gave up\"),\n+        }\n+    }\n+\n+    // -----------------------------------------------------------------\n+    // Closed forms: these pin the answer without trusting any solver.\n+    // -----------------------------------------------------------------\n+\n+    /// A tree is reversible, so detailed balance fixes π edge by edge:\n+    /// `π_i / π_j = a_ji / a_ij`. For a path at a constant ratio that is a\n+    /// geometric sequence, known exactly.\n+    #[test]\n+    fn path_matches_detailed_balance_closed_form() {\n+        for n in [2usize, 3, 8, 64, 300] {\n+            let chain = path_graph(n, 3.0);\n+            let mut want: Vec<f64> = (0..n).map(|i| 3f64.powi(-(i as i32))).collect();\n+            let sum: f64 = want.iter().sum();\n+            for w in want.iter_mut() {\n+                *w /= sum;\n+            }\n+            for sol in [\n+                force_sparse(&chain),\n+                dense_gth(&chain, TOL),\n+                dense_lu(&chain, TOL),\n+                solve(&chain, opts()),\n+            ] {\n+                assert!(\n+                    max_rel_error(&sol.pi, &want) < 1e-12,\n+                    \"{} wrong on path n={n}: {:?}\",\n+                    sol.method.label(),\n+                    &sol.pi[..4.min(n)]\n+                );\n+                assert!(sol.converged, \"{} not converged\", sol.method.label());\n+            }\n+        }\n+    }\n+\n+    /// The 3-cycle a>b>c>a is symmetric under rotation, so π must be uniform.\n+    /// This is the `test/fixtures/ranking/cycle.sorter` anchor.\n+    #[test]\n+    fn three_cycle_is_exactly_uniform() {\n+        let chain = chain_from_votes(3, &[(0, 1, 2.0, 1.0), (1, 2, 2.0, 1.0), (2, 0, 2.0, 1.0)]);\n+        let sol = solve(&chain, opts());\n+        for (i, &p) in sol.pi.iter().enumerate() {\n+            assert!(\n+                (p - 1.0 / 3.0).abs() < 1e-14,\n+                \"node {i} got {p}, expected 1/3\"\n+            );\n+        }\n+    }\n+\n+    /// Issue #146: a pure forward star at 2:1. The hub beat both spokes, so it\n+    /// must rank first — and by detailed balance on this tree, exactly 2:1:1.\n+    #[test]\n+    fn star_regression_issue_146() {\n+        let chain = chain_from_votes(3, &[(0, 1, 2.0, 1.0), (0, 2, 2.0, 1.0)]);\n+        let sol = solve(&chain, opts());\n+        assert!((sol.pi[0] - 0.5).abs() < 1e-14, \"hub: {}\", sol.pi[0]);\n+        assert!((sol.pi[1] - 0.25).abs() < 1e-14, \"spoke a: {}\", sol.pi[1]);\n+        assert!((sol.pi[2] - 0.25).abs() < 1e-14, \"spoke b: {}\", sol.pi[2]);\n+    }\n+\n+    /// A single comparison at ratio r:1 must give exactly r:1.\n+    #[test]\n+    fn single_pair_is_the_ratio() {\n+        for r in [1.0f64, 2.0, 3.0, 7.0] {\n+            let chain = chain_from_votes(2, &[(0, 1, r, 1.0)]);\n+            let sol = solve(&chain, opts());\n+            let want = r / (r + 1.0);\n+            assert!(\n+                (sol.pi[0] - want).abs() < 1e-15,\n+                \"ratio {r}: got {:?}\",\n+                sol.pi\n+            );\n+        }\n+    }\n+\n+    // -----------------------------------------------------------------\n+    // Cross-method equivalence\n+    // -----------------------------------------------------------------\n+\n+    #[test]\n+    fn all_methods_agree_on_well_conditioned_graphs() {\n+        let cases: Vec<(&str, RankChain)> = vec![\n+            (\"star-64\", star_graph(64)),\n+            (\"clique-40\", clique_graph(40)),\n+            (\"sparse-d6-120\", sparse_graph(120, 6)),\n+            (\"path-40\", path_graph(40, 2.0)),\n+        ];\n+        let tight = SolveOptions {\n+            tol: 1e-14,\n+            max_iters: 500_000,\n+            ..opts()\n+        };\n+        for (label, chain) in cases {\n+            let reference = dense_gth(&chain, TOL).pi;\n+            for sol in [\n+                force_sparse(&chain),\n+                dense_lu(&chain, TOL),\n+                sor(&chain, tight, 1.0),\n+                power(&chain, tight),\n+                bicgstab(&chain, tight),\n+                solve(&chain, opts()),\n+            ] {\n+                // On a path the scores span 2^39, so an iterative method that\n+                // has driven the *absolute* residual to 1e-14 is still far off\n+                // in relative terms at the tail. Only the direct methods get\n+                // full relative precision there; that gap is the whole point.\n+                let rel = max_rel_error(&sol.pi, &reference);\n+                let bound = if label == \"path-40\" && !sol.method.is_direct() {\n+                    1e-2\n+                } else {\n+                    1e-9\n+                };\n+                assert!(\n+                    rel < bound,\n+                    \"{label}: {} disagrees with dense GTH by {rel:e}\",\n+                    sol.method.label()\n+                );\n+            }\n+        }\n+    }\n+\n+    /// Every solver's own residual must match what an independent recomputation\n+    /// says, so `Solution::residual` can be trusted as the acceptance criterion.\n+    #[test]\n+    fn reported_residual_matches_recomputation() {\n+        let chain = sparse_graph(200, 6);\n+        for sol in [\n+            solve(&chain, opts()),\n+            power(&chain, opts()),\n+            sor(&chain, opts(), 1.0),\n+            dense_gth(&chain, TOL),\n+        ] {\n+            let recomputed = chain.residual(&sol.pi);\n+            assert!(\n+                (recomputed - sol.residual).abs() <= 1e-18 + recomputed * 1e-9,\n+                \"{}: reported {:e} vs recomputed {recomputed:e}\",\n+                sol.method.label(),\n+                sol.residual\n+            );\n+        }\n+    }\n+\n+    // -----------------------------------------------------------------\n+    // No silent non-convergence\n+    // -----------------------------------------------------------------\n+\n+    #[test]\n+    fn power_iteration_reports_its_own_failure() {\n+        let chain = path_graph(1500, 2.0);\n+        let capped = SolveOptions {\n+            max_iters: 50,\n+            ..opts()\n+        };\n+        let sol = power(&chain, capped);\n+        assert!(!sol.converged, \"50 sweeps cannot solve a 1500-node path\");\n+        assert_eq!(sol.iterations, 50, \"should have used the whole budget\");\n+        assert!(sol.residual > TOL);\n+    }\n+\n+    /// The same input the capped power iteration fails on must come back solved\n+    /// through the hybrid, with the flag set honestly.\n+    #[test]\n+    fn hybrid_converges_where_power_iteration_cannot() {\n+        for n in [1500usize, 4000] {\n+            let chain = path_graph(n, 2.0);\n+            let pow = power(&chain, opts());\n+            assert!(\n+                !pow.converged,\n+                \"n={n}: power iteration was expected to hit the 10k cap\"\n+            );\n+\n+            let sol = solve(&chain, opts());\n+            assert!(sol.converged, \"n={n}: hybrid failed\");\n+            assert_eq!(sol.method, Method::SparseGth);\n+            assert!(sol.residual < 1e-12, \"n={n}: residual {:e}\", sol.residual);\n+        }\n+    }\n+\n+    /// On a path the true ranking is the construction order. `pi` runs out of\n+    /// exponent range past ~1000 nodes, but `log_pi` must stay strictly ordered\n+    /// all the way down — that is what makes it the right sort key.\n+    #[test]\n+    fn long_path_ranking_is_strictly_ordered_in_log_space() {\n+        let n = 3000;\n+        let chain = path_graph(n, 2.0);\n+        let sol = solve(&chain, opts());\n+        assert!(sol.converged);\n+        assert!(\n+            sol.underflowed,\n+            \"a 3000-node path at 2:1 spans 10^903 and must report underflow\"\n+        );\n+        for i in 0..n - 1 {\n+            assert!(\n+                sol.log_pi[i] > sol.log_pi[i + 1],\n+                \"log scores tie or invert at {i}: {} vs {}\",\n+                sol.log_pi[i],\n+                sol.log_pi[i + 1]\n+            );\n+        }\n+        // Every adjacent step is exactly ln 2 apart (detailed balance).\n+        for i in 0..n - 1 {\n+            let step = sol.log_pi[i] - sol.log_pi[i + 1];\n+            assert!(\n+                (step - 2f64.ln()).abs() < 1e-9,\n+                \"step at {i} is {step}, expected ln 2\"\n+            );\n+        }\n+    }\n+\n+    // -----------------------------------------------------------------\n+    // Structure\n+    // -----------------------------------------------------------------\n+\n+    /// A chain over disconnected components has infinitely many stationary\n+    /// distributions. Power iteration from uniform picks the one where each\n+    /// component keeps its share of the nodes; the solver must pick the same\n+    /// one, or whole-group rankings would silently change meaning.\n+    #[test]\n+    fn disjoint_components_split_mass_by_size() {\n+        let chain = chain_from_votes(4, &[(0, 1, 3.0, 1.0), (2, 3, 2.0, 1.0)]);\n+        let sol = solve(&chain, opts());\n+        assert!(sol.pi.iter().all(|p| p.is_finite()));\n+        // Each pair holds 1/2 of the mass, split 3:1 and 2:1 internally.\n+        for (i, want) in [(0, 0.375), (1, 0.125), (2, 1.0 / 3.0), (3, 1.0 / 6.0)] {\n+            assert!(\n+                (sol.pi[i] - want).abs() < 1e-14,\n+                \"node {i}: got {} want {want}\",\n+                sol.pi[i]\n+            );\n+        }\n+\n+        // And that is what unbounded power iteration converges to.\n+        let reference = power(\n+            &chain,\n+            SolveOptions {\n+                tol: 1e-15,\n+                max_iters: 200_000,\n+                ..opts()\n+            },\n+        );\n+        assert!(max_rel_error(&sol.pi, &reference.pi) < 1e-9);\n+    }\n+\n+    /// A lone item with no comparisons keeps a `1/n` share, exactly as the\n+    /// uniform-start power iteration leaves it.\n+    #[test]\n+    fn isolated_nodes_keep_a_uniform_share() {\n+        let chain = chain_from_votes(3, &[(0, 1, 3.0, 1.0)]);\n+        let sol = solve(&chain, opts());\n+        assert!(\n+            (sol.pi[2] - 1.0 / 3.0).abs() < 1e-14,\n+            \"isolate: {}\",\n+            sol.pi[2]\n+        );\n+        assert!((sol.pi[0] - 0.5).abs() < 1e-14, \"winner: {}\", sol.pi[0]);\n+        assert!(\n+            (sol.pi[1] - 1.0 / 6.0).abs() < 1e-14,\n+            \"loser: {}\",\n+            sol.pi[1]\n+        );\n+    }\n+\n+    #[test]\n+    fn cyclomatic_number_counts_independent_cycles() {\n+        assert_eq!(\n+            cyclomatic_number(&path_graph(50, 2.0)),\n+            0,\n+            \"a path is a tree\"\n+        );\n+        assert_eq!(cyclomatic_number(&star_graph(50)), 0, \"a star is a tree\");\n+        let triangle = chain_from_votes(3, &[(0, 1, 2.0, 1.0), (1, 2, 2.0, 1.0), (2, 0, 2.0, 1.0)]);\n+        assert_eq!(cyclomatic_number(&triangle), 1);\n+        // K_n has n(n-1)/2 edges and n-1 spanning-tree edges.\n+        assert_eq!(cyclomatic_number(&clique_graph(6)), 15 - 5);\n+    }\n+\n+    #[test]\n+    fn hybrid_picks_direct_for_trees_and_iterative_for_dense_graphs() {\n+        assert_eq!(\n+            solve(&path_graph(4000, 2.0), opts()).method,\n+            Method::SparseGth\n+        );\n+        assert_eq!(solve(&clique_graph(60), opts()).method, Method::SparseGth);\n+        let big_sparse = solve(&sparse_graph(2000, 6), opts());\n+        assert_ne!(big_sparse.method, Method::SparseGth);\n+        assert!(big_sparse.converged);\n+    }\n+\n+    #[test]\n+    fn degenerate_inputs_do_not_panic() {\n+        let empty = RankChain::from_normalized(0, vec![], 0);\n+        assert!(solve(&empty, opts()).pi.is_empty());\n+\n+        let single = RankChain::from_normalized(1, vec![], 0);\n+        assert_eq!(solve(&single, opts()).pi, vec![1.0]);\n+\n+        // Nodes with no comparisons at all: uniform, and no division by zero.\n+        let isolated = RankChain::from_normalized(5, vec![], 0);\n+        let sol = solve(&isolated, opts());\n+        assert!(sol.pi.iter().all(|&p| (p - 0.2).abs() < 1e-15));\n+\n+        // A unanimous edge (the loser never scored) leaves one direction empty.\n+        let unanimous = chain_from_votes(2, &[(0, 1, 1.0, 0.0)]);\n+        let sol = solve(&unanimous, opts());\n+        assert!(sol.pi[0] > sol.pi[1], \"{:?}\", sol.pi);\n+        assert!(sol.pi.iter().all(|p| p.is_finite()));\n+    }\n+\n+    // -----------------------------------------------------------------\n+    // Determinism\n+    // -----------------------------------------------------------------\n+\n+    /// The same graph presented in a different edge order must give a\n+    /// bit-identical answer — `RankChain::from_normalized` sorts precisely so\n+    /// that `HashMap` iteration order cannot leak into the low bits.\n+    #[test]\n+    fn edge_input_order_does_not_change_the_result() {\n+        let n = 300;\n+        let base: Vec<((usize, usize), f64)> = {\n+            let mut v = Vec::new();\n+            for i in 0..n - 1 {\n+                v.push(((i + 1, i), 2.0 / 3.0));\n+                v.push(((i, i + 1), 1.0 / 3.0));\n+            }\n+            for i in 0..n / 3 {\n+                v.push(((i, (7 * i + 11) % n), 0.5));\n+                v.push((((7 * i + 11) % n, i), 0.5));\n+            }\n+            v\n+        };\n+        let reference = solve(&RankChain::from_normalized(n, base.clone(), 6), opts()).pi;\n+\n+        let mut state = 12345u64;\n+        for _ in 0..6 {\n+            let mut shuffled = base.clone();\n+            for i in (1..shuffled.len()).rev() {\n+                state ^= state >> 12;\n+                state ^= state << 25;\n+                state ^= state >> 27;\n+                let j = (state.wrapping_mul(0x2545_F491_4F6C_DD1D) % (i as u64 + 1)) as usize;\n+                shuffled.swap(i, j);\n+            }\n+            let got = solve(&RankChain::from_normalized(n, shuffled, 6), opts()).pi;\n+            assert_eq!(got, reference, \"shuffled edge order changed the scores\");\n+        }\n+    }\n+}\n","role":"user"}],"model":"~x-ai/grok-latest"}