From ad5555c3f465a960f3a12f85f3b41c545c9770e1 Mon Sep 17 00:00:00 2001
From: ChrisRackauckas-Claude <accounts@chrisrackauckas.com>
Date: Thu, 9 Jul 2026 05:49:23 -0400
Subject: [PATCH] Document public Integrals APIs

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
---
 docs/src/solvers/IntegralSolvers.md       |   3 +
 docs/src/tutorials/infinity_transforms.md |   6 +
 src/Integrals.jl                          |  53 +++++-
 src/algorithms.jl                         | 218 +++++++++++++++++-----
 src/algorithms_extension.jl               | 218 +++++++++++++++++++++-
 src/algorithms_meta.jl                    |  27 ++-
 src/algorithms_sampled.jl                 |  50 +++--
 src/common.jl                             |  66 +++++--
 src/infinity_handling.jl                  |  30 +++
 test/qa/public_api_docs.jl                |  46 +++++
 test/qa/qa.jl                             |   2 +
 11 files changed, 638 insertions(+), 81 deletions(-)
 create mode 100644 test/qa/public_api_docs.jl

diff --git a/docs/src/solvers/IntegralSolvers.md b/docs/src/solvers/IntegralSolvers.md
index 8178f52e..6e1bd96e 100644
--- a/docs/src/solvers/IntegralSolvers.md
+++ b/docs/src/solvers/IntegralSolvers.md
@@ -31,4 +31,7 @@ CubaCuhre
 GaussLegendre
 QuadratureRule
 ArblibJL
