[ChrisRackauckas-Claude]

Summary

The 1-D SampledIntegralProblem path in evalrule(::AbstractVector, …) reduces weights[i] * data[i] with a sequential @inbounds loop (~1.2 ns/element). That loop is already gated on ArrayInterface.fast_scalar_indexing(data). @simd is applied only for TrapezoidalUniformWeights (branch-free getindex via ifelse); all other weight types — including SimpsonUniformWeights and every NonuniformWeights — keep the sequential loop exactly as on master. Gating on the broader UniformWeights also covered branchy Simpson uniform weights and made SimpsonsRule-on-range ~1.3× slower on Julia 1.10 aarch64. Existing Core tolerances were left unchanged and still pass.

Reassociation (Trapezoidal uniform only; Simpson / nonuniform match master bitwise): vs a sequential reference on sin-like data over N ∈ {10³,10⁴,10⁵,10⁶}, max relative difference is ~9e-15 (1.12) / ~4e-15 (1.10). Vs master on random-sign Float64 data that cancels, differences reach up to 1.9e-13 relative on Julia 1.12 and 2.6e-14 on Julia 1.10. For Float32 at N=10⁶ the Float32 result can differ by up to ~9e-5 relative from master, but the PR is more accurate: trapezoid error vs a Float64 reference drops from 8.3e-5 on master to 6.2e-7.

Verification

Before/after @benchmark minima (JULIA_NUM_THREADS=1 taskset -c 69, amdci2 x86)

case	N	Julia 1.10 before → after	Julia 1.12 before → after
trapz range	10³	1.24 → 0.29 µs	1.23 → 0.31 µs
trapz range	10⁵	120.13 → 23.15 µs	120.12 → 25.15 µs
trapz range	10⁶	1227 → 236 µs	1201 → 250 µs
trapz vector	10⁵ / 10⁶	unchanged (~120 / ~1220 µs)	unchanged
simpson range / vector	all N	byte-identical to master (sequential path)	byte-identical to master

Simpson paths: reinterpret(UInt64, solve!(…).u) == reinterpret(UInt64, sequential_ref) for range and vector grids over N ∈ {8…10⁶} on both Julia 1.10.12 and 1.12.7. No Simpson speedup is claimed.

Correctness: existing rtol=1e-12 sequential-reference testset kept; full sampled_tests.jl on Julia 1.10.12 and 1.12.7 (46 pass).

Test-group tails

GROUP=Core (Julia 1.12.7):

Core/alt_transformation_tests.jl |   19     19  19.6s
Core/gaussian_quadrature_tests.jl |   37     37  5.9s
Core/hadaptiveintegration_tests.jl |    9      9  11.9s
Core/inf_integral_tests.jl |  746    746  3m39.9s
Core/interface_tests.jl |  871    871  8m58.8s
Core/precompile_workload.jl |    1      1  1.1s
Core/quadrule_tests.jl |    4      4  2.3s
Core/sampled_tests.jl |   46     46  9.5s
Core/verbosity_tests.jl |   25     25  1.3s
Testing Integrals tests passed

GROUP=QA (Julia 1.12.7):

QA/qa.jl |   27     27  2m18.6s
Testing Integrals tests passed

Also ran Runic (--inplace) and typos on the changed files.

Not verified

• Full GROUP=AD / ADBatch / Everything.
• Matrix / higher-rank evalrule paths (already had inner @simd; untouched).
• GPU / non-fast_scalar_indexing arrays (still use _evalrule_general).
• Making SimpsonUniformWeights branch-free with nested ifelse (review option b; deferred — would reclaim Simpson speedup on 1.10/1.12).
• Re-timing on aarch64 after this narrower gate (reviewer’s M2 Simpson regression should be gone because Simpson no longer uses @simd).

Reviewer pushback

• Prefer option (b): branch-free SimpsonUniformWeights getindex so Simpson-on-range can also @simd without the 1.10 aarch64 regression.
• Prefer materializing nonuniform weights into a dense Vector once, then @simd.
• Whether reassociation / Float32 accuracy notes belong in the public docs.

Please ignore until reviewed by @ChrisRackauckas.

Risk assessment

Independent review: pending

🤖 Generated with Cursor Agent 2026.10.01-14929f9 (model: unknown (Cursor auto)), transcript /home/crackauc/sandbox/goals/performance/jobs/integrals/int-simd/log.txt on amdci2; orchestrated by Claude Code (claude-opus-5-5[1m]) https://claude.ai/code/session_01LPHREnnonfLg1VcE1EJovv

