[ChrisRackauckas-Claude]

Summary

SimpsonsRule on an AbstractRange used the uniform 17/59/43/49 end corrections from both ends even when length(x) < 8, so those corrections overlapped and integrals of constants / x^2 on [0,1] were wrong (e.g. ∫1 ≈ 0.96875 for n=3). This keeps find_weights returning SimpsonUniformWeights for every AbstractRange length (cache-stable / @inferred), and implements correct composite Simpson weights for 3 ≤ n < 8 inside getindex: classic 1/3 for odd n, and 3/8 on the first three panels plus 1/3 on the remainder for even n (n == 4 is pure 3/8). For n ≥ 8 the previous 17/59/43/49 weights are bitwise unchanged. n < 3 on a range throws ArgumentError.

A 2-point uniform Range previously returned a wrong value (∫1 over [0,1] ≈ 0.708, absolute error ≈ −0.29) and now throws ArgumentError("The length of the grid must exceed 2 for simpsons rule.") from find_weights (Vector path still asserts in getindex).

Verification

• Original bug (Julia 1.12.4): unfixed src + small-n tests → failures for n = 3..7; with fix → pass (before.txt / after.txt in the job dir).
• F2 cache regression: with the intermediate type-split approach, init at n = 5 then cache.x length 11 threw TypeError in setfield!; with this commit, 5→11 and 11→5 reuse both return 1.0 and match a fresh solve.
• F2 checks: @inferred find_weights on a range; bitwise identity of SimpsonUniformWeights vs the old formula for n = 8..20; correctness for n = 3..10 on Float64/Float32/BigFloat for 1, x^2, and x^3 (142/142 in the focused suite). Julia 1.10.12: same 142/142.

Test-group tails (after F2)

GROUP=Core (Julia 1.12.4):

Core/alt_transformation_tests.jl |   19     19  19.7s
Core/gaussian_quadrature_tests.jl |   37     37  5.5s
Core/hadaptiveintegration_tests.jl |    9      9  11.9s
Core/inf_integral_tests.jl |  746    746  3m35.9s
Core/interface_tests.jl |  871    871  8m53.1s
Core/precompile_workload.jl |    1      1  1.2s
Core/quadrule_tests.jl |    4      4  2.4s
Core/sampled_tests.jl |  158    158  11.3s
Core/verbosity_tests.jl |   25     25  1.2s
Testing Integrals tests passed

GROUP=QA (Julia 1.12.4):

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

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

Not verified

• Full GROUP=AD / ADBatch / Everything.
• GPU / non-fast-scalar-indexing arrays with small ranges.
• Perf of small-n getindex vs the old incorrect uniform path.

Reviewer pushback

• Prefer 3/8 at the trailing end (or both ends) instead of leading 3/8 for even n < 8 (same polynomial degree; different weight layout).
• Whether n = 2 should keep a numeric (wrong) result for back-compat instead of ArgumentError.
• Docs / docstring for SimpsonsRule should mention the n ≥ 3 requirement and the small-n vs n ≥ 8 weight formulas.

Please ignore until reviewed by @ChrisRackauckas.

Risk assessment

Independent review: F1 MERGE (Claude Code, Opus 5.5); F2 addressed Devin CHANGE finding (cache type split) — re-review pending.

🤖 Generated with Cursor Agent 2026.10.01-14929f9 (model: unknown (Cursor auto)), transcript /home/crackauc/sandbox/goals/performance/jobs/integrals/int-simpson/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: #378 (comment)

Risk assessment

• Risk: low
• Blast radius: Uniform-range Simpson integration, chiefly grids of 3–7 samples; invalid ranges with fewer than three samples now throw ArgumentError. No public names, signatures, or runtime dependencies change.
• Evidence: All 15 reported PR checks pass. Locally, the head sampled test file passes 158/158 assertions. An independent weight-level polynomial probe passes 216/216 cases (Float32, Float64, BigFloat; n=3:20; degrees 0:3); the same probe on the PR base fails 60 assertions. The cache-crossing tests pass in both directions.
• Independent review: codex-cli gpt-6.1-sol rated it low, high confidence: correct small-grid coefficients, stable cache weight type, analytic regression coverage, and passing sampled tests and CI.
• Merge: auto-merge candidate — narrow internal correction with targeted regression tests and all reported checks passing; no blocking findings.

Independent risk assessment by codex-cli (gpt-6.1-sol), reviewed head 42e98d447f; dispatched by the fleet master, Claude Code (claude-opus-5-5), 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 medium.