+HAdaptiveIntegrationJL
+FastTanhSinhQuadratureJL
+ChangeOfVariables
 ```
diff --git a/docs/src/tutorials/infinity_transforms.md b/docs/src/tutorials/infinity_transforms.md
index a3ad2932..f3ba186c 100644
--- a/docs/src/tutorials/infinity_transforms.md
+++ b/docs/src/tutorials/infinity_transforms.md
@@ -23,6 +23,12 @@ with specific decay behavior.
 
 Integrals.jl provides three built-in transformations for handling infinite bounds:
 
+```@docs
+transformation_if_inf
+transformation_tan_inf
+transformation_cot_inf
+```
+
 ### Default: `transformation_if_inf`
 
 The default transformation uses rational functions to map infinite domains to finite intervals:
diff --git a/src/Integrals.jl b/src/Integrals.jl
index 107b0563..88800e0e 100644
--- a/src/Integrals.jl
+++ b/src/Integrals.jl
@@ -31,6 +31,12 @@ include("simpsons.jl")
     QuadSensitivityAlg
 
 Abstract type for quadrature sensitivity algorithms.
+
+## Interface
+
+Subtypes are used as `sensealg` values in `init` and `solve` for `IntegralProblem`s.
+Concrete algorithms must provide or select a vector-Jacobian product path compatible
+with the integrand calling convention and the active automatic differentiation backend.
 """
 abstract type QuadSensitivityAlg end
 """
@@ -38,9 +44,26 @@ abstract type QuadSensitivityAlg end
 
 Wrapper for custom vector-Jacobian product functions in automatic differentiation.
 
+## Arguments
+
+  - `vjp`: Vector-Jacobian product strategy or callable object.
+
 # Fields
 
-  - `vjp::V`: The vector-Jacobian product function
+  - `vjp::V`: Stored vector-Jacobian product strategy.
+
+## Returns
+
+Returns a `ReCallVJP` sensitivity wrapper for the `sensealg` keyword.
+
+## Example
+
+```julia
+using Integrals
+
+prob = IntegralProblem((x, p) -> p[1] * x, (0.0, 1.0), [2.0])
+cache = init(prob, QuadGKJL(); sensealg = Integrals.ReCallVJP(Integrals.ZygoteVJP()))
+```
 """
 struct ReCallVJP{V}
     vjp::V
@@ -51,12 +74,23 @@ end
 
 Abstract type for vector-Jacobian product (VJP) methods used in automatic differentiation
 of integrals.
+
+## Interface
+
+Subtypes identify an automatic differentiation backend. Extension methods implement the
+backend-specific `_compute_dfdp_and_f` dispatch used during differentiated solves. The
+method must return both the parameter derivative information and the integrand wrapper
+needed by the quadrature solve.
 """
 abstract type IntegralVJP end
 """
     ZygoteVJP <: IntegralVJP
 
 Uses Zygote.jl for vector-Jacobian products in automatic differentiation of integrals.
+
+## Returns
+
+Returns a `ZygoteVJP` backend marker for `ReCallVJP`.
 """
 struct ZygoteVJP <: IntegralVJP end
 """
@@ -66,7 +100,11 @@ Uses ReverseDiff.jl for vector-Jacobian products in automatic differentiation of
 
 # Fields
 
-  - `compile::Bool`: Whether to compile the tape for better performance
+  - `compile::Bool`: Whether to compile the ReverseDiff tape.
+
+## Returns
+
+Returns a `ReverseDiffVJP` backend marker for `ReCallVJP`.
 """
 struct ReverseDiffVJP <: IntegralVJP
     compile::Bool
@@ -76,6 +114,10 @@ end
     MooncakeVJP <: IntegralVJP
 
 Uses Mooncake.jl for vector-Jacobian products in automatic differentiation of integrals.
+
+## Returns
+
+Returns a `MooncakeVJP` backend marker for `ReCallVJP`.
 """
 struct MooncakeVJP <: IntegralVJP end
 
@@ -111,7 +153,12 @@ Exception thrown when unrecognized keyword arguments are passed to `solve`.
 
 # Fields
 
-  - `kwargs`: The keyword arguments that were passed
+  - `kwargs`: Keyword arguments that were passed to `solve`.
+
+## Returns
+
+Constructs an exception object. `showerror` reports the allowed common solver keywords
+and the unrecognized names.
 """
 struct CommonKwargError <: Exception
     kwargs::Any
diff --git a/src/algorithms.jl b/src/algorithms.jl
index 6a766f15..72c394e4 100644
--- a/src/algorithms.jl
+++ b/src/algorithms.jl
@@ -1,14 +1,38 @@
 """
     QuadGKJL(; order = 7, norm = norm, buffer = nothing)
 
-One-dimensional Gauss-Kronrod integration from QuadGK.jl.
-This method also takes the optional arguments `order` and `norm`.
-Which are the order of the integration rule
-and the norm for calculating the error, respectively.
-Lastly, the `buffer` keyword, if set (e.g. `buffer=true`), will allocate a buffer to reuse
-for multiple integrals and may require evaluating the integrand unless an
-`integrand_prototype` is provided. Unlike the `segbuf` keyword to `quadgk`, you do not
-allocate the buffer as this is handled automatically.
+One-dimensional adaptive Gauss-Kronrod integration from QuadGK.jl.
+
+## Keyword Arguments
+
+  - `order`: Order of the embedded Gauss-Kronrod rule. Higher values reduce the number
+    of subintervals for smooth integrands and increase work per subinterval.
+  - `norm`: Function used by QuadGK.jl to turn the integral estimate and error estimate
+    into scalar convergence criteria. This is most useful for array-valued integrands.
+  - `buffer`: If non-`nothing`, allocate and cache a QuadGK segment buffer during
+    `init`. This avoids repeated buffer allocation when `solve!` is called on the same
+    cache. Buffer construction may evaluate the integrand unless the problem supplies an
+    `integrand_prototype`.
+
+## Fields
+
+  - `order::Int`: Stored quadrature rule order.
+  - `norm`: Stored convergence norm.
+  - `buffer`: Stored buffer option.
+
+## Returns
+
+Returns a `QuadGKJL` algorithm object for use with
+`solve(prob::IntegralProblem, alg)`, `init`, and `solve!`.
+
+## Example
+
+```julia
+using Integrals
+
+prob = IntegralProblem((x, p) -> x^2, (0.0, 1.0))
+sol = solve(prob, QuadGKJL(order = 9); reltol = 1e-10)
+```
 
 ## References
 
@@ -34,15 +58,37 @@ QuadGKJL(; order = 7, norm = norm, buffer = nothing) = QuadGKJL(order, norm, buf
 """
     HCubatureJL(; norm=norm, initdiv=1, buffer = nothing)
 
-Multidimensional "h-adaptive" integration from HCubature.jl.
-This method also takes the optional arguments `initdiv` and `norm`.
-Which are the initial number of segments
-each dimension of the integration domain is divided into,
-and the norm for calculating the error, respectively.
-Lastly, the `buffer` keyword, if set (e.g. `buffer=true`), will allocate a buffer to reuse
-for multiple integrals and may require evaluating the integrand unless an
-`integrand_prototype` is provided. Unlike the `buffer` keyword to `hcubature/hquadrature`,
-you do not allocate the buffer as this is handled automatically.
+Multidimensional h-adaptive integration from HCubature.jl.
+
+## Keyword Arguments
+
+  - `initdiv`: Initial number of segments used to divide each integration dimension.
+  - `norm`: Function used by HCubature.jl to turn integral and error estimates into
+    scalar convergence criteria. This is most useful for array-valued integrands.
+  - `buffer`: If non-`nothing`, allocate and cache HCubature work buffers during
+    `init`. The buffer is managed by Integrals.jl; callers should not pass an HCubature
+    buffer directly.
+
+## Fields
+
+  - `initdiv::Int`: Stored initial subdivision count.
+  - `norm`: Stored convergence norm.
+  - `buffer`: Stored buffer option.
+
+## Returns
+
+Returns an `HCubatureJL` algorithm object for `IntegralProblem`s over scalar or vector
+box domains.
+
+## Example
+
+```julia
+using Integrals
+
+f(x, p) = x[1]^2 + x[2]^2
+prob = IntegralProblem(f, ([0.0, 0.0], [1.0, 1.0]))
+sol = solve(prob, HCubatureJL(initdiv = 2))
+```
 
 ## References
 
@@ -73,15 +119,39 @@ end
 
 Multidimensional adaptive Monte Carlo integration from MonteCarloIntegration.jl.
 Importance sampling is used to reduce variance.
-This method also takes three optional arguments `nbins`, `ncalls` and `debug`
-which are the initial number of bins
-each dimension of the integration domain is divided into,
-the number of function calls per iteration of the algorithm,
-and whether debug info should be printed, respectively.
+
+## Keyword Arguments
+
+  - `nbins`: Initial number of bins used for each integration dimension.
+  - `ncalls`: Number of integrand calls requested per Monte Carlo iteration.
+  - `debug`: Whether to request debug output from MonteCarloIntegration.jl.
+  - `seed`: Optional random seed passed to the underlying algorithm.
+
+## Fields
+
+  - `nbins::Int`: Stored bin count.
+  - `ncalls::Int`: Stored calls-per-iteration target.
+  - `debug::Bool`: Stored debug-output flag.
+  - `seed`: Stored seed value.
+
+## Returns
+
+Returns a `VEGAS` algorithm object for multidimensional scalar Monte Carlo
+integration.
+
+## Example
+
+```julia
+using Integrals
+
+f(x, p) = exp(-sum(abs2, x))
+prob = IntegralProblem(f, (zeros(3), ones(3)))
+sol = solve(prob, VEGAS(ncalls = 2_000); reltol = 1e-2)
+```
 
 ## Limitations
 
-This algorithm can only integrate `Float64`-valued functions
+This algorithm can only integrate scalar `Float64`-valued functions.
 
 ## References
 
@@ -109,22 +179,52 @@ function VEGAS(; nbins = 100, ncalls = 1000, debug = false, seed = nothing)
 end
 
 """
-    GaussLegendre{C, N, W}
-
-Struct for evaluating an integral via (composite) Gauss-Legendre quadrature.
-The field `C` will be `true` if `subintervals > 1`, and `false` otherwise.
-
-The fields `nodes::N` and `weights::W` are defined by
-`nodes, weights = gausslegendre(n)` for a given number of nodes `n`.
-
-The field `subintervals::Int64 = 1` (with default value `1`) defines the
-number of intervals to partition the original interval of integration
-`[a, b]` into, splitting it into `[xⱼ, xⱼ₊₁]` for `j = 1,…,subintervals`,
-where `xⱼ = a + (j-1)h` and `h = (b-a)/subintervals`. Gauss-Legendre
-quadrature is then applied on each subinterval. For example, if
-`[a, b] = [-1, 1]` and `subintervals = 2`, then Gauss-Legendre
-quadrature will be applied separately on `[-1, 0]` and `[0, 1]`,
-summing the two results.
+    GaussLegendre(; n = 250, subintervals = 1, nodes = nothing, weights = nothing)
+    GaussLegendre(nodes, weights, subintervals = 1)
+
+Fixed-order Gauss-Legendre quadrature, optionally applied on a composite partition of
+the integration interval.
+
+## Arguments
+
+  - `nodes`: Quadrature nodes on the standard interval `[-1, 1]`.
+  - `weights`: Quadrature weights corresponding to `nodes`.
+  - `subintervals`: Number of equally sized subintervals used for composite
+    Gauss-Legendre quadrature.
+
+## Keyword Arguments
+
+  - `n`: Number of quadrature nodes to construct when `nodes` or `weights` is
+    `nothing`.
+  - `subintervals`: Number of subintervals used for composite quadrature. Must be
+    positive.
+  - `nodes`: Optional precomputed nodes. If omitted, `gausslegendre(n)` is used.
+  - `weights`: Optional precomputed weights. If omitted, `gausslegendre(n)` is used.
+
+## Fields
+
+  - `nodes`: Stored quadrature nodes.
+  - `weights`: Stored quadrature weights.
+  - `subintervals::Int64`: Stored number of subintervals.
+
+The type parameter `C` is `true` when `subintervals > 1` and `false` otherwise.
+Composite mode splits `[a, b]` into `subintervals` pieces and applies the same rule to
+each piece.
+
+## Returns
+
+Returns a `GaussLegendre` algorithm object. `GaussLegendre(; n=...)` requires
+FastGaussQuadrature.jl to be loaded so that the `gausslegendre` extension method is
+available.
+
+## Example
+
+```julia
+using Integrals, FastGaussQuadrature
+
+prob = IntegralProblem((x, p) -> cos(x), (0.0, pi / 2))
+sol = solve(prob, GaussLegendre(n = 64))
+```
 """
 struct GaussLegendre{C, N, W} <: SciMLBase.AbstractIntegralAlgorithm
     nodes::N
@@ -151,12 +251,40 @@ end
 """
     QuadratureRule(q; n=250)
 
-Algorithm to construct and evaluate a quadrature rule `q` of `n` points computed from the
-inputs as `x, w = q(n)`. It assumes the nodes and weights are for the standard interval
-`[-1, 1]^d` in `d` dimensions, and rescales the nodes to the specific hypercube being
-solved. The nodes `x` may be scalars in 1d or vectors in arbitrary dimensions, and the
-weights `w` must be scalar. The algorithm computes the quadrature rule `sum(w .* f.(x))` and
-the caller must check that the result is converged with respect to `n`.
+Evaluate a user-supplied fixed quadrature rule.
+
+The rule function `q` must support `nodes, weights = q(n)` and return nodes and weights
+for the standard interval or hypercube `[-1, 1]^d`. Integrals.jl rescales the nodes to
+the problem domain before evaluating the integrand. Nodes may be scalars in one
+dimension or vectors in multiple dimensions; weights must be scalar.
+
+## Arguments
+
+  - `q`: Function that returns `(nodes, weights)` for a requested node count.
+
+## Keyword Arguments
+
+  - `n`: Number of quadrature nodes requested from `q`. Must be positive.
+
+## Fields
+
+  - `q`: Stored quadrature rule constructor.
+  - `n::Int`: Stored quadrature node count.
+
+## Returns
+
+Returns a `QuadratureRule` algorithm object. The method computes the fixed quadrature
+sum and reports success; callers are responsible for checking convergence by changing
+`n` or otherwise validating the chosen rule.
+
+## Example
+
+```julia
+using Integrals, FastGaussQuadrature
+
+prob = IntegralProblem((x, p) -> x^4, (-1.0, 1.0))
+sol = solve(prob, QuadratureRule(gausslegendre; n = 8))
+```
 """
 struct QuadratureRule{Q} <: SciMLBase.AbstractIntegralAlgorithm
     q::Q
diff --git a/src/algorithms_extension.jl b/src/algorithms_extension.jl
index 5231e479..acc417e0 100644
--- a/src/algorithms_extension.jl
+++ b/src/algorithms_extension.jl
@@ -4,12 +4,25 @@
     AbstractIntegralExtensionAlgorithm <: SciMLBase.AbstractIntegralAlgorithm
 
 Abstract type for integration algorithms provided through package extensions.
+
+## Interface
+
+Concrete subtypes are lightweight algorithm configuration objects. The package extension
+that owns the backend must implement `Integrals.__solvebp_call(cache, alg, sensealg,
+domain, p; kwargs...)` or another solver layer used by `solve!`. Constructors may check
+that the required extension is loaded and should store all backend options in fields.
 """
 abstract type AbstractIntegralExtensionAlgorithm <: SciMLBase.AbstractIntegralAlgorithm end
 """
     AbstractIntegralCExtensionAlgorithm <: AbstractIntegralExtensionAlgorithm
 
 Abstract type for integration algorithms that use C or C++ libraries through package extensions.
+
+## Interface
+
+Subtypes follow [`AbstractIntegralExtensionAlgorithm`](@ref). Backend wrappers are
+responsible for converting Julia integrands, domains, tolerances, and iteration limits to
+the C-compatible callback and option formats required by the wrapped library.
 """
 abstract type AbstractIntegralCExtensionAlgorithm <: AbstractIntegralExtensionAlgorithm end
 
@@ -17,15 +30,50 @@ abstract type AbstractIntegralCExtensionAlgorithm <: AbstractIntegralExtensionAl
     AbstractCubaAlgorithm <: AbstractIntegralCExtensionAlgorithm
 
 Abstract type for integration algorithms from the Cuba.jl package.
+
+## Interface
+
+Concrete Cuba algorithms store Cuba option fields and require `using Cuba` before
+construction. Their extension methods must accept multidimensional `IntegralProblem`s,
+forward common `solve` tolerances and iteration limits, and return a SciMLBase integral
+solution with the backend's estimate and residual/error estimate.
 """
 abstract type AbstractCubaAlgorithm <: AbstractIntegralCExtensionAlgorithm end
 
 """
-    CubaVegas()
+    CubaVegas(; flags = 0, seed = 0, minevals = 0, nstart = 1000,
+        nincrease = 500, gridno = 0)
 
 Multidimensional adaptive Monte Carlo integration from Cuba.jl.
 Importance sampling is used to reduce variance.
 
+## Keyword Arguments
+
+  - `flags`: Cuba flags bitmask.
+  - `seed`: Random seed passed to Cuba.
+  - `minevals`: Minimum number of integrand evaluations.
+  - `nstart`: Number of evaluations in the first iteration.
+  - `nincrease`: Increase in evaluations for later iterations.
+  - `gridno`: Cuba grid slot used for state reuse.
+
+## Fields
+
+The fields match the keyword arguments: `flags`, `seed`, `minevals`, `nstart`,
+`nincrease`, and `gridno`.
+
+## Returns
+
+Returns a `CubaVegas` algorithm object. `Cuba.jl` must be loaded before construction.
+
+## Example
+
+```julia
+using Integrals, Cuba
+
+prob = IntegralProblem((x, p) -> x[1] * x[2], (zeros(2), ones(2)))
+sol = solve(prob, CubaVegas(nstart = 2_000))
+```
+
 ## References
 
 ```tex
