jamestjsp · Oct 12, 2026 · Oct 12, 2026 · Oct 12, 2026
diff --git a/.ergo/journal.jsonl b/.ergo/journal.jsonl
@@ -0,0 +1 @@
+{"version":1,"task_id":"HHJNJ5","kind":"created","at":"2026-10-12T00:23:24.976736Z"}
diff --git a/.ergo/plans.jsonl b/.ergo/plans.jsonl
@@ -360,3 +360,4 @@
 {"type":"result","ts":"2026-05-11T00:15:38.967907Z","data":{"task_id":"HKA3XF","summary":"Final bb03ad/bb04ad report shows no row more than 1.5x slower than Fortran.","path":"benchmark_c_vs_fortran_report.md","sha256_at_attach":"be1164009b3b0151c7c9b584b48c0f98a9ec1d6cb36fac98f76a659735246aba","mtime_at_attach":"2026-05-11T00:15:20.784333433Z","git_commit_at_attach":"cca77d92b979c0d4d3a5fbe330a576cdeffcd834","ts":"2026-05-11T00:15:38.967907Z"}}
 {"type":"body","ts":"2026-05-11T00:15:38.97072Z","data":{"id":"HKA3XF","body":"## Goal\n- Run the final benchmark comparison and decide whether the C11 Lyapunov generator speed issue is closed.\n- Collect validation evidence for implementation review.\n\n## Result\n- Reinstalled the package after optimization and reran targeted BB03AD/BB04AD behavior tests.\n- Regenerated the combined `bb03ad bb04ad` C-vs-Fortran report after the BB04AD optimization and report-section work.\n- The generated `## Fortran Faster Rows` section shows only BB04AD ex4.4 rows at 1.11x, 1.07x, and 1.06x Fortran advantage, so no row remains more than 1.5x slower than Fortran and no follow-up task is required by this gate.\n\n## Validation Gates\n- `uv pip install .`\n- `.venv/bin/pytest tests/python/test_bb03ad.py tests/python/test_bb04ad.py -v` (30 passed)\n- `.venv/bin/meson compile -C build-bench-c`\n- `.venv/bin/python scripts/benchmark_c_vs_fortran.py --c-build-dir build-bench-c --output benchmark_c_vs_fortran_report.md -r bb03ad bb04ad`\n- `git status --short` checked; only `.ergo/plans.jsonl` and generated `benchmark_c_vs_fortran_report.md` remain outside code commits.\n","ts":"2026-05-11T00:15:38.97072Z"}}
 {"type":"state","ts":"2026-05-11T00:15:38.97072Z","data":{"id":"HKA3XF","state":"done","ts":"2026-05-11T00:15:38.97072Z"}}