Full review

VERDICT: CHANGES
RISK: medium

Review of #378 (head 390cb70). The PR links no issue. Repo root has no AGENTS.md or CLAUDE.md; I applied the global CLAUDE.md rules.

The numerical fix is correct. On master, a uniform AbstractRange with n < 8 gets overlapping 17/59/43/49 end corrections, and with this PR such ranges now give results bitwise-identical to the Vector path. The problem is how it gets there: find_weights now returns a different concrete weight type depending on the range's length at runtime. That breaks the documented cache-reuse interface.

Blocking findings

1. src/simpsons.jl:84-89: find_weights on an AbstractRange returns a different type depending on length. A reused SampledIntegralCache now throws TypeError when the grid crosses n = 8 (confirmed by running).

   • Why it breaks: SampledIntegralCache{…, Tc, …} stores cacheval::Tc as a concrete type fixed at init (src/common.jl:264-274). __solvebp_call then does cache.cacheval = find_weights(grid, alg) after cache.x is replaced (src/sampled.jl:149). The init docstring (src/common.jl:~305) documents this case: "Updating the sample locations marks the cached weights stale so they are recomputed before the next solve."

   • Probe: scratch/probe.jl and scratch/cache_master.jl. Each one calls init on range(0,1;length=n1), then sets cache.x = range(0,1;length=n2) and cache.y = ones(n2), then calls solve!:

     n1 → n2	master	PR
     5 → 11	1.0 (correct)	TypeError: in setfield!, expected SimpsonNonuniformWeights{StepRangeLen…}, got a value of type SimpsonUniformWeights…
     11 → 5	1.0156250000000002 (wrong, the bug being fixed)	TypeError: … expected SimpsonUniformWeights{Float64}, got … SimpsonNonuniformWeights…
     5 → 6	0.99166… (wrong)	1.0

     The 5 → 11 row is a strict regression: master returned the correct value and the PR throws.

   • Also a Union return type: Base.return_types(Integrals.find_weights, (StepRangeLen…, SimpsonsRule)) is now Union{SimpsonNonuniformWeights{StepRangeLen…}, SimpsonUniformWeights{Float64}}, so init(prob, SimpsonsRule()) for a range is inferred as a Union of two cache types. JET @report_opt still reports 0 because of union splitting, so QA does not catch it. Anyone storing the cache in a struct field typed from init will see this.

   • Required change: find_weights(::AbstractRange, ::SimpsonsRule) must return a single concrete type for every length. One way is to keep returning SimpsonUniformWeights and make its getindex correct for 3 ≤ n < 8 using closed-form uniform weights:

     • odd n: h/3·[1,4,2,…,4,1];
     • even n: either the weights SimpsonNonuniformWeights produces on a uniform grid (for n = 4 that is h·[1/3, 5/4, 1, 5/12], which I checked sums to 3h), or a 3/8 tail that is exact for cubics;
     • add an n < 3 check.

   • Required test: add a cache-reuse case to the "Caching interface" testset that crosses the threshold in both directions (for example 5 → 11 and 11 → 5) and compares against a fresh solve.

Non-blocking findings

1. Docstring out of date (src/algorithms_sampled.jl:50-51, suspected from reading). It still says AbstractRange grids "use the composite Simpson 1/3 and 3/8 rules", which no longer describes n < 8 under this PR. It also does not state the n ≥ 3 requirement. The PR body lists this as reviewer pushback but leaves it unfixed.
2. New behavior on a public path without a test. 2-point Range input now throws AssertionError (confirmed: master returns 0.7083… for ∫1 on range(0,1;length=2), the PR throws). Throwing is better than a silent wrong answer. Two problems:
   • nothing tests it;
   • input validation via @assert (pre-existing, inside getindex) can be compiled out, so an ArgumentError in find_weights would be more appropriate.
3. Accuracy drops for even n < 8 (by design, confirmed). For n = 4 and n = 6, ∫x³ has errors 3.1e-3 and 4.0e-4. That matches the Vector path exactly, so it is not a regression versus the stated goal. But a closed-form 3/8-based rule would be exact for cubics there, as the n ≥ 8 uniform rule already is. If you implement the fix in blocking finding 1, use that rule.
4. Test quality is fine. Expected values are analytic (1 and 1/3), not code output. Default ≈ tolerance is appropriate because both integrands are integrated exactly. The test fails on master (confirmed: 12 pass / 20 fail in the new testset) and passes with the fix. It only probes constants and x², so it cannot distinguish rules of different polynomial degree.
5. Interaction with sibling PRs:
   • SIMD the 1-D sampled-data weight×data reduction #381 adds @simd only when weights isa UniformWeights. With this PR, small ranges leave that branch (perf only, and moot if blocking finding 1 is fixed as suggested).
   • Fix sampled evalrule silent wrong results on OffsetArrays #379, Stabilize sampled eachslice integration with a function barrier #380 and SIMD the 1-D sampled-data weight×data reduction #381 all add testsets at the same spot in test/sampled_tests.jl (after line 32), so expect textual merge conflicts.
