add integrate_1d_gauss_kronrod and integrate_1d_double_exponential

[avehtari]

Closes #3000.

EDIT 2026-05-22: PR text updated to note variadic API

Summary

Adds two new explicit-name 1-D quadrature functions, both wrapping
existing Boost quadrature routines. They cover disjoint integrand
weaknesses and live side by side with the existing integrate_1d:

Function	Backend	Strength	Weakness
integrate_1d (existing)	tanh_sinh / exp_sinh / sinh_sinh (DE)	endpoint singularities	off-centre peaks, no abs_tol
integrate_1d_double_exponential (new)	same as integrate_1d	endpoint singularities	off-centre peaks
integrate_1d_gauss_kronrod (new)	Boost gauss_kronrod<double, 21>	smooth integrands, off-centre peaks	endpoint singularities

integrate_1d_double_exponential is a behavioural superset of
integrate_1d at default settings: same DE dispatch, same xc
semantics, plus two new optional arguments (absolute_tolerance,
max_refinements) and the variadic-args API (see Public API below).

integrate_1d_gauss_kronrod is a new function with the same
variadic-args API as integrate_1d_double_exponential, an adaptive
Gauss-Kronrod backend (K21, fixed order), and the same optional
arguments (absolute_tolerance, max_depth).

My use case is integrals of functions which are product of normal
density and a smooth function. Due to the light tails of the normal
density this functions goes to zero in tails and Gauss-Kronrod excels
compare to double exponential. I really needed this for some
experiments I'm running, and abs_tol was also required to get rid of
numerical problems.

Why add these alongside integrate_1d?

• Off-centre peaks. DE concentrates nodes near the integration
  endpoints; integrands with near-zero endpoints and a sharp peak
  in the middle get undersampled. Gauss-Kronrod distributes nodes across the whole
  interval by Legendre weights and picks up such peaks directly.
• QUADPACK-style mixed convergence. Both new functions accept
  an explicit absolute_tolerance argument; the convergence test
  becomes error <= max(rel_tol * L1, abs_tol). With
  absolute_tolerance = 0 (default) this reduces to the strict
  pure-relative test of integrate_1d. With abs_tol > 0 the
  user can escape the pathological regime in which the strict
  test is measuring accumulated floating-point round-off against
  itself (e.g. failure mode of nested
  integrate_1d_* in the deep tail of a Gaussian factor).
• Explicit max-iteration knob. Both new functions expose the
  per-class refinement / bisection cap (max_refinements for DE,
  max_depth for GK) as an optional argument. integrate_1d
  currently uses Boost's per-class defaults implicitly with no way
  to override.
• Variadic form. With new function names we can introduce
  variadic arguments following the existing ODE-functions.

Public API

Each new function exposes two overloads in Stan-language, mirroring
the ode_rk45 / ode_rk45_tol convention: a default-tolerance form
and a _tol form that takes the three tolerance / iteration knobs
explicitly. Both forms are variadic: any number of arguments of
any Stan type follow the fixed prefix, and are forwarded to the
integrand functor.

integrate_1d_gauss_kronrod(f, a, b, ...args);
integrate_1d_gauss_kronrod_tol(f, a, b, rel_tol, abs_tol, max_depth, ...args);

integrate_1d_double_exponential(f, a, b, ...args);
integrate_1d_double_exponential_tol(f, a, b, rel_tol, abs_tol, max_refinements, ...args);

C++ user-facing form for both:

template <typename F, typename... Args>
double integrate_1d_{gauss_kronrod|double_exponential}(
    const F& f, double a, double b,
    std::ostream* msgs,
    const Args&... args);

template <typename F, typename... Args>
double integrate_1d_{gauss_kronrod|double_exponential}_tol(
    const F& f, double a, double b,
    double relative_tolerance,
    double absolute_tolerance,
    int max_depth_or_refinements,
    std::ostream* msgs,
    const Args&... args);

The integrand functor signature is variadic, matching the call:

real f(real x, real xc, ...args);

This differs from integrate_1d, whose user functor is forced into a
fixed (x, xc, array[] real theta, array[] real x_r, array[] int x_i)
shape; callers had to pack parameters and data into those three
containers. The variadic form lets users pass real, int, vector,
matrix, array[] of any of these, etc., directly. xc is meaningful
under DE; under GK it is always passed as NaN (Boost's gauss_kronrod
has no distance-to-boundary concept).

Design notes

Gauss-Kronrod

