fix integrate_1d_gauss_kronrod gradient speed

[avehtari]

Making this PR has been assisted by Claude. I have checked and edited all new code, comments, and PR text

Summary

The gradient computation for integrate_1d_gauss_kronrod calls the Boost adaptive Gauss Kronrod again, but Boost computes each parameter's gradient integral with the same Boost adaptive Gauss–Kronrod call. Boost's bisection criterion is relative to the leaf estimate, so a gradient component that is tiny relative to f and dominated by round-off (e.g. the score of a saturated observation deep in a likelihood tail) never converges and bisects to max_depth (2¹⁵ × 21 evaluations) on every call, although its contribution to ∂ log I is negligible. This can be fixed by shifting the integral integrand as
g_i(x) = ∂f/∂θ_i (x) + c f(x)
and returning ∫ g_i − c I, where I is the value integral already computed and c != 0 is a constant
such that ∂f/∂θ_i ≈ −c·f does not occur for common integrands; c = 1 is a bad choice because log-density gradients of logit-type likelihoods saturate at ±1 in the tails, which turns the shifted integrand back into noise. The inverse golden ratio is far from every integer and is not a structural constant of any link.

As this is needed only for integrate_1d_gauss_kronrod, integrate_1d_adjoint gets an option whether the shift is used. The plain integrate_1d and integrate_1d_double_exponential are unchanged: their L1-relative termination handles ordinary gradient integrands, and enabling the shift there costs 7–32 % more gradient evaluations on 10 of 35 suite integrands with no gain.

Tests

Measured with the fix (same model, same point):

	gradient	vs integrate_1d	max |Δ grad| vs integrate_1d
GK, unpatched	7582 ms	134×	2.9e-11
GK, patched	156 ms	2.75×	2.9e-11

49× faster gradients, results unchanged to 1e-11 (patched vs unpatched
GK: 1.0e-11). Per subject the counter shows 668 850 → 2 478 gradient
evaluations. The residual 2.75× is GK's per-integral cost on this
integrand (value-only GQ is 2.1× DE), not a gradient effect.

Effect on every integrand in the existing GK reverse-mode test suite

Integrand evaluations (value alone; all gradient components together)
and wall for value + gradients (mean of 200, -O3), unpatched →
patched. Same arguments and tolerances as the tests. #var = number
of var inputs (parameters + var endpoints).

integrand (test)	#var	value	gradients unpatched → patched	wall µs
f1 exp(x)+θ, [0.2, 0.7]	1	21	21 → 21	1.3 → 1.2
f2 exp(θ cos 2πx)+θ, var endpoints	3	21	23 → 23	2.9 → 2.9
f3 exp(x)+θ₀^2.5+2θ₁³+2θ₂, var endpoints	5	21	65 → 65	11 → 7.6
f3, params only	3	21	63 → 63	11 → 5.7
f10 1/(1+x⁴/1.7), (−∞, ∞), no params	0	147	0 → 0	6.3 → 2.7
f11 beta(5, 3), [0, 1]	2	21	294 → 336	29 → 33
f11 beta(3, 5), [0, 1]	2	21	294 → 336	20 → 23
f12 N(5.7, 1), (−∞, ∞)	2	357	672 → 672	44 → 42
f13 x+θ₀+θ₁, endpoints = params	4	21	44 → 44	1.2 → 1.2
Cauchy pdf, (−∞, ∞)	2	231	546 → 462	38 → 33
Exponential pdf, [0, ∞)	1	105	231 → 147	11 → 7.6
Normal pdf, (−∞, ∞)	2	441	1008 → 924	57 → 53
Student-t pdf, (−∞, ∞)	3	651	3129 → 2205	323 → 241
Uniform pdf, endpoints = params	4	21	44 → 44	2.5 → 2.3
noisy gradient (new test), [−1, 1]	1	21	throws → 21	— → 2.6
odd, zero integral (new test)	1	21	(n/a) → 21	— → 0.8
negative integral (new test)	1	21	(n/a) → 21	— → 0.8

Unpatched, the noisy-gradient case does not merely hit max_depth: it
then fails Stan's post-check (error estimate 2.4e-13 exceeds max(rel_tol × L1, abs_tol) with L1 ≈ 0) and throws, i.e. in a
model block every proposal is rejected. Elsewhere the change is
neutral or a mild gain (heavy-tailed pdfs 15–35% fewer gradient
evaluations); the one regression is f11 (beta kernel: the gradient
integrands carry a log x factor that is easier than f on its own),
+14 % gradient evaluations at unchanged wall. All 17 pre-existing
tests pass at their tolerances.

A gradient component that is analytically zero but autodiffs to
round-off noise reproduces it without any likelihood machinery
(cos² + sin² = 1). Added to test/unit/math/rev/functor/integrate_1d_gauss_kronrod_test.cpp
as GradientShift_noisy_gradient (evaluation-count bound + accuracy;
fails on the unpatched code, passes on the patched) together with
GradientShift_guards (exactly-zero integral → shift disabled; negative
integral → shift active), and the existing 17 tests pass unchanged.

Side Effects

1. Cost when a gradient integrand is easier than f. `f = spike(x)
   • θ`: the gradient integrand is the constant 1 and took 21
     evaluations; shifted it re-resolves the spike, 189. Bound: gradient
     cost ≤ (p + 1) × value cost, against 2^15 × 21 per component before.
2. Absolute, not relative, gradient accuracy. The shifted gradient
   inherits the value integral's absolute error (≲ tol × L1(f)), so a
   component with |∂I/∂θ| ≲ tol × L1(f) loses relative precision.
   Measured (f = spike + ε θ x, spike width 0.01, tol = 1e-6):
   unshifted exact; shifted absolute error 3e-19 … 3e-18, i.e. relative
   6e-15 at ε = 1e-3 down to 5e-10 at ε = 1e-9 — round-off level in
   practice, because leaves that pass the value's test are usually far
   below tol. For HMC this is the right accuracy: value and gradient
   enter the Hamiltonian at the same absolute scale.

Release notes

integrate_1d_gauss_kronrod gradient computation was modified to have much improved worst performance timing

Checklist

• Copyright holder: (fill in copyright holder information)

Aki Vehtari

• the basic tests are passing

  • unit tests pass (to run, use: ./runTests.py test/unit)
  • header checks pass, (make test-headers)
  • dependencies checks pass, (make test-math-dependencies)
  • docs build, (make doxygen)
  • code passes the built in C++ standards checks (make cpplint)

• the code is written in idiomatic C++ and changes are documented in the doxygen

• the new changes are tested