Independent review: Devin Fusion (fusion-claude-opus-5-5-high-sidekick-swe-2-medium) rated it low, verdict MERGE: #381 (comment)

Independent review (added by the fleet master)

Verdict: MERGE (reviewed head 635644c)

• Risk: low
• Blast radius: CPU scalar-indexable one-dimensional sampled integrals with TrapezoidalUniformWeights. Floating-point reduction order changes; Simpson, nonuniform, matrix and GPU paths are unchanged. No new dependencies, exports, signatures, or documented accuracy guarantees change; no SemVer-breaking contract identified.
• Evidence: Initial diff-only rating was medium, assigned before reading the body; narrowed to low after local checks. Ran JULIA_NUM_THREADS=4 JULIA_NUM_PRECOMPILE_TASKS=4 JULIA_DEPOT_PATH=<local> julia +1 --project=review-env evidence/review381.jl on Julia 1.13.1: all 46 sampled tests and 11 extra Float32/Float64/ComplexF64/BigFloat/Rational/Unitful checks passed. Ran the identical PR sampled test file against base 63b17e8: 46/46 passed, including the added 16 assertions. This is a performance change, so failing-before correctness evidence is inapplicable. Timed the exact base vector method versus head at N=100000: uniform trapezoidal 120.61→25.70 μs; uniform Simpson 244.52→240.75 μs; nonuniform trapezoidal 121.00→120.83 μs; nonuniform Simpson 502.54→502.02 μs. Random-sign N=1000000 data showed the expected reassociation: Float64 absolute difference 1.01e-17; Float32 error against the Float64 sequential reference decreased from 2.05e-8 to 9.05e-9. All 15 reported PR checks pass; latest default-branch CI, downgrade, documentation and formatting runs pass. No skipped assertions, loosened existing tolerances, swallowed errors, removed coverage, or workflow changes. Full Core/QA/AD, LTS, other architectures and GPU tests were not rerun locally.
• Independent review: codex-cli gpt-6.1-sol rated it low, high confidence: narrow SIMD gate produces a measured speedup with passing targeted coverage and green CI.
• Merge: auto-merge candidate — low risk, green checks, independent harness from Cursor author, and no blocking findings.

🤖 Section added by an AI agent — harness: Claude Code 2.1.300 · model: claude-opus-5-5 (review by codex-cli, gpt-6.1-sol)
Conversation: https://claude.ai/code/session_01VZjwf7Ux9YtTzwHGHMqnsq

[ChrisRackauckas]

🤖 Automated comment from an AI agent running as @ChrisRackauckas — not written or reviewed by Chris.

Independent review (Devin CLI 3000.11.3, model fusion-claude-opus-5-5-high-sidekick-swe-2-medium): CHANGES, risk low.

Full review

VERDICT: CHANGES
RISK: low

Blocking findings

1. src/sampled.jl:57: the @simd gate makes SimpsonsRule on uniform grids (AbstractRange) about 1.32× slower on Julia 1.10 (LTS) on aarch64. Confirmed by running.
   The gate is weights isa UniformWeights, so it also covers SimpsonUniformWeights. That type's getindex (src/simpsons.jl:16-25) has early-return branches. This is the same kind of code the PR keeps off @simd for NonuniformWeights, because of a 1.65× regression on 1.10. The PR table says Simpson-on-range gets faster on 1.10, but it was only measured on amdci2 (x86). On Apple M2 it gets slower:

   • Command: julia +lts --project=<env> scratch/bench_simp.jl, run once with master and once with the PR. It takes @belapsed Integrals.evalrule(y, w, 1) three times per size; y = sin.(range(0,1,N)).
   • Julia 1.10.12, master: N=1e3 1.24 µs, 1e4 12.42 µs, 1e5 124.92 µs, 1e6 1250 µs.
   • Julia 1.10.12, PR: N=1e3 1.65 µs, 1e4 16.58 µs, 1e5 165.0 µs, 1e6 1651 µs.
   • All three runs agree to within <1%. On Julia 1.12.7 the same case improves (124.9 → 72.9 µs at 1e5). Trapezoidal-on-range improves on both versions (1.10: 123.8 → 41.6 µs; 1.12: 123.7 → 20.8 µs at 1e5).

   Change required: either

   • (a) restrict @simd to weight types whose getindex is branch-free, i.e. TrapezoidalUniformWeights today, or
   • (b) make SimpsonUniformWeights getindex branch-free with nested ifelse (keep the 17h/48-style expressions so the weights stay bitwise identical).

   For (b), a scratch probe (scratch/simp_variant.jl, a copy of the loop with an ifelse Simpson weight) gives 60.4 µs at N=1e5 on both 1.10 and 1.12, against 124.9 µs sequential. So (b) gets the gain without the LTS regression.

   Either way, re-measure on aarch64 and x86 for 1.10 and 1.12, and correct the PR table. Option (b) touches simpsons.jl, which open PR Fix SimpsonsRule on uniform grids with n < 8 #378 also edits (only find_weights, so it shouldn't conflict).