• K21 default Kronrod order. Matches scipy/MATLAB/QUADPACK
  conventions. Polynomial exactness of the K rule at N=21 is
  degree 31, comfortable for Gaussian-times-likelihood shapes
  that dominate Stan use. Higher N would tighten initial node
  spacing on sharply peaked integrands at 1.5x-3x constant cost
  on smooth integrals. The two robustness tools we already have
  (absolute_tolerance and informed bound choice) cover the
  failure modes we have encountered; exposing N is a sensible
  future extension but deliberately out of scope here.
• absolute_tolerance argument. Boost's adaptive
  Gauss-Kronrod accumulates error += error_local across
  bisection leaves, each carrying a 2 * eps * |K_local| floor.
  Accumulated round-off scales as ~2^max_depth * eps * L1. For
  integrals whose L1 falls below ~1e-8, the strict
  pure-relative test measures noise against itself; QUADPACK-style
  mixed convergence is the established fix.
• Deliberate weakness on endpoint singularities. Documented
  in tests (a endpoint_singularity_throws test asserts that GK
  throws on 1/sqrt(x) and beta-with-small-shapes integrands).
  Users with such integrands keep using DE.

Double-exponential

• absolute_tolerance applies per piece of the zero-crossing
  split. The existing integrate_1d worker splits integrals
  that cross zero into two pieces (per Boost's exp_sinh docs); the
  new convergence test applies independently to each piece. This
  is the simplest interpretation of "abs_tol is the absolute error
  floor we accept" and the right behaviour when one piece of the
  split has near-zero L1.
• max_refinements default 15. Matches Boost's tanh_sinh
  default and integrate_1d_gauss_kronrod's max_depth for
  symmetry. Boost's exp_sinh/sinh_sinh defaults are 9; the
  unified default of 15 here is mildly more conservative on
  infinite-interval cases.
• xc semantics preserved. DE computes a meaningful xc
  (distance to nearest boundary) and passes it to the user
  functor exactly as integrate_1d does. User code that
  exploits xc carries over unchanged.

Practical guidance

Integrand class	Recommended
Smooth on bounded interval	either (GK marginally faster)
Smooth on ±∞ tail	either
Algebraic/log endpoint singularity	integrate_1d_double_exponential
Gradient with endpoint log-singularity (Beta/Gamma/Weibull lpdf gradients)	integrate_1d_double_exponential
Off-centre peak with near-zero endpoints	integrate_1d_gauss_kronrod
Nested integration with deep-tail evaluations	either, with abs_tol > 0
Modest oscillation	either

Implementation layout

stan/math/prim/functor/integrate_1d_gauss_kronrod.hpp
stan/math/prim/functor/integrate_1d_double_exponential.hpp
stan/math/rev/functor/integrate_1d_gauss_kronrod.hpp
stan/math/rev/functor/integrate_1d_double_exponential.hpp
stan/math/fwd/functor/integrate_1d_gauss_kronrod.hpp
stan/math/fwd/functor/integrate_1d_double_exponential.hpp
stan/math/{prim,rev,fwd}/functor.hpp                # one #include each, two new
test/unit/math/{prim,rev,mix}/functor/              # 71 new tests total, all passing

stanc3 changes will be in a companion PR:

Testing

• 71 new unit tests across prim (15), rev (54), and mix (2),
  all passing. Tests for the new DE function are 1:1 clones of
  integrate_1d's tests with the function name renamed (the new
  function is a strict superset at abs_tol = 0); tests for GK
  are adapted from the same baseline with the endpoint-singular
  cases marked as expected-throw and the gradient-endpoint-
  log-singular PDFs (Beta, ChiSquare, Gamma, Weibull) omitted.
• End-to-end cmdstan smoke tests for both functions exercise
  both overloads on integrals with known closed forms. The GK
  model uses mu as a parameters-block variable to exercise the
  rev autodiff specialisation across 200 MCMC draws (returns
  exactly 1 every time).
• Real-world stress test for integrate_1d_gauss_kronrod: the
  nested random-effects marginal likelihood from a Student-t hierarchical
  linear model (144 observations, 4 chains x 100 draws) produces ELPD
  estimates that agree with independent bridge sampling to
  0.05 nats per fold at the 99th percentile. The existing integrate_1d was
  silently returning zeros.

Companion PRs

• stanc3: function signatures Add new variadic integrate_1d_* signatures stanc3#1626
• stan-dev/docs: function reference entries for both new
  functions (will be made after this PR has been approved).

Release notes

Two new explicit-name Boost 1D adaptive quadrature functions integrate_1d_double_exponential and integrate_1d_gauss_kronrod with three tolerance / max iterations arguments and variadic form.

Checklist

• Copyright holder: (fill in copyright holder information)

Stan Development Team. The code is duplication of existing integrate_1d code and there is no additional creativity.