@@ -51,12 +99,31 @@ struct CubaVegas <: AbstractCubaAlgorithm
 end
 
 """
-    CubaSUAVE()
+    CubaSUAVE(; flags = 0, seed = 0, minevals = 0, nnew = 1000,
+        nmin = 2, flatness = 25.0)
 
 Multidimensional adaptive Monte Carlo integration from Cuba.jl.
 Suave stands for subregion-adaptive VEGAS.
 Importance sampling and subdivision are thus used to reduce variance.
 
+## Keyword Arguments
+
+  - `flags`: Cuba flags bitmask.
+  - `seed`: Random seed passed to Cuba.
+  - `minevals`: Minimum number of integrand evaluations.
+  - `nnew`: Number of new samples per subdivision.
+  - `nmin`: Minimum samples required before subdivision.
+  - `flatness`: Flatness parameter used by SUAVE subdivision.
+
+## Fields
+
+The fields match the keyword arguments: `flags`, `seed`, `minevals`, `nnew`, `nmin`,
+and `flatness`.
+
+## Returns
+
+Returns a `CubaSUAVE` algorithm object. `Cuba.jl` must be loaded before construction.
+
 ## References
 
 ```tex
@@ -82,11 +149,36 @@ struct CubaSUAVE{R} <: AbstractCubaAlgorithm where {R <: Real}
 end
 
 """