+{"type":"transaction","version":1,"events":[{"type":"new_task","ts":"2026-10-12T00:23:24.976736Z","data":{"id":"HHJNJ5","uuid":"3810e153-8d5f-459f-8cd2-b66cd2969df2","epic_id":"","state":"todo","title":"AB13DD: FMA contraction diverges from Fortran reference","body":"Repro: clang -O2 (fp-contract=on) emits fmadd/fnmsub at src/AB/ab13dd.c lines ~213,405,428,454,581-587,635-639,984-985,1079-1080 (objdump -d -l on ab13dd.o); gfortran -O3 AB13DD.f.o has only 8 fused ops. Same class as AB13CD bug: C `gamma*gamma - x` contracted while reference hoists/rounds GAMMA**2 (AB13CD fixed on claude/intelligent-elgamal-797628 via separate `gamma2` statement; random sweep 3000 cases went 41 -\u003e 0 mismatches).\nNext: map each C fused site to the Fortran object's fused/unfused ops; split statements where reference rounds; random differential sweep of ab13dd vs Fortran (driver pattern: tests/python/test_ab13cd_fortran_reference.py).\nFiles: src/AB/ab13dd.c, SLICOT-Reference/src/AB13DD.f\n","created_at":"2026-10-12T00:23:24.976736Z"}}]}
diff --git a/python/wrappers/py_ab.c b/python/wrappers/py_ab.c
@@ -8227,7 +8227,7 @@ PyObject* py_ab13cd(PyObject* self, PyObject* args) {
     f64 *c_data = (f64*)PyArray_DATA(c_array);
     f64 *d_data = (f64*)PyArray_DATA(d_array);
 
-    f64 fpeak;
+    f64 fpeak = 0.0;
     i32 info;
 
     f64 hnorm = ab13cd(n, m, np, a_data, lda, b_data, ldb, c_data, ldc, d_data, ldd,

diff --git a/src/AB/ab13cd.c b/src/AB/ab13cd.c
@@ -246,15 +246,21 @@ f64 ab13cd(i32 n, i32 m, i32 np, const f64 *a, i32 lda, const f64 *b, i32 ldb,
     f64 gamma = zero;
     f64 gammau = zero;
 
-    while (iter < MAXIT) {
+    for (;;) {
         iter++;
+        if (iter > MAXIT) {
+            *info = 2;
+            return zero;
+        }
         gamma = (one + two * tol) * gammal;
+        // Separate statement blocks FMA contraction; gfortran -O3 reference hoists GAMMA**2 and rounds it.
+        f64 gamma2 = gamma * gamma;
 
         for (i32 j = 0; j < m; j++) {
             for (i32 i = 0; i <= j; i++) {
                 dwork[iw4_ + i + j * m] = -dwork[iw2 + i + j * m];
             }
-            dwork[iw4_ + j + j * m] = gamma * gamma - dwork[iw2 + j + j * m];
+            dwork[iw4_ + j + j * m] = gamma2 - dwork[iw2 + j + j * m];
         }
 
         SLC_DLACPY("Full", &m, &n, dwork, &m, &dwork[iw5_], &m);
@@ -276,7 +282,7 @@ f64 ab13cd(i32 n, i32 m, i32 np, const f64 *a, i32 lda, const f64 *b, i32 ldb,
             for (i32 i = 0; i <= j; i++) {
                 dwork[iw7 + i + j * np] = -dwork[iw3_ + i + j * np];
             }
-            dwork[iw7 + j + j * np] = gamma * gamma - dwork[iw3_ + j + j * np];
+            dwork[iw7 + j + j * np] = gamma2 - dwork[iw3_ + j + j * np];
         }
 
         SLC_DLACPY("Full", &np, &n, c, &ldc, &dwork[iw8], &np);
@@ -384,11 +390,6 @@ f64 ab13cd(i32 n, i32 m, i32 np, const f64 *a, i32 lda, const f64 *b, i32 ldb,
         }
     }
 
-    if (iter > MAXIT) {
-        *info = 2;
-        return zero;
-    }
-
     f64 hnorm = (gammal + gammau) / two;
 
     dwork[0] = (f64)lwamax;

diff --git a/tests/python/data/ab13cd_fortran_reference.npz b/tests/python/data/ab13cd_fortran_reference.npz
diff --git a/tests/python/test_ab13cd_fortran_reference.py b/tests/python/test_ab13cd_fortran_reference.py
@@ -29,16 +29,24 @@ def _ab13cd_driver(a, b, c, d) -> str:
 {fortran_matrix("a", a)}{fortran_matrix("b", b)}{fortran_matrix("c", c)}{fortran_matrix("d", d)}
   hnorm = AB13CD(n, m, p, a, n, b, n, c, p, d, p, 1.0d-9, iwork, dwork, &
        ldwork, cwork, lcwork, bwork, info)
-  print '(I0,1X,ES26.17E3)', info, hnorm
+  print '(I0,1X,ES26.17E3,1X,ES26.17E3)', info, hnorm, dwork(2)
 end program main
 """
 
 
 DATA = Path(__file__).parent / "data" / "ab13cd_fortran_reference.npz"
 
-# Pre-DLAPY2 C (hypot) returned exactly half the norm for half_norm_stall and
-# was off by 1.5% for lightly_damped_mode.
-CASES = ["half_norm_stall", "lightly_damped_mode"]
+# half_norm_stall, lightly_damped_mode: hypot (vs DLAPY2) changed the initial
+# OMEGA; half_norm_stall then hit the MAXIT exit (GAMMAU=0 -> hnorm/2).
+# maxit_info2: Fortran exhausts MAXIT -> INFO=2; C fell out of the loop with
+# GAMMAU=0 and INFO=0.
+# dlapy2_omega: hypot picked a different initial OMEGA.
+# gamma2_rounding, gamma2_dpotrf_info: C contracted GAMMA**2 - X into an FMA;
+# Fortran rounds GAMMA**2 first. Changed the Hamiltonian / DPOTRF outcome.
+CASES = [
+    "half_norm_stall", "lightly_damped_mode", "maxit_info2", "dlapy2_omega",
+    "gamma2_rounding", "gamma2_dpotrf_info",
+]
 
 
 @pytest.mark.parametrize("name", CASES)
@@ -47,14 +55,16 @@ def test_ab13cd_hnorm_matches_fortran(name, tmp_path):
 
     with np.load(DATA) as data:
         a, b, c, d = (data[f"{name}_{k}"] for k in "abcd")
-    info_f, hnorm_f = run_fortran_driver(_ab13cd_driver(a, b, c, d), tmp_path).split()
+    info_f, hnorm_f, fpeak_f = run_fortran_driver(_ab13cd_driver(a, b, c, d), tmp_path).split()
 
     n, m = b.shape
     p = c.shape[0]
-    hnorm, _, info = ab13cd(
+    hnorm, fpeak, info = ab13cd(
         n, m, p, np.asfortranarray(a), np.asfortranarray(b),
         np.asfortranarray(c), np.asfortranarray(d), 1e-9,
     )
 
-    assert info == int(info_f) == 0
-    np.testing.assert_allclose(hnorm, float(hnorm_f), rtol=1e-10, atol=0)
+    assert info == int(info_f)
+    if info == 0:
+        np.testing.assert_allclose(hnorm, float(hnorm_f), rtol=1e-10, atol=0)
+        np.testing.assert_allclose(fpeak, float(fpeak_f), rtol=1e-10, atol=0)
Original file line number	Diff line number	Diff line change
		@@ -0,0 +1 @@
		{"version":1,"task_id":"HHJNJ5","kind":"created","at":"2026-10-12T00:23:24.976736Z"}