The copyright holder is typically you or your assignee, such as a university or company. By submitting this pull request, the copyright holder is agreeing to the license the submitted work under the following licenses:
  - Code: BSD 3-clause (https://opensource.org/licenses/BSD-3-Clause)
  - Documentation: CC-BY 4.0 (https://creativecommons.org/licenses/by/4.0/)

• 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

AI Use Disclosure

The duplication of the existing integrate_1d code and tests were made using Claude AI Agent. The actual quadrature algorithm is in Boost library, and the code is just a wrapper to call Boost. Variadic code changes by @WardBrian, and test file function calls modified to follow new API by Claude.

[WardBrian]

Would it make sense to add a suffixed alias of the existing integrate_1d as part of this? I think we have done similar in the past when adding e.g. a second kind of algebra solver

[avehtari]

If we would add integrate_1d_double_exponential I would then also add abs_tol argument and it would not be just an alias.

[WardBrian]

That would also be good, then.

If I recall correctly these functions are also already ready to be variadic at the math level but just aren’t in the language, should we make these new versions variadic from the start?

[avehtari]

I' m confident I can get Claude to create integrate_1d_double_exponential with added abs_tol argument, but I'm not confident I could trust it it to do the change to variadic arguments reliably. variadic arguments would be nice as there is a lot of packing and unpacking going on now, like this

    int K = x_i[1];
    real nu     = theta[1];
    real sigma  = theta[2];
    real sd1_1  = theta[3];
    real sd1_2  = theta[4];
    real rho    = theta[5];
    real z1     = theta[6 + 3 * K];

[WardBrian]

In the existing code (and this PR, looking over it briefly) the _impl function already is variadic. So we would just rename that to drop the _impl, remove the adaptor struct wrapping F, and we should be good to go.

I can help with this (and the stanc stuff, since it gets slightly more complicated to add a variadic), but I think it’s worth doing if we’re introducing new names anyway

[avehtari]

I can add integrate_1d_double_exponential with the additional arguments, and the test also if adding variadic arguments is easy. Should include them both to this PR?

[avehtari]

Ah, didn't see your latest comment. I can add integrate_1d_double_exponential with additional arguments and then you can add variadic

[WardBrian]

Sounds good. As long as claude follows the lead on the existing code, the extra steps to make it variadic should be very easy. The main reason it hasn’t been so far is that we need to come up with a new name for that version, which you’re already doing!

[avehtari]

Ok, added integrate_1d_double_exponential

[WardBrian]

Ok. I will try to review tomorrow, but if not it will take a few days due to the long holiday weekend in the US.

If changes are necessary for the signatures to be variadic in Stanc, would you like me to comment on the diff or just go ahead and make them?

[avehtari]

Just go ahead and make them so you can test right away

[WardBrian]

@avehtari I updated the signatures and posted stan-dev/stanc3#1626. Do you mind updating the tests for the new signatures (variadic rather than forcing the use of the 3 arrays)? I'm guessing claude would be decent at it

[WardBrian]

Is signature compatibility really worth a useless argument? I guess it is nice for people who are switching back and forth?

[avehtari]

• Fixed two pre-existing developer bugs in the fwd headers: integrate_1d_*_impl references that needed renaming to _tol
• Fixed stan_math tests
• 71 unit tests + 8 cmdstan integrals, OK
• did not yet update the PR text (too late now)

[stan-buildbot]

---

Name	Old Result	New Result	Ratio	Performance change( 1 - new / old )
gp_regr/gp_regr.stan	0.1	0.1	1.03	3.11% faster
gp_regr/gen_gp_data.stan	0.03	0.03	1.0	0.39% faster
arK/arK.stan	1.87	1.81	1.03	3.24% faster
eight_schools/eight_schools.stan	0.06	0.06	1.04	3.83% faster
low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	9.26	9.06	1.02	2.17% faster
pkpd/one_comp_mm_elim_abs.stan	20.71	19.94	1.04	3.71% faster
pkpd/sim_one_comp_mm_elim_abs.stan	0.27	0.26	1.02	2.43% faster
sir/sir.stan	77.52	73.64	1.05	5.01% faster
gp_pois_regr/gp_pois_regr.stan	3.19	2.89	1.1	9.27% faster
low_dim_gauss_mix/low_dim_gauss_mix.stan	2.87	2.81	1.02	2.0% faster
irt_2pl/irt_2pl.stan	4.76	4.58	1.04	3.66% faster
arma/arma.stan	0.32	0.31	1.05	4.82% faster
garch/garch.stan	0.48	0.47	1.02	2.4% faster
low_dim_corr_gauss/low_dim_corr_gauss.stan	0.01	0.01	1.0	0.01% faster
performance.compilation	240.31	239.1	1.01	0.5% faster
Mean result: 1.0325789349696564

---

Jenkins Console Log
Blue Ocean
Commit hash: 17c32cddcf8ccbd9fb3ccb00ff530f90ec16062c

---

Machine information

No LSB modules are available. Distributor ID: Ubuntu Description: Ubuntu 20.04.3 LTS Release: 20.04 Codename: focal

CPU:
Architecture: x86_64
CPU op-mode(s): 32-bit, 64-bit
Byte Order: Little Endian
Address sizes: 46 bits physical, 48 bits virtual
CPU(s): 80
On-line CPU(s) list: 0-79
Thread(s) per core: 2
Core(s) per socket: 20
Socket(s): 2
NUMA node(s): 2
Vendor ID: GenuineIntel
CPU family: 6
Model: 85
Model name: Intel(R) Xeon(R) Gold 6148 CPU @ 2.40GHz
Stepping: 4
CPU MHz: 3395.183
CPU max MHz: 3700.0000
CPU min MHz: 1000.0000
BogoMIPS: 4800.00
Virtualization: VT-x
L1d cache: 1.3 MiB
L1i cache: 1.3 MiB
L2 cache: 40 MiB
L3 cache: 55 MiB
NUMA node0 CPU(s): 0,2,4,6,8,10,12,14,16,18,20,22,24,26,28,30,32,34,36,38,40,42,44,46,48,50,52,54,56,58,60,62,64,66,68,70,72,74,76,78
NUMA node1 CPU(s): 1,3,5,7,9,11,13,15,17,19,21,23,25,27,29,31,33,35,37,39,41,43,45,47,49,51,53,55,57,59,61,63,65,67,69,71,73,75,77,79
Vulnerability Gather data sampling: Mitigation; Microcode
Vulnerability Itlb multihit: KVM: Mitigation: Split huge pages
Vulnerability L1tf: Mitigation; PTE Inversion; VMX conditional cache flushes, SMT vulnerable
Vulnerability Mds: Mitigation; Clear CPU buffers; SMT vulnerable
Vulnerability Meltdown: Mitigation; PTI
Vulnerability Mmio stale data: Mitigation; Clear CPU buffers; SMT vulnerable
Vulnerability Reg file data sampling: Not affected
Vulnerability Retbleed: Mitigation; IBRS
Vulnerability Spec rstack overflow: Not affected
Vulnerability Spec store bypass: Mitigation; Speculative Store Bypass disabled via prctl
Vulnerability Spectre v1: Mitigation; usercopy/swapgs barriers and __user pointer sanitization
Vulnerability Spectre v2: Mitigation; IBRS; IBPB conditional; STIBP conditional; RSB filling; PBRSB-eIBRS Not affected; BHI Not affected
Vulnerability Srbds: Not affected
Vulnerability Tsx async abort: Mitigation; Clear CPU buffers; SMT vulnerable
Vulnerability Vmscape: Mitigation; IBPB before exit to userspace
Flags: fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush dts acpi mmx fxsr sse sse2 ss ht tm pbe syscall nx pdpe1gb rdtscp lm constant_tsc art arch_perfmon pebs bts rep_good nopl xtopology nonstop_tsc cpuid aperfmperf pni pclmulqdq dtes64 monitor ds_cpl vmx smx est tm2 ssse3 sdbg fma cx16 xtpr pdcm pcid dca sse4_1 sse4_2 x2apic movbe popcnt tsc_deadline_timer aes xsave avx f16c rdrand lahf_lm abm 3dnowprefetch cpuid_fault epb cat_l3 cdp_l3 invpcid_single pti intel_ppin ssbd mba ibrs ibpb stibp tpr_shadow vnmi flexpriority ept vpid ept_ad fsgsbase tsc_adjust bmi1 hle avx2 smep bmi2 erms invpcid rtm cqm mpx rdt_a avx512f avx512dq rdseed adx smap clflushopt clwb intel_pt avx512cd avx512bw avx512vl xsaveopt xsavec xgetbv1 xsaves cqm_llc cqm_occup_llc cqm_mbm_total cqm_mbm_local dtherm ida arat pln pts hwp hwp_act_window hwp_epp hwp_pkg_req pku ospke md_clear flush_l1d arch_capabilities

G++:
g++ (Ubuntu 9.4.0-1ubuntu1~20.04) 9.4.0
Copyright (C) 2019 Free Software Foundation, Inc.
This is free software; see the source for copying conditions. There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

Clang:
clang version 10.0.0-4ubuntu1
Target: x86_64-pc-linux-gnu
Thread model: posix
InstalledDir: /usr/bin

[avehtari]

I updated the PR text to mention the variadic API