-    CubaDivonne()
+    CubaDivonne(; flags = 0, seed = 0, minevals = 0, key1 = 47, key2 = 1,
+        key3 = 1, maxpass = 5, border = 0.0, maxchisq = 10.0,
+        mindeviation = 0.25, xgiven = zeros(Cdouble, 0, 0), nextra = 0,
+        peakfinder = C_NULL)
 
 Multidimensional adaptive Monte Carlo integration from Cuba.jl.
 Stratified sampling is used to reduce variance.
 
+## Keyword Arguments
+
+  - `flags`: Cuba flags bitmask.
+  - `seed`: Random seed passed to Cuba.
+  - `minevals`: Minimum number of integrand evaluations.
+  - `key1`, `key2`, `key3`: Divonne rule-selection keys.
+  - `maxpass`: Maximum number of partitioning passes.
+  - `border`: Border width excluded from partitioning.
+  - `maxchisq`: Maximum chi-square value used for consistency checks.
+  - `mindeviation`: Minimum relative deviation used for subdivision.
+  - `xgiven`: Matrix of user-supplied points.
+  - `nextra`: Number of extra points supplied through `peakfinder`.
+  - `peakfinder`: Cuba peak-finder callback pointer.
+
+## Fields
+
+The fields match the keyword arguments.
+
+## Returns
+
+Returns a `CubaDivonne` algorithm object. `Cuba.jl` must be loaded before construction.
+
 ## References
 
 ```tex
@@ -120,10 +212,26 @@ struct CubaDivonne{R1, R2, R3, R4} <:
 end
 
 """
-    CubaCuhre()
+    CubaCuhre(; flags = 0, minevals = 0, key = 0)
 
 Multidimensional h-adaptive integration from Cuba.jl.
 
+## Keyword Arguments
+
+  - `flags`: Cuba flags bitmask.
+  - `minevals`: Minimum number of integrand evaluations.
+  - `key`: Cuba Cuhre rule-selection key.
+
+## Fields
+
+  - `flags::Int`: Stored flags bitmask.
+  - `minevals::Int`: Stored minimum evaluation count.
+  - `key::Int`: Stored rule-selection key.
+
+## Returns
+
+Returns a `CubaCuhre` algorithm object. `Cuba.jl` must be loaded before construction.
+
 ## References
 
 ```tex