Non-blocking findings

• Accuracy (confirmed by running, scratch/probe.jl): I compared against a BigFloat sum built from the package's own weights. The PR is as accurate as master or more so in every case I tried:

  • Float64 sin data, N=1e6: relative error 4.3e-15 vs 6.8e-15 on master.
  • Float32, N=1e6: 4.8e-7 vs 1.76e-5.
  • Cancelling randn data: absolute error stays at ≤2.1e-17, same order as master.

  The PR's Float32 numbers (8.3e-5 → 6.2e-7) used different data, so I could not reproduce them exactly. The direction holds.

• Type coverage (confirmed by running): results match master for:

  • Unitful axis and data, Complex, BigFloat, SVector and Vector{Vector} data.
  • A strided view.
  • ForwardDiff through both the data and the range endpoint.
  • n = 1…6, and an empty input (ArgumentError, as before).

  The PR adds no new failure mode.

• The new testset only guards against regressions (test/sampled_tests.jl:35-55). It passes on master too, which is expected for a performance PR. The expected value is computed independently in the test (a sequential sum with the cached weights), not copied from solver output.

  • rtol = 1e-12 is not derived, but it is reasonable: the data are positive and N ≤ 1e4, so worst-case reassociation error is about N·eps ≈ 2e-12, and the observed error is about 1e-16.
  • The test checks against the package's own weights, not the analytic integral. Weight correctness is covered by the existing "Sampled Integration" testset.
  • It reads the internal cache.cacheval. That is acceptable inside the package's own tests.

• Pre-existing, not caused by this PR: on uniform grids with n < 8, SimpsonsRule gives wrong results. For example, n=4 with constant 1 gives 1.0556, identically on master and on the PR. Open PR Fix SimpsonsRule on uniform grids with n < 8 #378 fixes this.

• Sibling PR Fix sampled evalrule silent wrong results on OffsetArrays #379 edits src/sampled.jl:53 (the fast_scalar_indexing condition), 4 lines above this hunk. I expect it to merge cleanly or with a trivial conflict (suspected from reading; I did not test the merge).

• Attribution: the PR body has no 🚧 UNREVIEWED banner as its first line, which the governing CLAUDE.md requires. It has "Please ignore until reviewed" near the bottom instead. The commit trailers are present (Agent-Model: unknown (Cursor auto) is allowed). The repo has no AGENTS.md or CLAUDE.md at its root.

• Scope: one concern, no new dependencies, no new public names, no narrating comments. Runic --check and typos are clean on both changed files (confirmed by running).

What I ran

All runs were on an Apple M2 Max under Julia 1.12.7 and 1.10.12, using scratch envs that Pkg.develop either the PR checkout or a git archive origin/master copy.

• include("repo/test/sampled_tests.jl") on 1.12.7 and 1.10.12: 24 + 16 + 1 + 5 = 46 pass on each. This matches the PR's claim of 46.
• scratch/bench.jl on master and the PR, both Julia versions: trap/simp × range/vector × N ∈ {1e3, 1e5}, plus Float32 range. Vector (nonuniform) timings are unchanged, as the PR claims.
• scratch/bench_simp.jl: Simpson on a range, N = 1e3…1e6, three repeats, master vs PR, both versions (the blocking finding).
• scratch/simp_variant.jl: sequential vs @simd loop with the package's Simpson weights and with an ifelse version, both versions.
• scratch/probe.jl on master and the PR (Julia 1.12), then diff: accuracy against a BigFloat reference and the type-coverage cases above.
• julia +1.12 --project=@runic -m Runic --check --diff and typos on src/sampled.jl and test/sampled_tests.jl.
• gh pr view 381, gh pr list, and gh pr diff 378 / 379 for sibling work. The PR links no issue and no tracking issue.