6. PR body. It is accurate on the numbers I reproduced. It contains no command lines, only result tails. The before.txt/after.txt files it cites live on amdci2 and are not attached.

What I ran

All runs used Julia 1.11.9 with scratch envs in ../scratch (env = PR head, envm = git archive origin/master).

• scratch/probe.jl (PR head), per-n errors for n = 2..12 on Range vs Vector:
  • n < 8: Range and Vector are bitwise equal (checked with exp data on a shifted range).
  • Also correct: reversed ranges, 1:5, 1:4, Float32, BigFloat, Unitful ranges, and Matrix data with dim = 1 and 2.
  • The cache-crossing failure and the Union return type in blocking finding 1 come from this run.
• scratch/cache_master.jl (master): the cache results above, and master ∫1 for n = 2..7: 0.7083, 0.96875, 1.0556, 1.0156, 0.99167, 0.99653.
• The PR's test/sampled_tests.jl:
  • against master src: 36 pass / 20 fail; all failures are in the new testset (12/20);
  • against PR src: 56/56, plus "dimensional axes" 1/1 and "Caching interface" 5/5.
• scratch/jet.jl: JET @report_opt of solve and of solve!(init(…)) with SimpsonsRule on ranges of length 101 and 5 gives 0 reports. return_types(init, …) shows the Union of two cache types.
• Formatting and spelling:
  • julia +1.12 --project=@runic -m Runic --check --diff src/simpsons.jl test/sampled_tests.jl: clean.
  • typos on the diff: only a false positive on a git hash.
• Read gh pr diff for Fix sampled evalrule silent wrong results on OffsetArrays #379, Stabilize sampled eachslice integration with a function barrier #380 and SIMD the 1-D sampled-data weight×data reduction #381.
• No Julia processes of mine are left running.

What I did not verify

• I did not rerun full GROUP=Core (the PR claims a pass), GROUP=QA beyond the SimpsonsRule JET check, or AD/ADBatch.
• I did not test on Julia 1.10 or 1.12.
• I did not check GPU or non-fast-scalar-indexing arrays with small ranges (they go through _evalrule_general; suspected unaffected).
• I did not benchmark the small-n nonuniform fallback.

---

🤖 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: #378 (head 42e98d4, merge-base 63b17e8). No linked issue, and there is no AGENTS.md or CLAUDE.md at the repo root.

Blocking findings

None.

Non-blocking findings