@@ -188,17 +296,37 @@ end
     AbstractCubatureJLAlgorithm <: AbstractIntegralCExtensionAlgorithm
 
 Abstract type for integration algorithms from the Cubature.jl package.
+
+## Interface
+
+Subtypes require `using Cubature` before construction and store Cubature.jl option
+values. Extension methods must pass common `solve` tolerances and iteration limits to
+Cubature.jl and return a SciMLBase integral solution.
 """
 abstract type AbstractCubatureJLAlgorithm <: AbstractIntegralCExtensionAlgorithm end
 
 """
-    CubatureJLh(; error_norm=Cubature.INDIVIDUAL)
+    CubatureJLh(; error_norm = Cubature.INDIVIDUAL)
 
 Multidimensional h-adaptive integration from Cubature.jl.
 `error_norm` specifies the convergence criterion  for vector valued integrands.
 Defaults to `Cubature.INDIVIDUAL`, other options are
 `Cubature.PAIRED`, `Cubature.L1`, `Cubature.L2`, or `Cubature.LINF`.
 
+## Keyword Arguments
+
+  - `error_norm`: Cubature.jl error norm used for vector-valued integrands. The
+    constructor default `0` corresponds to `Cubature.INDIVIDUAL`.
+
+## Fields
+
+  - `error_norm::Int32`: Stored Cubature.jl error-norm code.
+
+## Returns
+
+Returns a `CubatureJLh` algorithm object. `Cubature.jl` must be loaded before
+construction.
+
 ## References
 
 ```tex
@@ -224,7 +352,7 @@ function CubatureJLh(; error_norm = 0)
 end
 
 """
-    CubatureJLp(; error_norm=Cubature.INDIVIDUAL)
+    CubatureJLp(; error_norm = Cubature.INDIVIDUAL)
 
 Multidimensional p-adaptive integration from Cubature.jl.
 This method is based on repeatedly doubling the degree of the cubature rules,
@@ -233,6 +361,20 @@ The used cubature rule is a tensor product of Clenshaw–Curtis quadrature rules
 `error_norm` specifies the convergence criterion  for vector valued integrands.
 Defaults to `Cubature.INDIVIDUAL`, other options are
 `Cubature.PAIRED`, `Cubature.L1`, `Cubature.L2`, or `Cubature.LINF`.
+
+## Keyword Arguments
+
+  - `error_norm`: Cubature.jl error norm used for vector-valued integrands. The
+    constructor default `0` corresponds to `Cubature.INDIVIDUAL`.
+
+## Fields
+
+  - `error_norm::Int32`: Stored Cubature.jl error-norm code.
+
+## Returns
+
+Returns a `CubatureJLp` algorithm object. `Cubature.jl` must be loaded before
+construction.
 """
 struct CubatureJLp <: AbstractCubatureJLAlgorithm
     error_norm::Int32
@@ -257,6 +399,22 @@ real- and complex-valued functions with both inplace and out-of-place forms. See
 documentation for additional details the algorithm arguments and on implementing
 high-precision integrands. Additionally, the error estimate is included in the return value
 of the integral, representing a ball.
+
+## Keyword Arguments
+
+  - `check_analytic`: Whether Arblib should check analyticity assumptions.
+  - `take_prec`: Whether to pass precision information through to the integrand.
+  - `warn_on_no_convergence`: Whether to warn when Arblib reports non-convergence.
+  - `opts`: Arblib option pointer or object passed to the backend.
+
+## Fields
+
+The fields match the keyword arguments: `check_analytic`, `take_prec`,
+`warn_on_no_convergence`, and `opts`.
+
+## Returns
+
+Returns an `ArblibJL` algorithm object. `Arblib.jl` must be loaded before construction.
 """
 struct ArblibJL{O} <: AbstractIntegralCExtensionAlgorithm
     check_analytic::Bool
@@ -276,12 +434,25 @@ end
 """
     VEGASMC(; kws...)
 
-Markov-chain based Vegas algorithm from MCIntegration.jl
+Markov-chain based Vegas algorithm from MCIntegration.jl.
 
 Refer to
 [`MCIntegration.integrate`](https://numericaleft.github.io/MCIntegration.jl/dev/lib/montecarlo/#MCIntegration.integrate-Tuple%7BFunction%7D)
 for documentation on the keywords, which are passed directly to the solver with a set of
 defaults that works for conforming integrands.
+
+## Keyword Arguments
+
+  - `kws...`: Backend keyword arguments forwarded to `MCIntegration.integrate`.
+
+## Fields
+
+  - `kws::NamedTuple`: Stored backend keyword arguments.
+
+## Returns
+
+Returns a `VEGASMC` algorithm object. `MCIntegration.jl` must be loaded before
+construction.
 """
 struct VEGASMC{K <: NamedTuple} <: AbstractIntegralExtensionAlgorithm
     kws::K
@@ -304,6 +475,19 @@ This algorithm supports integration over:
 
 Any keyword arguments are passed directly to `HAdaptiveIntegration.integrate`.
 
+## Keyword Arguments
+
+  - `kws...`: Backend keyword arguments forwarded to `HAdaptiveIntegration.integrate`.
+
+## Fields
+
+  - `kws::NamedTuple`: Stored backend keyword arguments.
+
+## Returns
+
+Returns an `HAdaptiveIntegrationJL` algorithm object. `HAdaptiveIntegration.jl` must be
+loaded before construction.
+
 ## Example
 
 ```julia
