Cross-platform golden-record drift in Lambda_A traced to adaptive quadrature + libm (env, not code)

[krystophny]

Summary

Cross-platform golden-record drift in Lambda_A (~4.6e-7 relative) is caused by adaptive quadrature in src/utils/integrate.f90 amplifying ULP-level libm differences. A three-version comparison confirms the drift is purely environment, not code, and reproduces the macOS-vs-Linux symptom already described in #88. The open question is whether to widen the tolerance (#88) or make the quadrature reproducible, which connects to the libneo/fortnum quadrature-replacement work.

What was measured

The golden-record corpus build (plasma/proj/golden/rabe) ran the landremanpaul2021_QH_lowres case three ways and compared with compare_golden.py:

• A: the committed historical test/golden/expected/rabe.nc (original binary, original build environment).
• B: the same source rebuilt at commit a8c3b21 with the current toolchain (gfortran 14.2.0).
• C: current main (1e5ef76), current toolchain.

Result:

• A vs B (compiler/environment): Lambda_A 4.6e-7, remainder 6.0e-8, Lambda_S 3.6e-9, lambda_LC_bB 2.8e-9, nu_star_crit 1.7e-10; s_tor, split_maxima, R exact.
• B vs C (code change a8c3b21 -> main): all nine variables bitwise identical.
• A vs C: identical to A vs B.

So all observed drift between the historical record and a current build is environment. No code change between a8c3b21 and main moves any output. The Lambda_A 4.6e-7 and the ~3e-9 floor match #88 exactly ("~4e-7 for Lambda_A through libm and adaptive-quadrature differences").

Mechanism

Lambda_A, nu_star_crit, and remainder are the variables that go through integrate_1d in src/utils/integrate.f90, which calls quadpack dqag (globally adaptive subdivision, jacobwilliams/quadpack via cmake/quadpack.cmake). The integrand is evaluated through libneo field/spline routines (pinned at libneo 752aae9c in cmake/libneo.cmake). Adaptive subdivision turns ULP-level libm and field-evaluation differences into a different sequence of panels, so the error lands two orders of magnitude above the 3e-9 floor of the non-quadrature outputs. The 3e-9 outputs are direct arithmetic; the 4.6e-7 outputs are the adaptive integrals.