What I did not verify

• The PR's GROUP=Core (about 13 min) and GROUP=QA results. I did not rerun them because of host limits. Only sampled.jl and its test file change.
• GROUP=AD / ADBatch. Sampled-integration AD paths do not call evalrule directly (grep of ext/ found no use besides a TODO). This is suspected from reading only.
• Timings on x86 (amdci2), and the PR's claim that Simpson-on-range timings on master are bimodal.
• Whether Fix sampled evalrule silent wrong results on OffsetArrays #379 actually merges cleanly with this PR.
• GPU arrays (they take _evalrule_general, which this PR does not change).

---

🤖 Posted by an AI agent — harness: Devin CLI 3000.11.3 (review), Claude Code 2.1.285 (posting) · model: fusion-claude-opus-5-5-high-sidekick-swe-2-medium
Conversation: local Claude session 3c311569-dda6-4687-bc9f-65febb212c33 (fleet master)

[ChrisRackauckas]

🤖 Automated comment from an AI agent running as @ChrisRackauckas — not written or reviewed by Chris.

Independent review (Devin CLI 3000.11.3, model fusion-claude-opus-5-5-high-sidekick-swe-2-medium): MERGE, risk low.

Full review

VERDICT: MERGE
RISK: low

PR: #381 (head 635644c, base origin/master e02a9f9). No linked issue in the PR body, and no AGENTS.md/CLAUDE.md at the repo root. I checked against ~/.claude/CLAUDE.md instead.

The change adds @simd to the 1-D evalrule reduction only when weights isa TrapezoidalUniformWeights (src/sampled.jl:57-65). It also adds one testset (test/sampled_tests.jl:35-55). Every other weight type runs the same loop as before.

Blocking findings

None.

Non-blocking findings

1. The new test cannot fail on master, and 3 of its 4 cases check nothing. (confirmed by running)

   • I ran test/sampled_tests.jl against master's src. The new testset passes there too (16/16). That is expected for a pure-performance PR. But it means the test only guards the SIMD loop against gross errors such as a wrong start index or a dropped term. It does not show that anything changed.
   • For (xv, Trapezoidal), (x, Simpson) and (xv, Simpson), the reference is the same sequential loop the code runs, so those assertions are true by construction. Only (x, TrapezoidalRule()) exercises the new path. It would be cleaner to keep just the range/trapezoid case.
   • rtol = 1e-12 is not derived from anything. The worst-case reassociation bound at N=10⁴ is N·eps ≈ 2.2e-12. It passes here because the data is positive and well-conditioned (observed differences are about 1e-15). A one-line justification would help, or use a tolerance scaled by N·eps.

2. Rounding now depends on the element type. (confirmed by running) On the PR, Julia 1.10 and 1.12:

   • Float64 y = x.^2, N=101: 0.33335000000000004.
   • The same integral with Complex, Unitful, Real[], ForwardDiff Dual or vector-of-vectors data: 0.33335000000000015. Those types don't vectorize and stay sequential.
   • On master every type gave ...015.
   • So the ForwardDiff primal/partials and the plain Float64 result now differ by about 1 ulp. This is harmless, but it is a behavior change worth one sentence in the PR body.

3. Float64 results on uniform trapezoid grids change by up to ~1e-13 relative; Float32 by up to ~2.5e-5. (confirmed by running) Against a BigFloat reference, the new results are more accurate everywhere I measured:

   case (Julia 1.12)	master rel. error	PR rel. error
   Float32, N=10⁶, sin	1.1e-4	1.0e-6
   Float32, N=10⁶, randn	9.4e-5	2.2e-5
   Float64, N=10⁶, randn	6.7e-14	1.1e-14

   One case went the other way: Float32, N=10⁴, sin went from 1.48e-7 to 1.72e-7 on 1.12 (on 1.10 it improved to 4.1e-8). That is the same order of magnitude. Any downstream user who compares sampled-trapezoid output with == will see diffs. That is acceptable for a patch release, and the PR body already says so.

4. Style. (from reading) weights isa TrapezoidalUniformWeights is a branch on a concrete type inside generic code. It is resolved at compile time, so it costs nothing at runtime. A small dispatched helper (e.g. _reduce1d(data, w::TrapezoidalUniformWeights)) would be more idiomatic and would make extending it to a branch-free SimpsonUniformWeights later trivial. Optional.