@@ -340,6 +524,26 @@ especially for integrands with endpoint singularities or infinite derivatives at
   - `atol`: Absolute tolerance for convergence (default: `0.0`)
   - `max_levels`: Maximum number of refinement levels in adaptive integration (default: `10`)
 
+## Fields
+
+  - `rtol`: Stored relative tolerance.
+  - `atol`: Stored absolute tolerance.
+  - `max_levels::Int`: Stored maximum number of refinement levels.
+
+## Returns
+
+Returns a `FastTanhSinhQuadratureJL` algorithm object. FastTanhSinhQuadrature.jl must be
+loaded before construction.
+
+## Example
+
+```julia
+using Integrals, FastTanhSinhQuadrature
+
+prob = IntegralProblem((x, p) -> sqrt(x), (0.0, 1.0))
+sol = solve(prob, FastTanhSinhQuadratureJL(rtol = 1e-10))
+```
+
 ## Limitations
 
   - Only supports 1D, 2D, and 3D integration
diff --git a/src/algorithms_meta.jl b/src/algorithms_meta.jl
index 35561395..b8b9ac9b 100644
--- a/src/algorithms_meta.jl
+++ b/src/algorithms_meta.jl
@@ -3,6 +3,18 @@
 
 Abstract type for meta-algorithms that wrap other integration algorithms,
 typically to apply transformations or preprocessing steps.
+
+## Interface
+
+Concrete subtypes must wrap an underlying `SciMLBase.AbstractIntegralAlgorithm` and
+must implement the same solver lifecycle as ordinary integral algorithms:
+
+  - `init_cacheval(alg, prob)` may allocate reusable cache state.
+  - `__solve(cache, alg, sensealg, domain, p; kwargs...)` or the lower
+    `__solvebp_call` layer must transform the problem and delegate to the wrapped
+    algorithm.
+  - returned solutions must be built for the user's original problem when the
+    transformation is transparent to callers.
 """
 abstract type AbstractIntegralMetaAlgorithm <: SciMLBase.AbstractIntegralAlgorithm end
 
@@ -10,8 +22,9 @@ abstract type AbstractIntegralMetaAlgorithm <: SciMLBase.AbstractIntegralAlgorit
     ChangeOfVariables(fu2gv, alg)
 
 Apply a change of variables from `∫ f(u,p) du` to an equivalent integral `∫ g(v,p) dv` using
-a helper function `fu2gv(f, u_domain) -> (g, v_domain)` where `f` and `g` should be
-integral functions. Acts as a wrapper to algorithm `alg`.
+a helper function `fu2gv(f, u_domain) -> (g, v_domain)`. The transformed integrand `g`
+must obey the same `IntegralFunction` or `BatchIntegralFunction` calling convention as
+`f`.
 
 This meta-algorithm allows users to apply custom or alternative transformations when
 integrating, particularly useful for handling infinite domains where different
@@ -26,6 +39,16 @@ transformations may provide better accuracy for specific integrand types.
   - [`transformation_cot_inf`](@ref): Cotangent transformation for semi-infinite domains
 - `alg`: The underlying integration algorithm to use (e.g., `QuadGKJL()`, `HCubatureJL()`)
 
+## Fields
+
+  - `fu2gv`: Stored transformation function.
+  - `alg`: Stored wrapped algorithm.
+
+## Returns
+
+Returns a `ChangeOfVariables` meta-algorithm. When solved, the result is reported as a
+solution of the original integral problem.
+
 ## Example
 
 ```julia
diff --git a/src/algorithms_sampled.jl b/src/algorithms_sampled.jl
index 565d9687..3d030317 100644
--- a/src/algorithms_sampled.jl
+++ b/src/algorithms_sampled.jl
@@ -3,24 +3,42 @@
 
 Abstract type for integration algorithms that work with sampled data points,
 such as the trapezoidal rule and Simpson's rule.
+
+## Interface
+
+Concrete subtypes are used with `SampledIntegralProblem` and must support
+`find_weights(x, alg)`, where `x` is the problem sampling grid and `alg` is the
+concrete algorithm. The returned weights must:
+
+  - have the same length as the integration axis of the sampled data,
+  - support `iterate`, `length`, `eltype`, `size`, and scalar indexing, and
+  - be valid inputs to `evalrule(data, weights, dim)`.
+
+The generic sampled solver calls `find_weights` during `init` and recomputes weights
+when the cache grid `x` is replaced.
 """
 abstract type AbstractSampledIntegralAlgorithm <: SciMLBase.AbstractIntegralAlgorithm end
 
 """
     TrapezoidalRule
 
-Struct for evaluating an integral via the trapezoidal rule.
+Composite trapezoidal rule for integrating sampled data.
+
+## Returns
+
+Returns a `TrapezoidalRule` algorithm object for
+`solve(prob::SampledIntegralProblem, alg)`.
 
-Example with sampled data:
+## Example
 
 ```julia
 using Integrals
+
 f = x -> x^2
 x = range(0, 1, length = 20)
 y = f.(x)
-problem = SampledIntegralProblem(y, x)
-method = TrapezoidalRule()
-solve(problem, method)
+prob = SampledIntegralProblem(y, x)
+sol = solve(prob, TrapezoidalRule())
 ```
 """
 struct TrapezoidalRule <: AbstractSampledIntegralAlgorithm end