This is the same class of non-reproducibility the libneo modularization and fortnum effort target (itpplasma/libneo#291): adaptive quadrature and special functions whose results depend on the math library and link configuration.

Why this is worth more than a tolerance bump

#88 widens the golden tolerance to 1e-6 and #86 reverted both that and the quadpack link change (#87). Widening to 1e-6 still catches algorithmic regressions (which start at 1e-3), but it hides whether Lambda_A's 4.6e-7 is benign rounding or an ill-conditioned integrand. The two-order gap between the quadrature outputs and the floor is the signal worth checking: a near-singular integrand plus a loose epsabs/epsrel in dqag would produce exactly this, and would also explain the macOS reproduction failures.

Concrete questions:

1. Is the Lambda_A integrand near-singular at the substitution boundaries (integrate_1d_substituted)? What epsabs/epsrel/limit does dqag use, and does tightening them shrink the cross-platform spread?
2. Does pinning vs bumping libneo (the field/spline side of the integrand) change Lambda_A independently of the quadrature? Test with LIBNEO_BRANCH against the pinned 752aae9c.
3. Would a fixed-order or reproducible-adaptive quadrature (a candidate for the fortnum quadrature library, Umbrella: replace avoidable numerical dependencies with AD-ready fortnum libneo#291/#288) make the golden record bit-reproducible across platforms, removing the need for the 1e-6 band-aid?

Suggested split

Use the golden-record corpus to separate the two regimes instead of one tolerance:

• Same-machine ref-vs-cur (the gr-runner build_ref vs build_cur model) stays strict (1e-10): catches real code changes to machine precision.
• Cross-platform / cross-compiler check runs at the looser tolerance, documented as the libm + adaptive-quadrature floor, until the quadrature is made deterministic.

References

• Drift first noted and band-aided: Widen golden-record tolerance for cross-platform runs #88 (cross-platform tolerance), Link quadpack publicly into utils_lib #87 (quadpack link), reverted in Revert direct-to-main quadpack link and golden tolerance changes #86; Build rabe with LLVM Fortran compilers #76 builds rabe with LLVM Fortran.
• Code: src/utils/integrate.f90 (dqag), cmake/quadpack.cmake, cmake/libneo.cmake (libneo 752aae9c), src/coefficients.f90, src/deviation.f90.
• Related: Umbrella: replace avoidable numerical dependencies with AD-ready fortnum libneo#291 (fortnum: replace avoidable numerical dependencies), #288 (quadrature).
• Found via the rabe golden-record corpus (gitlab.tugraz.at:plasma/proj/golden/rabe), three-version env-vs-code comparison.

[krystophny]

Follow-up with measured isolation (not just hypothesis):

FMA ruled out. Rebuilt current main with -ffp-contract=off: byte-identical to the default build, and the same drift vs the expected record. There is no -march=native/-mfma in cmake/compiler.cmake (-O2 -cpp), so no hardware FMA is emitted anyway. FP contraction is not the cause.

libneo source ruled out. Pinned at 752aae9c via cmake/libneo.cmake; FetchContent fetches that exact commit.

Code ruled out. Legacy rebuild at a8c3b21 is bitwise identical to current main.

Remaining cause: compiler-version + system libm. The expected record was produced on ubuntu-22.04 (gfortran 11, glibc 2.35); the drift appears on gfortran 14.2 / glibc 2.40 (and macOS ARM). ULP-level differences in libm and gfortran codegen, amplified by the adaptive dqag subdivision, give the 4.6e-7 in Lambda_A (vs the 3e-9 floor of the non-integral outputs).

Why CI stays green: make golden runs on push/PR, but always on ubuntu-22.04 with the same toolchain that made expected.nc, so it is byte-identical there. CI pins the platform implicitly via the runner image and never varies the compiler/libc, so it cannot observe cross-toolchain drift.

Options, weighed:

• nix flake (pin gfortran + glibc + libneo + quadpack): byte-identical builds and golden output on any Linux machine and CI; the right immediate fix for the "green CI, unreproducible record" gap. Does not bridge Linux/macOS, since nix-on-Darwin uses Apple libc/libm and x86/ARM FP differ at the ISA level.
• deterministic quadrature (fixed-node, the fortnum direction libneo#288): removes the adaptive amplification, dropping the drift from 4.6e-7 to ~1e-13 and likely bringing macOS under 1e-10. This is the root fix.
• two-tier test: strict same-toolchain ref-vs-cur at 1e-10 plus a documented cross-platform floor.

Recommendation: add a nix flake and run the golden under it on every PR (keep 1e-10); pursue deterministic quadrature via fortnum as the durable cure rather than widening the tolerance (#88).

[krystophny]

Design observation that explains why this surfaces here.

The golden record in rabe is a stored oracle: test/golden/expected/rabe.nc is a committed file, produced once on one toolchain, and compare_golden.py checks each fresh run against it at rtol=1e-10. Codes like SIMPLE use a recompute model instead: compare_golden_results.sh REF CUR builds two versions on the same machine in the same job (e.g. main vs a branch) and compares their fresh outputs, with no committed reference.

The two models answer different questions:

• Recompute (ref vs cur): "does this change move the answer relative to current main?" Both sides run through the identical compiler/libm on one runner, so the libm and adaptive-quadrature ULP noise is common to both and cancels. The comparison is bit-clean unless the code changed, on any platform.
• Stored oracle: "does today's build still match the bytes we committed?" This is inherently toolchain-bound. A run on a different gfortran/glibc differs from the stored bytes even when the code is unchanged, which is the 4.6e-7 in Lambda_A.

This is why CI stays green: make golden runs on ubuntu-22.04 with the toolchain that produced expected.nc, so it is byte-identical there; CI never varies the compiler/libc, so it cannot observe cross-toolchain drift. macOS ARM and newer-gfortran Linux fall outside that single pinned point.

Measured isolation, for completeness:

• FMA is not the cause: rebuilding current main with -ffp-contract=off is byte-identical to the default build and shows the same drift vs the stored record (no -march=native/-mfma in cmake/compiler.cmake, so no hardware FMA is emitted anyway).
• libneo source is not the cause: pinned at 752aae9c in cmake/libneo.cmake.
• Code is not the cause: a legacy rebuild at a8c3b21 is bitwise identical to current main.
• The remaining variable is the compiler version plus system libm (gfortran 11 / glibc 2.35 for the record, vs gfortran 14.2 / glibc 2.40 and macOS ARM), amplified by the adaptive dqag subdivision.

We are testing on our side whether pinning the whole toolchain with a nix flake (gfortran, glibc, cmake, netcdf; libneo stays at the existing 752aae9c FetchContent pin) gives byte-identical builds and golden output across machines. Results to follow.

[krystophny]

Automated triage pass (Claude), grounded in the repository at 794bdfd. Verify before acting.

Keeping this open, but narrowing it. Option 1 from the first comment (nix flake) has landed: commit 3628819 (PR #91) added flake.nix pinning nixos-22.11 (gfortran 11 + glibc 2.35, with the stated reason "newer compilers/libm drift past rtol=1e-10"), and the build/test steps of both .github/workflows/test.yml and .github/workflows/golden.yml now run through nix develop .#default. What remains of #90 is the root fix — the three concrete questions in the body — and the investigation below answers Question 1 from the code, which changes the shape of the remaining work.

Purpose

test/golden/expected/rabe.nc is the only regression oracle for the nine output variables (s_tor, Lambda_A, Lambda_B, nu_star_crit, Lambda_S, split_maxima, R, lambda_LC_bB, remainder). Today it is toolchain-bound: it passes only inside the nix shell, so a physicist building rabe on a laptop, or anyone bumping the flake pin, gets a red golden test with no way to tell rounding from a real regression. Deterministic quadrature removes the largest amplifier at source, bringing the quadrature outputs down to the same drift class as the direct-arithmetic outputs.

Current state

The entire adaptive path is 75 lines in src/utils/integrate.f90. The requested accuracy is the headline finding:

! src/utils/integrate.f90:27-31
real(quadkind), parameter :: abs_error_tol_quadkind = 0.0_dp
real(quadkind), parameter :: rel_error_tol_quadkind = 1.0e-6
integer, parameter :: order_key = 6
integer, parameter :: limit = 500
integer, parameter :: lenw = limit*4

With epsabs = 0, dqag terminates on a purely relative criterion at 1e-6, using the 61-point Gauss–Kronrod rule (key = 6). The golden comparison in test/golden/compare_golden.py:34-39 asserts rtol=1e-10, atol=0.0. The golden record demands four orders of magnitude more accuracy than the integrator was ever asked to deliver. The measured Lambda_A drift of 4.6e-7 is inside dqag's own requested tolerance — the integrator is doing exactly what it was told. (Minor, but fix it while you are there: 1.0e-6 on line 28 is a default-real literal, so the actual real(8) value is 9.99999997475e-07, not 1e-6.)

Question 1 of the body — "is the Lambda_A integrand near-singular at the substitution boundaries?" — is answerable without running anything. The four field-line integrals (src/fieldline/fieldline_integrals.f90:23-41) run from fieldline%phi_max(1) to phi_max(2) with

! src/fieldline/fieldline_integrands.f90:34
lambda_over_B_squared = sqrt(calc_lambda_squared(B, eta))/B**2.0_dp
! src/fieldline/fieldline_integrands.f90:52
lambda_squared = 1.0_dp - B*eta

and eta = fieldline%eta_b, which is set globally at src/fieldline/make_fieldline.f90:184-185:

flock%eta_b = (1.0_dp - 2.0_dp*machine_eps) &
              /get_global_B_max(flock%fieldlines)

machine_eps = epsilon(1.0_dp) (src/utils/constants.f90:6). nudge_maxima_inward (make_fieldline.f90:381-403) runs first (called at line 158) and recomputes fieldline%B_max at the nudged endpoints, so get_global_B_max (make_fieldline.f90:405-415, a maxval over all lines' B_max) returns a value that is actually attained at an interval endpoint. On the one field line carrying that global maximum, 1 - B*eta_b at the endpoint is therefore 2*machine_eps ~ 4.4e-16 — the nudge buys nothing there by construction.

Under the t**2 = x - a substitution (integrate_1d_substituted, integrate.f90:172-207, integrands at integrate.f90:209-227) the integrand near t = 0 behaves as 2t*sqrt(delta + c*t**4)/B**2 with delta ~ 4.4e-16. The transition sits at t* ~ (delta/c)**0.25 ~ 5e-4 for c = O(1) — a genuine boundary layer three to four decades inside the panel, whose resolution decision is made by dqag's error estimator on quantities that are pure roundoff. Answer to Question 1: yes, and specifically on the global-maximum field line. That is the amplifier.

There is a second candidate amplifier the issue does not separate out. refine_maximum (src/utils/find_extrema.f90:96-134) ends at line 133 with loc = x_sub(i_max) where i_max = maxloc(val_sub, dim=1), and its exit conditions cannot_resolve / cannot_resolve_at_edge (find_extrema.f90:138-165) use tol = epsilon(values(1))*2.0_dp as an absolute tolerance on the |B| samples. So the loop terminates precisely when the neighbouring |B| values are within 2*epsilon of each other, and the final maxloc is then decided at ULP level. A one-ULP libm change flips that maxloc, moving phi_max by the last bracket width; get_global_B_max may then pick a different line's maximum, shifting eta_b for every field line at once (max_n_fieldlines = 200 in test/golden/input/rabe.in). The two mechanisms must be separated before committing to a fix, because deterministic quadrature does not cure the second one.

What deterministic quadrature can and cannot buy

This bounds the whole design, so state it before writing code. The issue's own three-version measurement gives, for A vs B (environment only):

variable	drift	path
Lambda_A	4.6e-7	adaptive quadrature
remainder	6.0e-8	adaptive quadrature
Lambda_S	3.6e-9	direct arithmetic
lambda_LC_bB	2.8e-9	direct arithmetic
nu_star_crit	1.7e-10	adaptive quadrature

Lambda_S and lambda_LC_bB drift by ~3e-9 across toolchains without going through dqag at all. Removing the adaptive amplification therefore cannot make an unpinned make golden pass at rtol=1e-10: the floor is set by libm and field/spline evaluation, is ~30x above 1e-10, and is untouched by this work (it is what #92 addresses). The honest, checkable target for this issue is:

the quadrature outputs (Lambda_A, remainder, nu_star_crit) drop from the 1e-7 class into the same ~3e-9 class as the direct-arithmetic outputs, and the cross-toolchain golden job runs at a documented floor rather than at 1e-10.

The pinned job keeps rtol=1e-10. Do not gate this PR on an unpinned 1e-10 pass; it is not achievable by changing the integrator.

Reproduction / evidence

Baseline (this is the state the flake guarantees):

cd /path/to/rabe
nix develop .#default --command make
nix develop .#default --command make golden      # PASS, all nine variables

Drift, outside the pin (gfortran >= 13, glibc >= 2.38):

make            # system toolchain
make golden     # FAIL: Lambda_A rel ~4.6e-7 vs rtol=1e-10

Discriminating experiment — do this first, it is one afternoon and it decides the design. Both runs use test/golden/input/rabe.in, which names wout_LandremanPaul2021_QH_reactorScale_lowres_reference.nc (copied into the golden build dir from test/integration/vmec/input/ by test/golden/CMakeLists.txt:2):

1. In src/utils/integrate.f90, print the subdivision count last (already returned by dqag, integrate.f90:66) on every integrate_1d call. Run under the nix pin and under the system toolchain. If the last sequences differ, the panel structure is the amplifier.
2. Independently, print fieldline%phi_max(1:2), fieldline%B_max, and flock%eta_b after make_fieldline.f90:186 for all field lines. If eta_b or any phi_max differs between toolchains, the extremum finder is an amplifier in its own right, independently of the quadrature.
3. Tighten rel_error_tol_quadkind to 1.0e-12_dp and raise limit to 2000, rebuild both ways, rerun make golden. Expected if the quadrature dominates: Lambda_A's cross-toolchain spread collapses from 4.6e-7 toward the ~3e-9 floor of Lambda_S/lambda_LC_bB. If it does not move, mechanism 2 dominates.

Step 3 is also the cheap partial answer to Question 1 of the body, and is worth reporting here even if the full fix is deferred.

Behaviour (BDD)

Scenario: Quadrature outputs reach the direct-arithmetic drift floor
  Given rabe built under `nix develop .#default` and again with a system gfortran >= 14 / glibc >= 2.39
  When both builds run test/golden/input/rabe.in and the two rabe.nc are compared
  Then the relative difference in Lambda_A is below 1e-8
  And the relative difference in remainder is below 1e-8
  And no variable exceeds 1e-8

Scenario: Quadrature is reproducible for a fixed integrand
  Given test/unit/test_quadrature_gauss.f90 integrating 2*t*sqrt(4.4e-16 + t**4) on [0, 1]
  When the same source is built with -O0 and with -O2
  Then both results agree with each other to rtol < 1e-13
  And both agree with the reference value quoted in the test to rtol < 1e-13

Scenario: Deterministic rule reproduces the adaptive answer on the physics case
  Given the QH lowres reference wout and the deterministic integrator enabled
  When rabe.x runs inside `nix develop .#default`
  Then Lambda_A, Lambda_B, nu_star_crit and remainder agree with the current
       test/golden/expected/rabe.nc to rtol < 1e-6
  And the record is regenerated with `make golden_update` and the per-variable
      shift is tabulated in the PR body

Scenario: Existing unit tests still pass on the new rule
  Given the six analytic cases in test/unit/test_integrate_1d_substituted.f90
  When `ctest --test-dir build/test -R TestIntegrate1dSubstituted` runs
  Then all six pass at the file's existing reltol = 1.0e-8_dp with abstol = 0.0_dp
       (test/unit/test_integrate_1d_substituted.f90:8)

Implementation plan

1. src/utils/integrate.f90 — split the file. It holds two modules (integrate at line 1, integrate_substituted at line 158) and is the only file under src/ that does; CLAUDE.md:97 asks for "one focused responsibility per module" but does not mandate one module per file, so treat this as optional tidying, not a rule violation. If you do it, move integrate_substituted into a new src/utils/integrate_substituted.f90 unchanged and add it to the add_library(utils_lib STATIC ...) list in src/utils/CMakeLists.txt. No behaviour change; keep it as its own commit so the numerical diff below is readable.
2. New src/utils/quadrature_gauss.f90 — module quadrature_gauss, public gauss_legendre_nodes(n, x, w) and integrate_gauss_fixed(f, a, b, n_panels, n_order, result). Hard-code the n = 32 abscissae/weights to 17 significant digits rather than computing them by Golub–Welsch: it removes the eigensolver and makes the rule literally identical on every platform. (cmake/libneo.cmake:11 forces LIBNEO_BUILD_NUMERICS OFF so that "BLAS/LAPACK and fortnum are not required here", so an eigensolver would mean adding a dependency rabe currently avoids.) Sum panel contributions in a fixed order with compensated (Kahan–Babuška) summation so the result does not depend on FMA contraction or vectorisation width.
3. src/utils/integrate.f90 — add subroutine integrate_1d_deterministic(f, a, b, result) with the identical signature to integrate_1d (integrate.f90:22-25), delegating to integrate_gauss_fixed with defaults n_panels = 4, n_order = 32. Keep integrate_1d/dqag in place, selected by a module-level logical, save :: use_adaptive_quadrature = .false. with a public setter, so the adaptive path stays available for cross-checking. Do not change the integrand module-pointer idiom (integrate.f90:18, doc note at 98-112) — the thread-safety note there still applies.
4. integrate_substituted — switch the two integrate_1d calls at integrate.f90:199-200 to the deterministic entry point. The t**2 substitution stays: it is what makes a fixed-order rule legitimate, because it turns the endpoint sqrt behaviour into a polynomial.
5. New test/unit/test_integrate_deterministic.f90 — reuse the six analytic integrands already in test/unit/test_integrate_1d_substituted.f90 (checks at lines 32, 44, 56, 68, 80, 93) against the new rule at reltol = 1.0e-10_dp, plus one boundary-layer case 2*t*sqrt(4.4e-16 + t**4) on [0,1] with the reference value computed by mpmath and quoted in a comment. Register in test/unit/CMakeLists.txt following the four-line pattern at lines 16-19 (add_executable, target_link_libraries(... utils_lib), add_test, set_tests_properties(... LABELS "quick")).
6. New test/unit/test_quadrature_gauss.f90 — assert exactness of the n = 32 rule on x**k for k = 0 .. 63 (degree 2n-1) to rtol 1e-14, and assert node/weight symmetry to bitwise equality. Same CMake registration.
7. test/golden/compare_golden.py — add an optional third argument or --rtol/--atol flag. Today the script takes exactly two positional paths and hardcodes rtol=1e-10, atol=0.0 at lines 34-39, with no way to compare at another tolerance; both the record regeneration in step 8 and the cross-toolchain job in step 9 need one. Default must stay 1e-10 so GoldenRecordCompare (test/golden/CMakeLists.txt:13-17) is unaffected.
8. Regenerate the record. Run nix develop .#default --command make golden_run, compare the fresh build/test/golden/rabe.nc against test/golden/expected/rabe.nc using the relaxed tolerance added in step 7, record the shift per variable in the PR description, then make golden_update (which copies build/test/golden/rabe.nc over the expected file, Makefile:48-49). The shift is an algorithmic change and must be justified, not silently absorbed.
9. Add a cross-toolchain golden job. In .github/workflows/golden.yml, add a second job golden-unpinned running make golden without nix develop, at the documented floor (start at rtol=1e-8) rather than at 1e-10 — see "What deterministic quadrature can and cannot buy" above. It must run on a newer image than ubuntu-22.04: benchmark.yml:14,30-32 shows the apt-gfortran installation pattern to copy, but it pins ubuntu-22.04, whose gfortran 11 is the same version as the flake, so reusing that runner would not exercise the drift at all. Mark the job continue-on-error: true until step 8 lands, then make it required.
10. Upstream. Once stable in rabe, offer quadrature_gauss.f90 to the fortnum M4 "integrate" milestone (lazy-fortran/fortnum#34, tracked from the fan-out table in Umbrella: replace avoidable numerical dependencies with AD-ready fortnum libneo#291). Do not block this PR on fortnum: per that table fortnum is private and M4 is gated behind M0 and M1, neither of which has closed.

If step 2 of the discriminating experiment shows eta_b/phi_max also move across toolchains, split that off into its own issue ("refine_maximum argmax is decided at ULP level; global B_max selection is not reproducible") — it belongs with #92 (spline floor of the field-line |B| maxima) and #110 (exact-zero tolerance semantics), not here. #90 should stay scoped to the quadrature.

Acceptance criteria

• nix develop .#default --command make test_all passes, including the new TestIntegrateDeterministic and TestQuadratureGauss.
• TestIntegrate1dSubstituted still passes at its existing reltol = 1.0e-8_dp (test/unit/test_integrate_1d_substituted.f90:8).
• make golden still passes at the unchanged rtol=1e-10, atol=0.0 inside nix develop, against the regenerated record.
• Built with a system gfortran >= 14 / glibc >= 2.39 and compared against the pinned build's fresh output, Lambda_A and remainder are within 1e-8 relative — i.e. down from 4.6e-7 and 6.0e-8 into the same class as the ~3e-9 Lambda_S/lambda_LC_bB floor. No variable exceeds 1e-8.
• The regenerated test/golden/expected/rabe.nc differs from the current one by less than 1e-6 relative on every variable, and the per-variable shift is tabulated in the PR body.
• compare_golden.py accepts an explicit tolerance and still defaults to rtol=1e-10, atol=0.0.
• The golden-unpinned CI job is green at its documented tolerance and no longer continue-on-error.
• The slow-labelled integration tests show no new failures.

Explicitly not an acceptance criterion: an unpinned or macOS make golden passing at rtol=1e-10. The ~3e-9 direct-arithmetic floor makes that unreachable by changing the integrator, and the first comment above notes that nix-on-Darwin still uses Apple libm with x86/ARM ISA-level FP differences. A macOS measurement is worth reporting in the PR body, but wheel.yml builds wheels on macos-14 and does not run the golden test, so wiring one up is separate work.

References

• Piessens, de Doncker-Kapenga, Überhuber, Kahaner, QUADPACK: A Subroutine Package for Automatic Integration, Springer (1983) — the dqag error estimator and the key/limit semantics used at integrate.f90:54-68; key = 6 selects the 61-point Gauss–Kronrod rule.
• Golub & Welsch, "Calculation of Gauss quadrature rules", Math. Comp. 23, 221–230 (1969) — node/weight generation, and the reason step 2 prefers tabulated constants over computing them.
• Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM (2002), ch. 4 (Summation) — compensated summation, the basis for the fixed-order accumulation in step 2.
• Helander, Geiger, Maaßberg, Phys. Plasmas 18, 092505 (2011) — cited in-tree at src/offset/coefficients.f90:14 as the source for the finite-boundary-layer correction Lambda_S.
• Landreman & Catto, Phys. Plasmas 19, 056103 (2012) — cited in-tree at src/shaing_callen/shaing_callen.f90:14 as the source for the lambda_LC formula.
• Landreman & Paul, "Magnetic fields with precise quasisymmetry for plasma confinement", Phys. Rev. Lett. 128, 035001 (2022) — the QH configuration in the golden fixture (wout_LandremanPaul2021_QH_reactorScale_lowres_reference.nc).
• Direct prior art in this project: itpplasma/pyrabe#6 states that the JAX port uses fixed-order Gauss–Legendre n = 32 on the same sqrt-substituted integrands and that "GL n=32 matches QAGS to rel 2e-16 for smooth substituted integrands on a single field line"; itpplasma/pyrabe#23 repeats the measurement for all four integrands. That is the empirical support for the n = 32 choice in step 2, and it also means the deterministic rule has already been exercised against this exact integrand family.
• The SIMPLE compare_golden_results.sh recompute model (third comment above) is the complementary answer at the test-harness level — worth adding alongside, but it does not make a committed record portable.

Scope and risk

One PR for steps 1–9 if the discriminating experiment confirms the quadrature is the dominant amplifier; step 1 (the module split) should be a separate no-op commit inside it. Main risks: a fixed 4-panel/32-node rule under-resolving the t ~ 5e-4 boundary layer on the global-maximum field line (the contribution there is O(1e-16) of the integral, so almost certainly benign, but verify per-field-line rather than only on the surface averages), and the record regeneration in step 8 masking an unrelated regression — hence the requirement to tabulate the per-variable shift. Related open issues: #92 (spline floor of the field-line |B| maxima — this is what sets the residual floor left after the present fix), #110 (toolchain-sensitive exact-zero tolerance), itpplasma/pyrabe#6 and #23 (the GL(32) measurement), itpplasma/libneo#291 (fortnum umbrella).