5. Interaction with sibling PRs. (from reading the diffs, not merge-tested)

   • Fix sampled evalrule silent wrong results on OffsetArrays #379 changes the if ArrayInterface.fast_scalar_indexing(data) line of the same function (line 53). It does not touch the lines changed here.
   • Fix SimpsonsRule on uniform grids with n < 8 #378 and Fix sampled evalrule silent wrong results on OffsetArrays #379 both add testsets next to where this PR inserts its testset in test/sampled_tests.jl. Whichever lands second may need a trivial rebase.
   • Fix SimpsonsRule on uniform grids with n < 8 #378 makes SimpsonUniformWeights.getindex more branchy. That fits this PR's decision to keep Simpson out of the @simd path.

6. No narrating comments, no new dependencies, no API change. (confirmed by reading the diff) Runic and typos are clean.

What I ran

All on an Apple M2 Max (aarch64), JULIA_NUM_THREADS=1. "Master" means a separate clone at origin/master e02a9f9, loaded through its own scratch environment.

• Sampled tests: include("test/sampled_tests.jl") on the PR, Julia 1.12.7 and 1.10.12: 24 + 16 + 1 + 5 = 46 pass, matching the PR body's 46. The same file against master's src on 1.12.7: all pass, including the new testset.
• Probe script (scratch/probe.jl), run on PR and master under both Julia versions:

  • Bitwise outputs: 168 results over N ∈ {8, 9, 100, 1001, 10⁴, 10⁵, 10⁶} × {Float64, Float32} × 3 datasets (sin; randn; 1e8·randn+1, which cancels) × {range, vector} × {Trapezoidal, Simpson}. Every non-range/trap result is bit-identical to master on both 1.12 and 1.10, which confirms the PR's "Simpson / nonuniform match master bitwise" claim. All 34 range/trap results differ. Largest Float32 relative difference from master: 2.5e-5.

  • Accuracy vs BigFloat: the trapezoid errors summarized in finding 3.

  • Generic element types on the SIMD path: Unitful (m·s), Complex, BigFloat, Int data, Real[], vector-of-vectors and ForwardDiff Dual (Float64 and Float32) all work on the PR and give the same values as master. The only difference is the 1-ulp Float64 primal from finding 2.

  • Zygote: gradient through solve(...)[1] fails identically on master and the PR ("Need an adjoint for constructor IntegralSolution"). That failure was already there and has nothing to do with this PR.

  • Timings (@belapsed solve!(cache)):

    case	Julia 1.12 master → PR	Julia 1.10 master → PR
    trap range, N=10³	1.16 → 0.23 µs	1.16 → 0.39 µs
    trap range, N=10⁵	123.8 → 20.9 µs	124.9 → 41.6 µs
    trap range, N=10⁶	1250 → 215 µs	1239 → 417 µs
    trap vector, simp range, simp vector (all N)	unchanged	unchanged

    That is roughly 5.8× on 1.12 and 3× on 1.10. No Simpson regression on aarch64, so the regression the PR mentions for its earlier UniformWeights gate is gone.

• Formatting and spelling:
  • julia +1.12 --project=@runic -m Runic --check --diff src/sampled.jl test/sampled_tests.jl gives exit 0.
  • git diff origin/master...HEAD | typos - gives exit 0.

What I did not verify

• The x86 (amdci2) timings in the PR body. I only timed on aarch64, where the direction and size are consistent.
• The GROUP=Core and GROUP=QA tails in the PR body. I ran only sampled_tests.jl, the one file that touches this code. QA (Aqua/JET/ExplicitImports) cannot plausibly be affected by a loop annotation.
• AD test groups, GPU / non-fast_scalar_indexing arrays (that path is unchanged), and Julia 1.11.
• Whether the PR merges cleanly on top of Fix SimpsonsRule on uniform grids with n < 8 #378, Fix sampled evalrule silent wrong results on OffsetArrays #379 or Stabilize sampled eachslice integration with a function barrier #380. I read their diffs but did not merge them.

---

🤖 Posted by an AI agent — harness: Devin CLI 3000.11.3 (review), Claude Code 2.1.285 (posting) · model: fusion-claude-opus-5-5-high-sidekick-swe-2-medium
Conversation: local Claude session 3c311569-dda6-4687-bc9f-65febb212c33 (fleet master)