Answer

All checks pass and the markdown is fully self-contained (extracting its code blocks into a clean directory and running them gives 57/57). Here is the diagnosis and solution.

Building a Sound Subsampled-RDP Privacy Accountant

Deliverables (all verified in this environment):

Reproduce: python3 test_rdp_accountant.py

Root-cause analysis

The naive approach fails in four independent ways:

1. Floating-point overflow (unsound if swallowed). The Poisson-subsampled Gaussian RDP is a sum of terms containing exp(k(k−1)/(2σ²)). With σ=0.5, k=20 the exponent is 760 > log(DBL_MAX)≈709, so a direct sum(exp(...)) overflows to inf for ordinary parameters. If the overflow is clamped, the reported δ is too small.
2. Invalid amplification (the big one). Fixed-size sampling (without replacement) is not Poisson. For the same q=m/n the without-replacement RDP is strictly larger. Substituting Poisson under-states δ by 36×–239× on a concrete replace-one run — a guarantee that claims far more privacy than exists.
3. Unsound interpolation. Interpolating ε(α) itself need not upper-bound the truth. One must interpolate the convex log-moment g(α) = (α−1)ε(α) = log E[L^α], whose chord lies above the curve.
4. Adjacency / independence. The without-replacement bound is a replace-one result; the Poisson bound is add/remove. Order-wise RDP addition also requires independent subsampling across steps.

The exact fix

• Integer-order RDP computed entirely in log space with a stable log_sum_exp — no overflow.
• Separate formula per sampler: Poisson binomial sum vs. Wang–Balle–Kasiviswanathan forward-difference formula (log-domain) for fixed-size/WOR.
• Accumulate g(α) = (α−1)ε(α) across steps (valid for adaptive composition).
• Interpolate g (convex) between queried Rényi orders, then convert δ = min_α exp((α−1)(ε_α−ε))·(α−1)^(α−1)/α^α.
• Refuse to certify on σ≤0, bad ranges, invalid amplification/adjacency, dependent batching, non-finite RDP.

Key pieces (full source in rdp_accountant.py / SOLUTION.md):

def poisson_subsampled_gaussian_rdp(alpha, q, sigma):      # log-space
    log_terms = [ (lgamma(a+1)-lgamma(k+1)-lgamma(a-k+1))
                  + (a-k)*log1p(-q) + k*log(q)
                  + k*(k-1)/(2*sigma**2) for k in range(a+1) ]
    return log_sum_exp(log_terms) / (alpha - 1)

def fixed_size_subsampled_gaussian_rdp(alpha, q, sigma):   # WOR, alpha int
    return fixed_size_subsampled_gaussian_rdp_int(q, sigma, int(alpha)) / (alpha-1)
    # ... forward differences of exp(cgf) in log space (W-B-K Theorem 27)

def renyi_to_delta(alpha, rdp_alpha, epsilon):             # log-space closed form
    log_delta = (alpha-1)*(rdp_alpha-epsilon) + (alpha-1)*log(alpha-1) - alpha*log(alpha)
    return min(1.0, exp(log_delta)) if log_delta < 0 else 1.0

Usage:

acc = RDPAccountant()
acc.add_step(GaussianMechanism(sigma=1.1, sampling="poisson", sample_rate=0.01), times=1000)
acc.add_step(GaussianMechanism(sigma=1.5, sampling="fixed_size", batch_size=256,
                               dataset_size=50000, adjacency="replace_one"), times=200)
print(acc.sound, acc.get_delta(2.0), acc.get_epsilon(1e-5))

Verification

• Per-order RDP matches Google dp-accounting to 3.6e-16 (Poisson) and 4.4e-16 (WOR); composed ε matches (Poisson 1.711770091 vs 1.711770166; WOR identical to printed precision).
• Conversion is a valid upper bound of the analytic Gaussian δ; the tight binary-divergence conversion is provably ≥ analytic and ≤ the classic bound.
• Interpolation soundness checked against high-precision mpmath quadrature.
• Composition linearity, order-independence, and heterogeneous direct sums checked.
• Unsoundness detected: WOR δ(ε=4)=0.9117 vs buggy Poisson δ=0.0251 (36×), sound=False, message names the failed assumption; naive RDP overflows while log-space yields 124.9567.

Full root-cause discussion, the complete 764-line module, the 394-line test suite, and all captured outputs are in ~/SOLUTION.md.