1. src/algorithms_sampled.jl:47-70 (SimpsonsRule docstring) and the SimpsonUniformWeights docstring at src/simpsons.jl:1-11 are unchanged. They don't say that a range grid needs at least 3 points, or that for n < 8 the weights are classic 1/3 (odd n) or a leading 3/8 block plus 1/3 (even n). The PR body raises this itself. One sentence in the docstring would cover it. (From reading.)
2. Error types are inconsistent for short grids. An AbstractRange with n < 3 now throws ArgumentError (src/simpsons.jl:101-102), but a Vector with n = 2 still throws AssertionError, a behaviour that predates this PR. Confirmed by running: solve(SampledIntegralProblem(ones(2), [0.0, 1.0]), SimpsonsRule()) raises AssertionError. Turning a 2-point range from a wrong number (34h/48, i.e. ≈0.708 for ∫1 on [0,1]) into an error is a fix, not a breaking change. Ideally both paths would throw the same exception type, but that can be a follow-up.
3. For even n < 8, the range path puts the 3/8 block at the start. The Vector path puts its odd-panel correction at the end. Both are exact for cubics, so the test u_range ≈ u_vec (which only uses 1 and x^2) passes. For general smooth data the two paths now give slightly different O(h^4) results. That is not a bug, just a design choice worth knowing about.
4. In test/sampled_tests.jl:34, the comment "end corrections overlap for n < 8 without small-n weights" describes the bug being fixed rather than what the test checks. This is minor.
5. The new tests use the default ≈ tolerance (rtol = √eps) even though these rules are exact for degree ≤ 3. That isn't fitted to observed error, and the assertions do fail on the unfixed code (see below). They could be tightened, but they aren't wrong.
6. In test/sampled_tests.jl, PR 381 (SIMD the 1-D sampled-data weight×data reduction #381) inserts code at the same place this PR does. git merge-tree reports no conflict with 379, 381, or current origin/master.

What I ran

All of these are confirmed by running, on Julia 1.12.7 unless noted. The PR environment is repo dev'd; the base environment is a copy with src/ at the merge-base.

• PR's test/sampled_tests.jl with the PR source: Sampled Integration 146/146, dimensional axes 1/1, Caching interface 11/11.
• Same test file on the base source (fails without the fix, as expected):
  • n=3: ∫1 gives 0.96875, ∫x² gives 0.33073, ∫x³ gives 0.25391.
  • n=4: ∫1 gives 1.0556.
  • n=5: ∫1 gives 1.0156.
  • n=6: ∫1 gives 0.99167.
  • So the new tests do fail before the fix and pass after it.
• Same test file on Julia 1.10.12 (LTS) with the PR source: 146/146, 1/1, 11/11.
• scratch/probe.jl, an independent check I wrote:
  • Exact weights. I built SimpsonUniformWeights(n, big(1)//1) for n=3..30 and compared it against reference weights I assembled panel by panel (odd n: 1/3 panels; even n<8: 3/8 on panels 1–3 then 1/3 panels; n≥8: 17/59/43/49 over 48). It matches exactly in rational arithmetic. Σw = n−1, and Σw·xᵈ equals the exact integral for d=0..3. collect(w) (the iterate path) equals getindex. 196/196.
  • n ≥ 8 unchanged. Weights are bitwise identical (===) to the old formula for n=8..200 in Float64 and Float32. 386/386.
  • Odd n. Range and vector results agree to rtol 1e-14 for ∫sin on [0,π] (same classic rule).
  • Smooth-integrand accuracy. Errors for ∫sin on [0,π] are 0.094, 0.041, 0.0046, 0.0034, 0.00086 for n=3..7. They shrink monotonically and are the expected sizes (n=3 gives exactly 2π/3, as classic Simpson should).
  • Other grid and data shapes (42/42): decreasing range (∫₂⁻¹x³ = −3.75), Unitful range, UnitRange{Int}, matrix data with dim=2 and dim=1, 3-D data with dim=3, all for n=3..9.
  • Short ranges. ArgumentError for range length 0, 1, 2.
  • Types. @inferred solve at n=5, 6, 101 passes, and a Float32 grid keeps a Float32 result.
  • My first probe run had a wrong reference value for the decreasing range (I wrote −4.25 where −3.75 is correct). That was my error and I fixed it; the code's value was correct.
• scratch/jet.jl: @report_opt target_modules=(Integrals,) on solve(…, SimpsonsRule()) gives 0 reports at n = 3, 4, 5, 6, 7, 8, 101. This repeats the QA check at small n too.
• scratch/bench.jl (-O3, best of 7, two alternating runs per environment), base vs PR: n=5 22.6–23.2 ns vs 23.4–23.7 ns; n=101 135 vs 135–136 ns; n=1e4 12.53 µs for both; n=1e6 1.25 ms for both. No measurable regression from the n >= 8 branch.
• Formatting and spelling:
  • julia +1.12 --project=@runic -m Runic --check --diff on both changed files: clean.
  • typos on the diff: one hit, inside a git blob hash (4ba5641), which is a false positive.
• Sibling PRs: git merge-tree of this PR against 379, 381, and origin/master reports no conflicts. Neither sibling touches src/simpsons.jl. PR 381's SIMD path applies only to TrapezoidalUniformWeights, so it doesn't interact with this change.

What I did not verify

• I did not rerun the PR body's GROUP=Core (871-test interface suite and others) or GROUP=QA tails. I re-ran only the affected sampled_tests.jl file plus the SimpsonsRule JET check from QA. The rest of Core doesn't call SimpsonUniformWeights: per grep, the only other SimpsonsRule use is a BigFloat case in interface_tests.jl with a long grid.
• AD / ADBatch groups. The ext/ code has no reference to Simpson weights (grep), so I don't expect any effect, but I didn't run them.
• GPU and other non-fast-scalar-indexing arrays. These go through _evalrule_general → iterate. I checked that iterate matches getindex on CPU only.
• Docs build. The PR changes no docstrings.

---

🤖 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)