@@ -28,20 +46,26 @@ struct TrapezoidalRule <: AbstractSampledIntegralAlgorithm end
 """
     SimpsonsRule
 
-Struct for evaluating an integral via the Simpson's composite 1/3-3/8
-rule over `AbstractRange`s (evenly spaced points) and
-Simpson's composite 1/3 rule for non-equidistant grids.
+Composite Simpson rule for integrating sampled data.
+
+For evenly spaced `AbstractRange` grids this uses the composite Simpson 1/3 and 3/8
+rules. For non-equidistant grids it uses a composite Simpson 1/3 construction.
+
+## Returns
 
-Example with equidistant data:
+Returns a `SimpsonsRule` algorithm object for
+`solve(prob::SampledIntegralProblem, alg)`.
+
+## Example
 
 ```julia
 using Integrals
+
 f = x -> x^2
-x = range(0, 1, length = 20)
+x = range(0, 1, length = 21)
 y = f.(x)
-problem = SampledIntegralProblem(y, x)
-method = SimpsonsRule()
-solve(problem, method)
+prob = SampledIntegralProblem(y, x)
+sol = solve(prob, SimpsonsRule())
 ```
 """
 struct SimpsonsRule <: AbstractSampledIntegralAlgorithm end
diff --git a/src/common.jl b/src/common.jl
index 751bf8b3..96eb0dc6 100644
--- a/src/common.jl
+++ b/src/common.jl
@@ -116,19 +116,41 @@ function SciMLBase.solve(::IntegralProblem; kwargs...)
 end
 
 """
-```julia
-solve(prob::IntegralProblem, alg::SciMLBase.AbstractIntegralAlgorithm; kwargs...)
-```
+    solve(prob::IntegralProblem, alg::SciMLBase.AbstractIntegralAlgorithm; kwargs...)
+
+Solve an integral problem with the specified integration algorithm.
+
+## Arguments
+
+  - `prob`: Integral problem containing the integrand, domain, parameters, and problem
+    keyword options.
+  - `alg`: Integration algorithm, such as `QuadGKJL()`, `HCubatureJL()`, or another
+    subtype of `SciMLBase.AbstractIntegralAlgorithm`.
 
 ## Keyword Arguments
 
 The arguments to `solve` are common across all of the quadrature methods.
 These common arguments are:
 
-  - `maxiters` (the maximum number of iterations)
-  - `abstol` (absolute tolerance in changes of the objective value)
-  - `reltol` (relative tolerance  in changes of the objective value)
-  - `verbose` (verbosity control via [`IntegralVerbosity`](@ref); defaults to `IntegralVerbosity()`)
+  - `maxiters`: Maximum number of algorithm iterations or backend evaluations.
+  - `abstol`: Absolute tolerance for the integral residual/error estimate.
+  - `reltol`: Relative tolerance for the integral residual/error estimate.
+  - `verbose`: Verbosity control via [`IntegralVerbosity`](@ref). Defaults to
+    `IntegralVerbosity()`.
+
+## Returns
+
+Returns a SciMLBase integral solution whose `u` field is the integral estimate and whose
+`resid` field is the backend residual or error estimate when available.
+
+## Example
+
+```julia
+using Integrals
+
+prob = IntegralProblem((x, p) -> sin(x), (0.0, pi))
+sol = solve(prob, QuadGKJL(); reltol = 1e-10, abstol = 1e-10)
+```
 """
 function SciMLBase.solve(
         prob::IntegralProblem,
@@ -200,13 +222,35 @@ function SciMLBase.init(
 end
 
 """
-```julia
-solve(prob::SampledIntegralProblem, alg::SciMLBase.AbstractIntegralAlgorithm; kwargs...)
-```
+    solve(prob::SampledIntegralProblem, alg::SciMLBase.AbstractIntegralAlgorithm; kwargs...)
+
+Integrate sampled data with a sampled-data quadrature algorithm.
+
+## Arguments
+
+  - `prob`: Sampled integral problem containing sampled values `y`, sampling points
+    `x`, and the integration dimension.
+  - `alg`: Sampled-data integration algorithm, such as `TrapezoidalRule()` or
+    `SimpsonsRule()`.
 
 ## Keyword Arguments
 
-There are no keyword arguments used to solve `SampledIntegralProblem`s
+There are no algorithm-independent keyword arguments used to solve
+`SampledIntegralProblem`s.
+
+## Returns
+
+Returns a SciMLBase integral solution whose `u` field is the sampled integral estimate.
+
+## Example
+
+```julia
+using Integrals
+
+x = range(0, 1, length = 21)
+prob = SampledIntegralProblem(x .^ 2, x)
+sol = solve(prob, SimpsonsRule())
+```
 """
 function SciMLBase.solve(
         prob::SampledIntegralProblem,
diff --git a/src/infinity_handling.jl b/src/infinity_handling.jl
index 74883bf7..07364218 100644
--- a/src/infinity_handling.jl
+++ b/src/infinity_handling.jl
@@ -137,6 +137,16 @@ intervals using the following substitutions:
 This is the default transformation applied by algorithms like `QuadGKJL` and `HCubatureJL`
 when encountering infinite integration bounds.
 
+## Arguments
+
+  - `f`: `IntegralFunction` or `BatchIntegralFunction` to transform.
+  - `domain`: Original integration domain, usually `(lb, ub)`.
+
+## Returns
+
+Returns `(g, tdomain)`, where `g` is the transformed integrand and `tdomain` is the
+finite transformed domain.
+
 ## Example
 
 ```julia
@@ -239,6 +249,16 @@ using the transformation:
 This transformation can provide better accuracy than the default rational transformation
 for some integrands, particularly those that decay like `1/(1+x²)`.
 
+## Arguments
+
+  - `f`: `IntegralFunction` or `BatchIntegralFunction` to transform.
+  - `domain`: Original integration domain, usually `(lb, ub)`.
+
+## Returns
+
+Returns `(g, tdomain)`, where `g` is the transformed integrand and `tdomain` is the
+finite transformed domain.
+
 ## Example
 
 ```julia
@@ -357,6 +377,16 @@ For `(-\\infty, a]`:
 s = -\\cot\\left[\\frac{(\\pi + 2\\arctan(a))(\\xi+1)}{4}\\right], \\quad \\xi \\in [-1, 1]
 ```
 
+## Arguments
+
+  - `f`: `IntegralFunction` or `BatchIntegralFunction` to transform.
+  - `domain`: Original integration domain, usually `(lb, ub)`.
+
+## Returns
+
+Returns `(g, tdomain)`, where `g` is the transformed integrand and `tdomain` is the
+finite transformed domain.
+
 ## Example
 
 ```julia
diff --git a/test/qa/public_api_docs.jl b/test/qa/public_api_docs.jl
new file mode 100644
index 00000000..ba29cf14
--- /dev/null
+++ b/test/qa/public_api_docs.jl
@@ -0,0 +1,46 @@
+using Integrals, Test
+
+function _has_docstring(mod::Module, name::Symbol)
+    return haskey(Docs.meta(mod), Docs.Binding(mod, name))
+end
+
+function _rendered_docs_entries()
+    docs_root = normpath(joinpath(@__DIR__, "..", "..", "docs", "src"))
+    entries = Set{Symbol}()
+    for (root, _, files) in walkdir(docs_root)
+        for file in files
+            endswith(file, ".md") || continue
+            in_docs_block = false
+            for line in eachline(joinpath(root, file))
+                stripped = strip(line)
+                if stripped == "```@docs"
+                    in_docs_block = true
+                    continue
+                elseif startswith(stripped, "```")
+                    in_docs_block = false
+                    continue
+                end
+                in_docs_block || continue
+                isempty(stripped) && continue
+                startswith(stripped, "#") && continue
+                identifier = replace(stripped, r"\(.*" => "")
+                identifier = split(identifier, ".")[end]
+                isempty(identifier) || push!(entries, Symbol(identifier))
+            end
+        end
+    end
+    return entries
+end
+
+@testset "Integrals-owned public API docs" begin
+    public_names = filter(sort!(collect(names(Integrals; all = false)))) do name
+        name !== :Integrals && parentmodule(getfield(Integrals, name)) === Integrals
+    end
+
+    undocumented = filter(name -> !_has_docstring(Integrals, name), public_names)
+    @test isempty(undocumented)
+
+    rendered_entries = _rendered_docs_entries()
+    missing_rendered = filter(name -> name ∉ rendered_entries, public_names)
+    @test isempty(missing_rendered)
+end
diff --git a/test/qa/qa.jl b/test/qa/qa.jl
index a7aa633c..0c11f729 100644
--- a/test/qa/qa.jl
+++ b/test/qa/qa.jl
@@ -1,6 +1,8 @@
 using SciMLTesting, Integrals, Test
 using JET
 
+include("public_api_docs.jl")
+
 run_qa(
     Integrals;
     explicit_imports = true,