integrate_1d_gauss_kronrod gradient computation can be very slow

Summary

The gradient computation for integrate_1d_gauss_kronrod call the Boost adaptive Gauss Kronrod again, but Boost allows only relative tolerance and that can lead the adaptation to spend much more time than what the quadrature for the actual integrand takes.

Problem

integrate_1d_gauss_kronrod_tol computes parameter gradients through
internal::integrate_1d_adjoint, which integrates ∂f/∂θ_i separately
for every var argument with the same integrator and tolerance as the
value (rev/functor/integrate_1d_adjoint.hpp, integrate_grad).
Boost's gauss_kronrod::integrate decides whether to bisect a leaf
against tol × |estimate of that leaf| (gauss_kronrod.hpp,
recursive_adaptive_integrate).

This is a purely relative rule. A gradient component whose integrand
is tiny relative to f and dominated by round-off can never satisfy
it: the K21 − G10 disagreement on noise is O(1) relative, so that
component bisects to max_depth (2^15 leaves × 21 nodes by default) on
every call, although its contribution to ∂ log I/∂θ_i = (∫∂f/∂θ_i)/I
is negligible. Stan's own post-check, error ≤ max(rel_tol × L1, abs_tol), is applied only after Boost returns, and Boost's public
integrate hard-codes the recursion's absolute budget to 0, so Stan's
absolute_tolerance cannot reach the bisection.

Where this arises in practice: a marginal likelihood
∫ N(z) ∏_k BetaBinomial(y_k | T_k, φ·logit⁻¹(η_k), φ·logit⁻¹(−η_k)) dz
with saturated observations (y_k = 0). In a tail panel where
f ~ 1e-12, the φ-score of a saturated observation is a difference of
large digamma terms that cancels catastrophically, so ∂f/∂φ ~ 1e-23
and noisy. Measured per subject with an evaluation counter (unpatched):
value 273 evaluations; gradients 668 850, of which 666 393 are that one
component in that one panel; the other seven components converge in 21.

integrate_1d (double-exponential) and integrate_1d_double_exponential are not affected: tanh_sinh /
sinh_sinh terminate on err ≤ tol × L1 at every level.

Measured on a beta-binomial varying-intercept model (105 per-subject
marginals, ~9 var arguments each; one gradient of the whole model at
the same unconstrained point, n = 10; CmdStan 2.39, -O3):

In generated quantities (no gradients) the same GK integral costs 2.1×
DE, so the 134× is entirely the gradient path.

Fix

The quantity consumed downstream is ∂ log I/∂θ_i = (∫ ∂f/∂θ_i) / I,
so each gradient integral only needs absolute accuracy rel_tol × I.
Give Boost's relative rule that absolute scale by integrating

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 is constant (not 0 or 1). Every leaf estimate of g_i is then ≈ the leaf's value
mass, so a noisy component that is small relative to f is accepted
as soon as the value's leaf is, and the bisection stops where the value
integral's did. The subtraction's cancellation error is rel_tol × I,
the accuracy the value has anyway. Guarded on I finite and non-zero
(otherwise unchanged behaviour), so a signed f with zero integral is
no worse than today.

Start in rev/functor/integrate_1d_adjoint.hpp at integrate_grad, then read Boost's gauss_kronrod.hpp and recursive_adaptive_integrate behavior. Compare the gradient path with integrate_1d using the saturated beta-binomial case described in the issue; done means avoiding excessive evaluations while preserving the reported gradient accuracy and existing behavior for non-finite or zero value integrals.