diff --git a/.github/workflows/README.md b/.github/workflows/README.md
index c58ba2d..d6a54ab 100644
--- a/.github/workflows/README.md
+++ b/.github/workflows/README.md
@@ -8,11 +8,12 @@ This folder contains the repository's GitHub Actions workflows.
 - `ci-full.yml`: broader matrix CI for pull requests to `develop` and `master`.
 - `docs.yml`: Documenter build and deployment workflow.
 - `TagBot.yml`: release-tag automation.
+- `nightly.yml`: scheduled benchmarking on non-Linux platforms.
 
 ### Disabled
 
 - `format.yml.disabled`: disabled formatter workflow kept for reference.
 - `CompatHelper.yml.disabled`: dependency update automation. 
-- `nightly.yml.disabled`: scheduled regression sweep on non-Linux platforms.
+
 
 See [`../../docs/src/ci-testing.md`](../../docs/src/ci-testing.md) for the user-facing workflow overview.
diff --git a/.github/workflows/ci-full.yml b/.github/workflows/ci-full.yml
index ae5236a..a706975 100644
--- a/.github/workflows/ci-full.yml
+++ b/.github/workflows/ci-full.yml
@@ -14,7 +14,7 @@ jobs:
     name: Julia ${{ matrix.julia-version }} – ${{ matrix.os }}
     runs-on: ${{ matrix.os }}
     strategy:
-      fail-fast: false
+      fail-fast: true
       matrix:
         julia-version: ['1.10', '1.11']
         os: [ubuntu-latest, macos-latest, windows-latest]
@@ -46,7 +46,7 @@ jobs:
     runs-on: ${{ matrix.os }}
     needs: test
     strategy:
-      fail-fast: false
+      fail-fast: true
       matrix:
         os: [ubuntu-latest, macos-latest, windows-latest]
 
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 40ecf2d..8d07813 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -15,7 +15,7 @@ jobs:
     name: Julia ${{ matrix.julia-version }} – Linux
     runs-on: ubuntu-latest
     strategy:
-      fail-fast: false
+      fail-fast: true
       matrix:
         julia-version: ['1.12']
 
diff --git a/.github/workflows/nightly.yml.disabled b/.github/workflows/nightly.yml
similarity index 79%
rename from .github/workflows/nightly.yml.disabled
rename to .github/workflows/nightly.yml
index f7f3a71..032c3d5 100644
--- a/.github/workflows/nightly.yml.disabled
+++ b/.github/workflows/nightly.yml
@@ -14,12 +14,10 @@ jobs:
     name: Nightly – Julia ${{ matrix.julia-version }} – ${{ matrix.os }}
     runs-on: ${{ matrix.os }}
     strategy:
-      fail-fast: false
+      fail-fast: true
       matrix:
         include:
-          - os: macos-latest
-            julia-version: '1.12'
-          - os: windows-latest
+          - os: ubuntu-latest
             julia-version: '1.12'
 
     steps:
@@ -41,5 +39,5 @@ jobs:
       - name: Install dependencies
         run: julia --project=. -e 'using Pkg; Pkg.instantiate()'
 
-      - name: Run unit tests
-        run: julia --project=. -e 'using Pkg; Pkg.test()'
+      - name: Run benchmark
+        run: make ci
diff --git a/Makefile b/Makefile
index 3e8c929..50b648e 100644
--- a/Makefile
+++ b/Makefile
@@ -15,7 +15,9 @@ QMCPY_PYTHON_AUTO := $(shell \
 	done)
 PYTHON ?= $(if $(QMCPY_PYTHON_AUTO),$(QMCPY_PYTHON_AUTO),python)
 BENCH_COVERAGE ?= 0
+BENCH_BLAS_THREADS ?= 1
 JULIA_BENCH_COVERAGE_FLAG := $(if $(filter 1,$(BENCH_COVERAGE)),--code-coverage=user,)
+BENCH_THREAD_ENV := QMC_BENCH_BLAS_THREADS=$(BENCH_BLAS_THREADS) OPENBLAS_NUM_THREADS=$(BENCH_BLAS_THREADS) MKL_NUM_THREADS=$(BENCH_BLAS_THREADS) OMP_NUM_THREADS=$(BENCH_BLAS_THREADS) NUMEXPR_NUM_THREADS=$(BENCH_BLAS_THREADS)
 
 define RUN_TIMED
 	@start=$$(date +%s); \
@@ -136,7 +138,7 @@ notebook-%:
 # Override output label with: make bench LABEL=gaus or make bench gaus
 # Saves results to benchmark/results/latest.json (default) or <LABEL>.json
 bench:
-	$(call RUN_TIMED,BENCH_COVERAGE=$(BENCH_COVERAGE) julia $(JULIA_BENCH_COVERAGE_FLAG) benchmark/runbenchmarks.jl $(LABEL),bench)
+	$(call RUN_TIMED,$(BENCH_THREAD_ENV) BENCH_COVERAGE=$(BENCH_COVERAGE) julia $(JULIA_BENCH_COVERAGE_FLAG) benchmark/runbenchmarks.jl $(LABEL),bench)
 
 check-qmcpy-python:
 	@$(PYTHON) -c "import qmcpy" >/dev/null 2>&1 || \
@@ -152,7 +154,7 @@ check-qmcpy-python:
 REV ?= HEAD
 LABEL ?=
 bench-compare:
-	$(call RUN_TIMED,BENCH_COVERAGE=$(BENCH_COVERAGE) julia $(JULIA_BENCH_COVERAGE_FLAG) benchmark/compare.jl $(REV) $(LABEL),bench-compare)
+	$(call RUN_TIMED,$(BENCH_THREAD_ENV) BENCH_COVERAGE=$(BENCH_COVERAGE) julia $(JULIA_BENCH_COVERAGE_FLAG) benchmark/compare.jl $(REV) $(LABEL),bench-compare)
 
 # Side-by-side Julia vs QMCPy comparison.
 # Runs both the Julia and QMCPy benchmark harnesses for the requested labels.
@@ -163,35 +165,40 @@ bench-compare:
 # If LABEL is set, it also becomes the default input label for both sides.
 JL_LABEL ?= $(if $(LABEL),$(LABEL),latest)
 PY_LABEL ?= $(JL_LABEL)
+BENCH_COMPARE_OUT := $(if $(LABEL),benchmark/results/compare_$(LABEL).md,benchmark/results/compare_head.md)
+BENCH_COMPARE_PY_OUT := $(if $(LABEL),benchmark/results/compare_python_$(LABEL).md,benchmark/results/compare_python.md)
 bench-compare-py: bench check-qmcpy-python
-	$(call RUN_TIMED,$(PYTHON) benchmark/benchmark_qmcpy.py $(PY_LABEL) && BENCH_COVERAGE=$(BENCH_COVERAGE) julia $(JULIA_BENCH_COVERAGE_FLAG) benchmark/compare_py.jl $(JL_LABEL) $(PY_LABEL) $(LABEL),bench-compare-py)
+	$(call RUN_TIMED,$(BENCH_THREAD_ENV) $(PYTHON) benchmark/benchmark_qmcpy.py $(PY_LABEL) && $(BENCH_THREAD_ENV) BENCH_COVERAGE=$(BENCH_COVERAGE) julia $(JULIA_BENCH_COVERAGE_FLAG) benchmark/compare_py.jl $(JL_LABEL) $(PY_LABEL) $(LABEL),bench-compare-py)
 
 # Explicit labeled Julia-vs-QMCPy comparison flow.
 # Usage: make bench-compare-py-label LABEL=base
 bench-compare-py-label: check-qmcpy-python
-	$(call RUN_TIMED,BENCH_COVERAGE=$(BENCH_COVERAGE) julia $(JULIA_BENCH_COVERAGE_FLAG) benchmark/runbenchmarks.jl $(LABEL) && $(PYTHON) benchmark/benchmark_qmcpy.py $(LABEL) && BENCH_COVERAGE=$(BENCH_COVERAGE) julia $(JULIA_BENCH_COVERAGE_FLAG) benchmark/compare_py.jl $(LABEL) $(LABEL) $(LABEL),bench-compare-py-label)
+	$(call RUN_TIMED,$(BENCH_THREAD_ENV) BENCH_COVERAGE=$(BENCH_COVERAGE) julia $(JULIA_BENCH_COVERAGE_FLAG) benchmark/runbenchmarks.jl $(LABEL) && $(BENCH_THREAD_ENV) $(PYTHON) benchmark/benchmark_qmcpy.py $(LABEL) && $(BENCH_THREAD_ENV) BENCH_COVERAGE=$(BENCH_COVERAGE) julia $(JULIA_BENCH_COVERAGE_FLAG) benchmark/compare_py.jl $(LABEL) $(LABEL) $(LABEL),bench-compare-py-label)
 
 # Run the labeled Julia-only comparison and Julia-vs-QMCPy comparison in one task.
+# This target does not define a third ratio; inspect:
+#   - $(BENCH_COMPARE_OUT)      with ratio = reference ÷ local
+#   - $(BENCH_COMPARE_PY_OUT)   with ratio = Python ÷ Julia
 # Usage: make bench-all LABEL=base
 bench-all:
-	$(call RUN_TIMED,$(MAKE) bench-compare BENCH_COVERAGE=$(BENCH_COVERAGE) LABEL=$(LABEL) && $(MAKE) bench-compare-py BENCH_COVERAGE=$(BENCH_COVERAGE) LABEL=$(LABEL),bench-all)
+	$(call RUN_TIMED,$(MAKE) bench-compare BENCH_COVERAGE=$(BENCH_COVERAGE) BENCH_BLAS_THREADS=$(BENCH_BLAS_THREADS) LABEL=$(LABEL) && $(MAKE) bench-compare-py BENCH_COVERAGE=$(BENCH_COVERAGE) BENCH_BLAS_THREADS=$(BENCH_BLAS_THREADS) LABEL=$(LABEL) && printf '\n[bench-all] wrote %s (ratio = reference ÷ local) and %s (ratio = Python ÷ Julia)\n' "$(BENCH_COMPARE_OUT)" "$(BENCH_COMPARE_PY_OUT)",bench-all)
 
 # Run the benchmark suite with coverage enabled and produce an lcov report over
 # both src/ and benchmark/ coverage files.
 bench-coverage:
-	$(call RUN_TIMED,find src benchmark -name '*.cov' -delete && rm -f lcov.info && $(MAKE) bench BENCH_COVERAGE=1 && julia --project=. devtools/process_coverage.jl src benchmark && find src benchmark -name '*.cov' -delete,bench-coverage)
+	$(call RUN_TIMED,find src benchmark -name '*.cov' -delete && rm -f lcov.info && $(MAKE) bench BENCH_COVERAGE=1 BENCH_BLAS_THREADS=$(BENCH_BLAS_THREADS) && julia --project=. devtools/process_coverage.jl src benchmark && find src benchmark -name '*.cov' -delete,bench-coverage)
 
 # Run the Julia-vs-Julia comparison with coverage enabled and produce an lcov report.
 bench-compare-coverage:
-	$(call RUN_TIMED,find src benchmark -name '*.cov' -delete && rm -f lcov.info && $(MAKE) bench-compare BENCH_COVERAGE=1 REV=$(REV) LABEL=$(LABEL) && julia --project=. devtools/process_coverage.jl src benchmark && find src benchmark -name '*.cov' -delete,bench-compare-coverage)
+	$(call RUN_TIMED,find src benchmark -name '*.cov' -delete && rm -f lcov.info && $(MAKE) bench-compare BENCH_COVERAGE=1 BENCH_BLAS_THREADS=$(BENCH_BLAS_THREADS) REV=$(REV) LABEL=$(LABEL) && julia --project=. devtools/process_coverage.jl src benchmark && find src benchmark -name '*.cov' -delete,bench-compare-coverage)
 
 # Run the Julia-vs-QMCPy comparison with coverage enabled and produce an lcov report.
 bench-compare-py-coverage:
-	$(call RUN_TIMED,find src benchmark -name '*.cov' -delete && rm -f lcov.info && $(MAKE) bench-compare-py BENCH_COVERAGE=1 LABEL=$(LABEL) JL_LABEL=$(JL_LABEL) PY_LABEL=$(PY_LABEL) && julia --project=. devtools/process_coverage.jl src benchmark && find src benchmark -name '*.cov' -delete,bench-compare-py-coverage)
+	$(call RUN_TIMED,find src benchmark -name '*.cov' -delete && rm -f lcov.info && $(MAKE) bench-compare-py BENCH_COVERAGE=1 BENCH_BLAS_THREADS=$(BENCH_BLAS_THREADS) LABEL=$(LABEL) JL_LABEL=$(JL_LABEL) PY_LABEL=$(PY_LABEL) && julia --project=. devtools/process_coverage.jl src benchmark && find src benchmark -name '*.cov' -delete,bench-compare-py-coverage)
 
 # Run the full labeled benchmark workflow with coverage enabled and produce an lcov report.
 bench-all-coverage:
-	$(call RUN_TIMED,find src benchmark -name '*.cov' -delete && rm -f lcov.info && $(MAKE) bench-all BENCH_COVERAGE=1 LABEL=$(LABEL) && julia --project=. devtools/process_coverage.jl src benchmark && find src benchmark -name '*.cov' -delete,bench-all-coverage)
+	$(call RUN_TIMED,find src benchmark -name '*.cov' -delete && rm -f lcov.info && $(MAKE) bench-all BENCH_COVERAGE=1 BENCH_BLAS_THREADS=$(BENCH_BLAS_THREADS) LABEL=$(LABEL) && julia --project=. devtools/process_coverage.jl src benchmark && find src benchmark -name '*.cov' -delete,bench-all-coverage)
 
 # Compare two saved Julia benchmark-result labels and decide which one is better.
 # Usage: make bench-compare-labels LABEL_A=a LABEL_B=b [OUT_LABEL=report]
@@ -203,7 +210,7 @@ bench-compare-labels:
 # Combination of above targets
 # ============================================================================
 
-# Run the local full-check pipeline: format code, collect unit-test,
+# Run the full-check pipeline: format code, collect unit-test,
 # then collect full benchmark. Pass LABEL=... through to the benchmark step.
-local-ci:
-	$(call RUN_TIMED,$(MAKE) format && $(MAKE) test && $(MAKE) bench-all LABEL=$(LABEL),local-ci)
\ No newline at end of file
+ci:
+	$(call RUN_TIMED,$(MAKE) format && $(MAKE) test && $(MAKE) bench-all LABEL=$(LABEL),local-ci)
diff --git a/README.md b/README.md
index 8a38328..7b42647 100644
--- a/README.md
+++ b/README.md
@@ -25,18 +25,23 @@ println("Exact:    $(keister_exact(3))")
 
 ## Components
 
-QMC.jl provides four building blocks that snap together:
+QMC.jl provides four core building blocks, plus related kernel types:
 
 | Component | Types |
 |---|---|
 | **Discrete Distribution** | `IIDStdUniform`, `Lattice`, `DigitalNetB2`, `Halton`, `Kronecker`, `DigitalNetAnyBases`, `Faure` |
-| **True Measure** | `Uniform`, `Gaussian`, `BrownianMotion`, `Lebesgue`, `GeometricBrownianMotion`, `StudentT`, `Triangular`, `Kumaraswamy`, `JohnsonsSU`, `BernoulliCont`, `AcceptanceRejection`, `DistributionsWrapper`, `MaternGP`, `UniformTriangle`, `ZeroInflatedExpUniform` |
+| **True Measure** | `Uniform`, `Gaussian`, `BrownianMotion`, `Lebesgue`, `GeometricBrownianMotion`, `StudentT`, `Triangular`, `Kumaraswamy`, `JohnsonsSU`, `BernoulliCont`, `AcceptanceRejection`, `AcceptanceRejectionReal`, `DistributionsWrapper`, `MaternGP`, `UniformTriangle`, `ZeroInflatedExpUniform` |
 | **Integrand** | `Keister`, `Genz`, `AsianOption`, `FinancialOption`, `FinancialOptionML`, `BoxIntegral`, `Linear0`, `Sin1D`, `Ishigami`, `Hartmann6D`, `Multimodal2D`, `FourBranch2D`, `SensitivityIndices`, `BayesianLRCoeffs`, `UMBridgeWrapper`, `CustomFun` |
 | **Stopping Criterion** | `CubMCCLT`, `CubMCCLTVec`, `CubMCG`, `CubQMCLatticeG`, `CubQMCNetG`, `CubQMCNetGRep`, `CubQMCBayesLatticeG`, `CubQMCBayesNetG`, `CubQMCRepStudentT`, `CubMLMC`, `CubMLMCCont`, `CubMLQMC`, `CubMLQMCCont`, `PFGPCI` |
-| **Kernel** | `KernelShiftInvar`, `KernelDigShiftInvar`, `KernelMatern12`, `KernelMatern32`, `KernelMatern52`, `KernelGaussian`, `KernelMultiTask`, `SumKernel`, `ProductKernel` |
+| **Kernel** | `KernelShiftInvar`, `KernelDigShiftInvar`, `KernelMatern12`, `KernelMatern32`, `KernelMatern52`, `KernelGaussian`, `KernelRationalQuadratic`, `KernelSquaredExponential`, `KernelMultiTask`, `SumKernel`, `ProductKernel` |
 
 Additional utilities include periodization transforms, iteration diagnostics (`IterationLog`), resume/checkpoint support, Walsh-Hadamard and BRO-FFT transforms, and a LatNetBuilder linker.
 
+Advanced status notes: `PFGPCI` is currently exported as a parity placeholder,
+but `integrate(::PFGPCI)` intentionally errors until a Julia GP backend is
+implemented. Likewise, `gpu_fwht!` is currently a CPU-fallback stub rather than
+a real GPU acceleration path.
+
 ## Installation
 
 ### Prerequisites
diff --git a/benchmark/README.md b/benchmark/README.md
index fe9109b..21104be 100644
--- a/benchmark/README.md
+++ b/benchmark/README.md
@@ -43,6 +43,14 @@ This writes:
 benchmark/results/latest.json
 ```
 
+By default the benchmark targets now pin BLAS-style thread env vars to `1` for
+more stable dense-transform timings on both the Julia and QMCPy sides. Override
+that when needed:
+
+```bash
+make bench BENCH_BLAS_THREADS=4
+```
+
 To save to a labeled file instead:
 
 ```bash
@@ -151,6 +159,12 @@ benchmark/results/qmcpy_base.json
 benchmark/results/compare_python_base.md
 ```
 
+The Julia memory sidecar and the QMCPy JSON now also record lightweight
+benchmark metadata such as generation time, thread configuration, and the
+integrate timing sample/repeat counts. `compare_python*.md` surfaces that
+metadata in its header so stale or mismatched artifact pairings are easier to
+spot.
+
 ### Coverage Caveat
 
 Do not interpret coverage-enabled benchmark runs as representative performance comparisons.
@@ -195,6 +209,19 @@ If one memory ratio is `< 1` and the other is `> 1`, that is normal rather than
 
 If `<label>_memory.json` is missing because the Julia benchmarks were generated before this feature was added, `compare_python*.md` will show Julia RSS delta as `n/a` until that label is rerun with `make bench` or `make bench-compare-py`.
 
+The heaviest adaptive integrate rows are intentionally timed with more samples
+than the leaf transform/evaluate rows (`samples=9` in Julia, `repeat=9` in
+QMCPy) because those rows were the main source of weighted-ratio swings.
+
+The QMCPy `StudentT` transform rows now also use a dedicated higher-stability
+setting (`repeat=21`, `warmup=3`), and the Julia `StudentT` rows use
+`samples=9`. Those settings are recorded in the report header so unusually large
+cross-run changes in that family are easier to audit. The aggregate summary in
+`compare_python*.md` now also breaks out `StudentT` timing separately
+(`all matched rows`, `excluding StudentT`, and `StudentT only`) so the headline
+cross-language ratio is easier to interpret when that transform family dominates
+the total.
+
 The Python harness also mirrors most of the newer Julia-only benchmark rows:
 
 - large-`d` `Gaussian(diag)` and `Gaussian(dense)` transforms
@@ -274,6 +301,14 @@ benchmark/results/qmcpy_base.json
 benchmark/results/compare_python_base.md
 ```
 
+`bench-all` does not create a third combined ratio. It just runs both report
+generators, so you need to read each report with its own convention:
+
+- `compare_<label>.md`: `ratio = reference ÷ local`
+- `compare_python_<label>.md`: `ratio = Python ÷ Julia`
+
+So in both of those reports, `ratio > 1` means the local Julia side is faster.
+
 ## Output Files
 
 Common generated files:
diff --git a/benchmark/benchmark_qmcpy.py b/benchmark/benchmark_qmcpy.py
index aed43da..4d7abc5 100644
--- a/benchmark/benchmark_qmcpy.py
+++ b/benchmark/benchmark_qmcpy.py
@@ -44,6 +44,7 @@
 import tracemalloc
 import warnings
 import statistics
+from datetime import datetime
 from pathlib import Path
 
 import numpy as np
@@ -55,6 +56,17 @@
 LARGE_DIMS = [50, 200]
 LARGE_N = [1024, 4096]
 SEED = 42
+DEFAULT_REPEAT = 7
+INTEGRATE_REPEAT = 9
+STUDENT_T_REPEAT = 21
+STUDENT_T_WARMUP_RUNS = 3
+THREAD_ENV_KEYS = (
+    "QMC_BENCH_BLAS_THREADS",
+    "OPENBLAS_NUM_THREADS",
+    "MKL_NUM_THREADS",
+    "OMP_NUM_THREADS",
+    "NUMEXPR_NUM_THREADS",
+)
 
 
 def current_rss_kib():
@@ -90,11 +102,17 @@ def measure_memory(fn):
     return metrics
 
 
-def bench(make_call, *, repeat=7, warmup=True):
+def bench(make_call, *, repeat=DEFAULT_REPEAT, warmup=True):
     """Return median seconds plus approximate memory metrics per call."""
+    if isinstance(warmup, bool):
+        warmup_runs = 1 if warmup else 0
+    else:
+        warmup_runs = int(warmup)
+        if warmup_runs < 0:
+            raise ValueError("warmup must be bool or non-negative integer")
     fn = make_call()
-    if warmup:
-        fn()  # discard first call (cache/allocation warm-up)
+    for _ in range(warmup_runs):
+        fn()  # discard first call(s) to reduce cache/allocation warm-up noise
     mem = measure_memory(fn)
     # Pick an inner count so each timing sample is long enough to be stable.
     timer = timeit.Timer(fn)
@@ -134,12 +152,25 @@ def record(results, group, name, make_call, **kw):
         print(f"  {name:<45s}  ERR: {type(e).__name__}: {e}")
 
 
+def active_thread_env():
+    return {key: os.environ[key] for key in THREAD_ENV_KEYS if key in os.environ}
+
+
 def main():
     label = sys.argv[1] if len(sys.argv) > 1 else "latest"
     results = {"gen_samples": {}, "transform": {}, "evaluate": {}, "integrate": {}}
 
     print("\n\nQMCPy Benchmarks (qmcpy %s)" % getattr(qp, "__version__", "?"))
     print("=" * 70)
+    thread_env = active_thread_env()
+    if thread_env:
+        print(
+            "Thread env: "
+            + ", ".join(f"{key}={value}" for key, value in sorted(thread_env.items()))
+        )
+    print(
+        f"StudentT config: repeat={STUDENT_T_REPEAT}, warmup_runs={STUDENT_T_WARMUP_RUNS}"
+    )
 
     # 1. Discrete distribution sampling --------------------------------------
     # [C] = qmctoolscl C kernel (not a language comparison)
@@ -191,7 +222,14 @@ def make_student_t(n=n):
             tm = qp.StudentT(dd, loc=np.zeros(dim), shape=np.eye(dim), df=2.0)
             x = dd(n)
             return lambda: tm._transform(x)
-        record(results, "transform", f"StudentT d=10 n={n}", make_student_t)
+        record(
+            results,
+            "transform",
+            f"StudentT d=10 n={n}",
+            make_student_t,
+            repeat=STUDENT_T_REPEAT,
+            warmup=STUDENT_T_WARMUP_RUNS,
+        )
 
         def make_johnsons_su(n=n):
             dim = 10
@@ -336,7 +374,7 @@ def _scalar(x):
 
     for name, make_sc, warm in integrate_cases:
         record(results, "integrate", name, integrate_call(make_sc),
-               repeat=3, warmup=warm)
+               repeat=INTEGRATE_REPEAT, warmup=warm)
         # One-shot accuracy measurement: solution value + tolerances used. Kept
         # separate from the timed runs and never aborts the run on failure.
         try:
@@ -358,6 +396,16 @@ def _scalar(x):
     outfile = resdir / f"qmcpy_{label}.json"
     payload = {
         "qmcpy_version": getattr(qp, "__version__", "?"),
+        "python_version": sys.version.split()[0],
+        "label": label,
+        "generated_at": datetime.now().isoformat(timespec="seconds"),
+        "thread_env": thread_env,
+        "benchmark_config": {
+            "default_repeat": DEFAULT_REPEAT,
+            "integrate_repeat": INTEGRATE_REPEAT,
+            "student_t_repeat": STUDENT_T_REPEAT,
+            "student_t_warmup_runs": STUDENT_T_WARMUP_RUNS,
+        },
         "results": results,
     }
     outfile.write_text(json.dumps(payload, indent=2))
diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl
index 73fae96..0d8a293 100644
--- a/benchmark/benchmarks.jl
+++ b/benchmark/benchmarks.jl
@@ -12,6 +12,7 @@ using BenchmarkTools
 using QMC
 import QMC: Uniform
 using Logging
+using LinearAlgebra
 
 # Silence expected, non-actionable warnings (notably CubMCCLT's single-pass
 # non-convergence notice) for EVERY harness that loads this suite — the standalone
@@ -24,6 +25,20 @@ Logging.disable_logging(Logging.Warn)
 
 const SAMPLES = [256, 1024, 4096, 16384]
 const DIMS = [3, 10]
+const BENCH_SAMPLES = 5
+const INTEGRATE_BENCH_SAMPLES = 9
+const STUDENT_T_BENCH_SAMPLES = 9
+
+const BENCH_BLAS_THREADS = let raw = get(ENV, "QMC_BENCH_BLAS_THREADS", "1")
+    threads = try
+        parse(Int, raw)
+    catch err
+        throw(ArgumentError("QMC_BENCH_BLAS_THREADS must be an integer, got $(repr(raw))"))
+    end
+    threads > 0 || throw(ArgumentError("QMC_BENCH_BLAS_THREADS must be ≥ 1, got $threads"))
+    threads
+end
+BLAS.set_num_threads(BENCH_BLAS_THREADS)
 
 # Large-d cases for the Gaussian transform only (see block 2b). These dimensions
 # are where the diagonal fast path's O(n·d²)→O(n·d) saving becomes visible; the
@@ -35,11 +50,11 @@ const LARGE_N = [1024, 4096]
 
 function bench_gen_samples(dd_constructor, dim, n; kwargs...)
     dd = dd_constructor(dim; seed=42, kwargs...)
-    @benchmarkable gen_samples($dd, $n) evals=3 samples=5
+    @benchmarkable gen_samples($dd, $n) evals=3 samples=BENCH_SAMPLES
 end
 
 function bench_integrate(make_sc; kwargs...)
-    @benchmarkable integrate(sc) evals=1 samples=3 setup=(sc = $make_sc())
+    @benchmarkable integrate(sc) evals=1 samples=INTEGRATE_BENCH_SAMPLES setup=(sc = $make_sc())
 end
 
 # ── Benchmark Groups ─────────────────────────────────────────────────────
@@ -70,7 +85,7 @@ for dim in DIMS, n in SAMPLES
     tm_gauss = Gaussian(dd)
     x = gen_samples(dd, n)
     SUITE["transform"]["Gaussian d=$dim n=$n"] =
-        @benchmarkable transform($tm_gauss, $x) evals=3 samples=5
+        @benchmarkable transform($tm_gauss, $x) evals=3 samples=BENCH_SAMPLES
 end
 
 # 2b. Large-d transform cases — exercise the Gaussian diagonal fast path at
@@ -87,7 +102,7 @@ for dim in LARGE_DIMS, n in LARGE_N
 
     tm_diag = Gaussian(dd)                       # identity Σ ⇒ diagonal A ⇒ fast path
     SUITE["transform"]["Gaussian(diag) d=$dim n=$n"] =
-        @benchmarkable transform($tm_diag, $x) evals=3 samples=5
+        @benchmarkable transform($tm_diag, $x) evals=3 samples=BENCH_SAMPLES
 
     # Dense positive-definite Σ (unit diagonal, 0.5 off-diagonal) ⇒ dense A ⇒ GEMM.
     cov = fill(0.5, dim, dim)
@@ -96,7 +111,7 @@ for dim in LARGE_DIMS, n in LARGE_N
     end
     tm_dense = Gaussian(dd; covariance=cov)
     SUITE["transform"]["Gaussian(dense) d=$dim n=$n"] =
-        @benchmarkable transform($tm_dense, $x) evals=3 samples=5
+        @benchmarkable transform($tm_dense, $x) evals=3 samples=BENCH_SAMPLES
 end
 
 # 3. Integrand evaluation (pure Julia)
@@ -108,13 +123,13 @@ for n in SAMPLES
 
     f_keister = Keister(tm)
     SUITE["evaluate"]["Keister n=$n"] =
-        @benchmarkable evaluate($f_keister, $x) evals=3 samples=5
+        @benchmarkable evaluate($f_keister, $x) evals=3 samples=BENCH_SAMPLES
 
     # NOTE: QMCPy's Genz supports only OSCILLATORY / CORNER PEAK. For a matched
     # cross-language comparison use :oscillatory here (see benchmark_qmcpy.py).
     f_genz = Genz(tm; kind=:oscillatory)
     SUITE["evaluate"]["Genz(oscillatory) n=$n"] =
-        @benchmarkable evaluate($f_genz, $x) evals=3 samples=5
+        @benchmarkable evaluate($f_genz, $x) evals=3 samples=BENCH_SAMPLES
 end
 
 # 3b. Allocation-audit coverage (opt2). These are the transform/evaluate bodies the
@@ -136,18 +151,18 @@ for n in SAMPLES
     tm_j = JohnsonsSU(dd)
 
     SUITE["evaluate"]["BoxIntegral d=10 n=$n"] =
-        @benchmarkable evaluate($f_box, $xu) evals=3 samples=5
+        @benchmarkable evaluate($f_box, $xu) evals=3 samples=BENCH_SAMPLES
     SUITE["evaluate"]["Linear0 d=10 n=$n"] =
-        @benchmarkable evaluate($f_lin, $xu) evals=3 samples=5
+        @benchmarkable evaluate($f_lin, $xu) evals=3 samples=BENCH_SAMPLES
     SUITE["evaluate"]["Genz(gaussian_peak) d=10 n=$n"] =
-        @benchmarkable evaluate($f_gp, $xu) evals=3 samples=5
+        @benchmarkable evaluate($f_gp, $xu) evals=3 samples=BENCH_SAMPLES
     SUITE["evaluate"]["Genz(continuous) d=10 n=$n"] =
-        @benchmarkable evaluate($f_cont, $xu) evals=3 samples=5
+        @benchmarkable evaluate($f_cont, $xu) evals=3 samples=BENCH_SAMPLES
 
     SUITE["transform"]["StudentT d=10 n=$n"] =
-        @benchmarkable transform($tm_t, $xu) evals=3 samples=5
+        @benchmarkable transform($tm_t, $xu) evals=3 samples=STUDENT_T_BENCH_SAMPLES
     SUITE["transform"]["JohnsonsSU d=10 n=$n"] =
-        @benchmarkable transform($tm_j, $xu) evals=3 samples=5
+        @benchmarkable transform($tm_j, $xu) evals=3 samples=BENCH_SAMPLES
 end
 
 # Large-d variants for the reductions / GEMV kinds, where the cache-friendly
@@ -160,11 +175,11 @@ for dim in LARGE_DIMS, n in LARGE_N
     f_lin = Linear0(tm)
     f_gp = Genz(tm; kind=:gaussian_peak)
     SUITE["evaluate"]["BoxIntegral d=$dim n=$n"] =
-        @benchmarkable evaluate($f_box, $xu) evals=3 samples=5
+        @benchmarkable evaluate($f_box, $xu) evals=3 samples=BENCH_SAMPLES
     SUITE["evaluate"]["Linear0 d=$dim n=$n"] =
-        @benchmarkable evaluate($f_lin, $xu) evals=3 samples=5
+        @benchmarkable evaluate($f_lin, $xu) evals=3 samples=BENCH_SAMPLES
     SUITE["evaluate"]["Genz(gaussian_peak) d=$dim n=$n"] =
-        @benchmarkable evaluate($f_gp, $xu) evals=3 samples=5
+        @benchmarkable evaluate($f_gp, $xu) evals=3 samples=BENCH_SAMPLES
 end
 
 # 4. End-to-end integration (mixed: C-backed generators + pure-Julia compute)
diff --git a/benchmark/compare_py.jl b/benchmark/compare_py.jl
index cf03d03..bc2376e 100644
--- a/benchmark/compare_py.jl
+++ b/benchmark/compare_py.jl
@@ -34,6 +34,7 @@ Pkg.resolve()
 Pkg.instantiate()
 
 using BenchmarkTools
+using Dates
 using JSON3
 using Statistics
 using Printf
@@ -42,8 +43,62 @@ resdir = joinpath(@__DIR__, "results")
 compare_py_outfile(label::AbstractString) =
     isempty(label) ? joinpath(resdir, "compare_python.md") :
     joinpath(resdir, "compare_python_$(label).md")
+const ARTIFACT_SKEW_WARNING_SECONDS = 10 * 60
 
 is_c_kernel_row(name::AbstractString) = occursin(r"Lattice|DigitalNetB2|Halton", name)
+is_student_t_row(name::AbstractString) = occursin("StudentT", name)
+
+function maybe_get(obj, key, default=nothing)
+    obj === nothing && return default
+    for candidate in (key, Symbol(key))
+        try
+            haskey(obj, candidate) && return obj[candidate]
+        catch
+        end
+    end
+    return default
+end
+
+function artifact_mtime(path::AbstractString)
+    secs = floor(Int, stat(path).mtime)
+    return Dates.format(Dates.unix2datetime(secs), dateformat"yyyy-mm-ddTHH:MM:SS")
+end
+
+function format_seconds(seconds::Real)
+    total = round(Int, seconds)
+    if total < 60
+        return "$(total) s"
+    elseif total < 3600
+        mins, secs = divrem(total, 60)
+        return "$(mins) m $(secs) s"
+    end
+    hours, rems = divrem(total, 3600)
+    mins, secs = divrem(rems, 60)
+    return "$(hours) h $(mins) m $(secs) s"
+end
+
+function format_thread_env(thread_env)
+    thread_env === nothing && return "n/a"
+    entries = String[]
+    for (key, value) in pairs(thread_env)
+        push!(entries, string(key) * "=" * string(value))
+    end
+    isempty(entries) && return "n/a"
+    return join(sort(entries), ", ")
+end
+
+function wait_for_artifact(
+    path::AbstractString;
+    timeout_s::Real=10.0,
+    poll_interval_s::Real=0.1,
+)
+    deadline = time() + timeout_s
+    while true
+        isfile(path) && return true
+        time() >= deadline && return false
+        sleep(poll_interval_s)
+    end
+end
 
 function lookup_py_entry(py_results, group, name)
     py_group = get(py_results, Symbol(group), nothing)
@@ -89,6 +144,7 @@ function collect_comparison_rows(jl_results, jl_memory_results, py_results)
                                          Float64(py_entry["rss_delta_kib"]) : nothing,
                         py_error=nothing,
                         c_kernel=is_c_kernel_row(name),
+                        student_t=is_student_t_row(name),
                     ),
                 )
             else
@@ -108,6 +164,7 @@ function collect_comparison_rows(jl_results, jl_memory_results, py_results)
                         py_error=py_entry === nothing ? "no Python data" :
                                  string(py_entry["error"]),
                         c_kernel=is_c_kernel_row(name),
+                        student_t=is_student_t_row(name),
                     ),
                 )
             end
@@ -116,8 +173,9 @@ function collect_comparison_rows(jl_results, jl_memory_results, py_results)
     return rows
 end
 
-function summary_metrics(rows)
-    matched = filter(row -> row.py_ms !== nothing, rows)
+function summary_metrics(rows; include=(row -> true))
+    scoped_rows = filter(include, rows)
+    matched = filter(row -> row.py_ms !== nothing, scoped_rows)
     total_jl_ms = sum((row.jl_ms for row in matched); init=0.0)
     total_py_ms = sum((row.py_ms for row in matched); init=0.0)
     peak_rows = filter(row -> row.py_peak_kib !== nothing, matched)
@@ -131,7 +189,7 @@ function summary_metrics(rows)
     total_py_rss_delta_kib = sum((row.py_rss_delta_kib for row in rss_rows); init=0.0)
     return (
         matched=length(matched),
-        total=length(rows),
+        total=length(scoped_rows),
         jl_ms=total_jl_ms,
         py_ms=total_py_ms,
         time_ratio=total_jl_ms > 0 ? total_py_ms / total_jl_ms : NaN,
@@ -142,7 +200,7 @@ function summary_metrics(rows)
         rss_rows=length(rss_rows),
         jl_rss_kib=total_jl_rss_kib,
         py_rss_delta_kib=total_py_rss_delta_kib,
-        rss_ratio=total_jl_rss_kib > 0 ? total_py_rss_delta_kib / total_jl_rss_kib : NaN,
+        rss_ratio=total_jl_rss_kib > 0 ? total_py_rss_delta_kib / total_jl_rss_kib : nothing,
     )
 end
 
@@ -207,12 +265,12 @@ jl_file = joinpath(resdir, "$(jl_label).json")
 jl_mem_file = joinpath(resdir, "$(jl_label)_memory.json")
 py_file = joinpath(resdir, "qmcpy_$(py_label).json")
 
-isfile(jl_file) || error(
-    "Julia results not found: $jl_file\nRun: make bench" *
+wait_for_artifact(jl_file) || error(
+    "Julia results not found after waiting 10s: $jl_file\nRun: make bench" *
     (jl_label == "latest" ? "" : " then julia benchmark/runbenchmarks.jl $jl_label"),
 )
-isfile(py_file) || error(
-    "QMCPy results not found: $py_file\n" *
+wait_for_artifact(py_file) || error(
+    "QMCPy results not found after waiting 10s: $py_file\n" *
     "Run: python benchmark/benchmark_qmcpy.py $py_label",
 )
 
@@ -223,8 +281,26 @@ py_data = JSON3.read(read(py_file, String))
 jl_memory_results = jl_mem_data === nothing ? nothing : jl_mem_data["results"]
 py_results = py_data["results"]
 py_version = get(py_data, "qmcpy_version", "?")
+py_python_version = maybe_get(py_data, "python_version", "?")
+report_generated_at = Dates.format(Dates.now(), dateformat"yyyy-mm-ddTHH:MM:SS")
+jl_generated_at = maybe_get(jl_mem_data, "generated_at", "n/a")
+jl_blas_threads = maybe_get(jl_mem_data, "blas_threads", "n/a")
+jl_config = maybe_get(jl_mem_data, "benchmark_config", nothing)
+jl_integrate_samples = maybe_get(jl_config, "integrate_samples", "n/a")
+jl_student_t_samples = maybe_get(jl_config, "student_t_samples", "n/a")
+py_generated_at = maybe_get(py_data, "generated_at", "n/a")
+py_config = maybe_get(py_data, "benchmark_config", nothing)
+py_integrate_repeat = maybe_get(py_config, "integrate_repeat", "n/a")
+py_student_t_repeat = maybe_get(py_config, "student_t_repeat", "n/a")
+py_student_t_warmup_runs = maybe_get(py_config, "student_t_warmup_runs", "n/a")
+py_thread_env = format_thread_env(maybe_get(py_data, "thread_env", nothing))
+jl_file_mtime = artifact_mtime(jl_file)
+py_file_mtime = artifact_mtime(py_file)
+artifact_skew_seconds = abs(stat(jl_file).mtime - stat(py_file).mtime)
 rows = collect_comparison_rows(jl_results, jl_memory_results, py_results)
 summary = summary_metrics(rows)
+summary_non_student_t = summary_metrics(rows; include=row -> !row.student_t)
+summary_student_t = summary_metrics(rows; include=row -> row.student_t)
 
 jl_sol_file = joinpath(resdir, "$(jl_label)_solutions.json")
 jl_sol_data = isfile(jl_sol_file) ? JSON3.read(read(jl_sol_file, String)) : nothing
@@ -235,9 +311,34 @@ n_flagged = count(r -> r.check !== nothing && r.check.flagged, accuracy_rows)
 println("Julia vs QMCPy benchmark comparison")
 println("  Julia label   : $jl_label")
 println("  QMCPy label   : $py_label  (qmcpy $py_version)")
+println(
+    "  Julia artifact: $(basename(jl_file))  (mtime $jl_file_mtime, generated $jl_generated_at)",
+)
+println(
+    "  QMCPy artifact: $(basename(py_file))  (mtime $py_file_mtime, generated $py_generated_at)",
+)
+println("  Artifact skew : $(format_seconds(artifact_skew_seconds))")
+println(
+    "  Julia config  : BLAS threads=$(jl_blas_threads), integrate samples=$(jl_integrate_samples)",
+)
+println(
+    "  Python config : python $py_python_version, integrate repeat=$(py_integrate_repeat), thread env=$py_thread_env",
+)
+println(
+    "  StudentT config: Julia samples=$(jl_student_t_samples), Python repeat=$(py_student_t_repeat), warmup runs=$(py_student_t_warmup_runs)",
+)
 println()
 println("  NOTE: C-kernel rows [C] (Lattice/DigitalNetB2/Halton gen_samples) use the")
 println("  same qmctoolscl library on both sides and are NOT a language comparison.")
+if summary_student_t.total > 0
+    println("  NOTE: StudentT rows are also summarized separately because they can dominate")
+    println("  the weighted cross-language time ratio.")
+end
+if artifact_skew_seconds > ARTIFACT_SKEW_WARNING_SECONDS
+    println(
+        "  NOTE: input artifacts are more than 10 minutes apart; rerun `make bench-compare-py` if that was not intentional.",
+    )
+end
 println()
 
 header = @sprintf("  %-48s  %11s  %11s  %7s", "benchmark", "Julia (ms)", "Python (ms)", "ratio")
@@ -278,6 +379,24 @@ println(
     summary.matched,
     summary.total
 )
+if summary_student_t.total > 0
+    @printf(
+        "  weighted time ratio excluding StudentT = %.3f  (Python total %.3f ms vs Julia total %.3f ms across %d/%d matched rows)\n",
+        summary_non_student_t.time_ratio,
+        summary_non_student_t.py_ms,
+        summary_non_student_t.jl_ms,
+        summary_non_student_t.matched,
+        summary_non_student_t.total
+    )
+    @printf(
+        "  weighted time ratio StudentT only = %.3f  (Python total %.3f ms vs Julia total %.3f ms across %d/%d matched rows)\n",
+        summary_student_t.time_ratio,
+        summary_student_t.py_ms,
+        summary_student_t.jl_ms,
+        summary_student_t.matched,
+        summary_student_t.total
+    )
+end
 if summary.peak_rows > 0
     @printf(
         "weighted tracemalloc ratio = %.3f  (Python peak total %.1f KiB vs Julia alloc total %.1f KiB across %d/%d rows)\n",
@@ -292,7 +411,7 @@ else
         "weighted tracemalloc ratio = n/a  (QMCPy results do not record Python memory metrics)",
     )
 end
-if summary.rss_rows > 0
+if summary.rss_rows > 0 && summary.rss_ratio !== nothing
     @printf(
         "weighted RSS delta ratio   = %.3f  (Python RSS Δ total %.1f KiB vs Julia RSS Δ total %.1f KiB across %d/%d rows)\n",
         summary.rss_ratio,
@@ -301,6 +420,10 @@ if summary.rss_rows > 0
         summary.rss_rows,
         summary.total
     )
+elseif summary.rss_rows > 0
+    println(
+        "weighted RSS delta ratio   = n/a  (Julia RSS Δ total is 0.0 KiB across matched rows, so the weighted ratio is undefined)",
+    )
 else
     println("weighted RSS delta ratio   = n/a  (missing Julia or QMCPy RSS delta sidecar data)")
 end
@@ -350,7 +473,22 @@ open(outfile, "w") do io
     println(io, "|---|---|---|")
     println(io, "| Julia label | `$(jl_label)` | — |")
     println(io, "| Python label | — | `$(py_label)` |")
+    println(io, "| Julia artifact mtime | `$(jl_file_mtime)` | — |")
+    println(io, "| Python artifact mtime | — | `$(py_file_mtime)` |")
+    println(io, "| Julia generated at | `$(jl_generated_at)` | — |")
+    println(io, "| Python generated at | — | `$(py_generated_at)` |")
+    println(io, "| Julia BLAS threads | `$(jl_blas_threads)` | — |")
+    println(io, "| Python version | — | `$(py_python_version)` |")
     println(io, "| qmcpy version | — | $(py_version) |")
+    println(
+        io,
+        "| integrate timing | `samples=$(jl_integrate_samples)` | `repeat=$(py_integrate_repeat)` |",
+    )
+    println(
+        io,
+        "| StudentT timing | `samples=$(jl_student_t_samples)` | `repeat=$(py_student_t_repeat)`, `warmup=$(py_student_t_warmup_runs)` |",
+    )
+    println(io, "| thread env | — | `$(py_thread_env)` |")
     println(io, "")
     println(io, "**`ratio = Python time ÷ Julia time`**  ")
     println(
@@ -366,6 +504,13 @@ open(outfile, "w") do io
         io,
         "> `qmctoolscl` library on both sides and are **not** a Julia vs Python comparison.",
     )
+    if summary_student_t.total > 0
+        println(
+            io,
+            "> `StudentT` rows are also summarized separately below because they can dominate",
+        )
+        println(io, "> the weighted cross-language time ratio.")
+    end
     println(
         io,
         "> Julia `alloc KiB` is allocated bytes from BenchmarkTools. Julia `RSS Δ` and Python",
@@ -378,45 +523,81 @@ open(outfile, "w") do io
         io,
         "> peak is Python-managed temporary memory. These are related but not interchangeable.",
     )
+    println(
+        io,
+        "> Report generated at `$(report_generated_at)`. Input artifact skew: `$(format_seconds(artifact_skew_seconds))`.",
+    )
+    if artifact_skew_seconds > ARTIFACT_SKEW_WARNING_SECONDS
+        println(
+            io,
+            "> ℹ️ The Julia and QMCPy input artifacts are more than 10 minutes apart. If that was not intentional, rerun `make bench-compare-py` to refresh both sides together.",
+        )
+    end
     println(io, "")
     println(io, "## Aggregate Summary\n")
-    println(io, "| metric | ratio | Julia total | Python total | rows |")
-    println(io, "|:-------|------:|------------:|-------------:|-----:|")
+    println(io, "| metric | scope | ratio | Julia total | Python total | rows |")
+    println(io, "|:-------|:------|------:|------------:|-------------:|-----:|")
     println(
         io,
-        "| matched benchmarks | $(summary.matched)/$(summary.total) | — | — | $(summary.matched) |",
+        "| matched benchmarks | all matched rows | $(summary.matched)/$(summary.total) | — | — | $(summary.matched) |",
     )
     @printf(
         io,
-        "| weighted time ratio | %.3f | %.3f ms | %.3f ms | %d |\n",
+        "| weighted time ratio | all matched rows | %.3f | %.3f ms | %.3f ms | %d |\n",
         summary.time_ratio,
         summary.jl_ms,
         summary.py_ms,
         summary.matched
     )
+    if summary_student_t.total > 0
+        @printf(
+            io,
+            "| weighted time ratio | excluding StudentT | %.3f | %.3f ms | %.3f ms | %d |\n",
+            summary_non_student_t.time_ratio,
+            summary_non_student_t.jl_ms,
+            summary_non_student_t.py_ms,
+            summary_non_student_t.matched
+        )
+        @printf(
+            io,
+            "| weighted time ratio | StudentT only | %.3f | %.3f ms | %.3f ms | %d |\n",
+            summary_student_t.time_ratio,
+            summary_student_t.jl_ms,
+            summary_student_t.py_ms,
+            summary_student_t.matched
+        )
+    end
     if summary.peak_rows > 0
         @printf(
             io,
-            "| weighted tracemalloc ratio | %.3f | %.1f KiB | %.1f KiB | %d |\n",
+            "| weighted tracemalloc ratio | all matched rows | %.3f | %.1f KiB | %.1f KiB | %d |\n",
             summary.peak_ratio,
             summary.jl_peak_kib,
             summary.py_peak_kib,
             summary.peak_rows
         )
     else
-        println(io, "| weighted tracemalloc ratio | n/a | n/a | n/a | 0 |")
+        println(io, "| weighted tracemalloc ratio | all matched rows | n/a | n/a | n/a | 0 |")
     end
-    if summary.rss_rows > 0
+    if summary.rss_rows > 0 && summary.rss_ratio !== nothing
         @printf(
             io,
-            "| weighted RSS delta ratio | %.3f | %.1f KiB | %.1f KiB | %d |\n\n",
+            "| weighted RSS delta ratio | all matched rows | %.3f | %.1f KiB | %.1f KiB | %d |\n\n",
             summary.rss_ratio,
             summary.jl_rss_kib,
             summary.py_rss_delta_kib,
             summary.rss_rows
         )
+    elseif summary.rss_rows > 0
+        @printf(
+            io,
+            "| weighted RSS delta ratio | all matched rows | n/a | %.1f KiB | %.1f KiB | %d |\n\n",
+            summary.jl_rss_kib,
+            summary.py_rss_delta_kib,
+            summary.rss_rows
+        )
     else
-        println(io, "| weighted RSS delta ratio | n/a | n/a | n/a | 0 |\n")
+        println(io, "| weighted RSS delta ratio | all matched rows | n/a | n/a | n/a | 0 |\n")
     end
     println(io, "")
     println(
diff --git a/benchmark/results/README.md b/benchmark/results/README.md
deleted file mode 100644
index 15a6607..0000000
--- a/benchmark/results/README.md
+++ /dev/null
@@ -1,14 +0,0 @@
-# Benchmark Results
-
-This directory stores generated benchmark outputs and comparison reports.
-
-## Typical files
-
-- `<label>.json`: Julia `BenchmarkTools` results
-- `<label>_memory.json`: Julia retained-RSS sidecar
-- `<label>_solutions.json`: Julia integrate-case solution and exactness sidecar
-- `qmcpy_<label>.json`: QMCPy benchmark results
-- `compare_*.md`: Julia-vs-Julia comparison reports
-- `compare_python*.md`: Julia-vs-QMCPy comparison reports
-
-The JSON and Markdown outputs here are machine-specific benchmark artifacts and are intentionally not tracked in git.
diff --git a/benchmark/runbenchmarks.jl b/benchmark/runbenchmarks.jl
index 572e370..481972c 100644
--- a/benchmark/runbenchmarks.jl
+++ b/benchmark/runbenchmarks.jl
@@ -24,6 +24,8 @@ Pkg.resolve()
 Pkg.instantiate()
 
 using BenchmarkTools
+using Dates
+using LinearAlgebra
 using Statistics
 using Printf
 using JSON3
@@ -31,12 +33,30 @@ using JSON3
 include(joinpath(@__DIR__, "benchmarks.jl"))   # defines SUITE (and silences expected warnings)
 
 label = isempty(ARGS) ? "latest" : ARGS[1]
+generated_at = Dates.format(Dates.now(), dateformat"yyyy-mm-ddTHH:MM:SS")
 
-function current_rss_kib()
-    out = read(`ps -o rss= -p $(getpid())`, String)
-    return parse(Float64, strip(out))
+function benchmark_metadata(label::AbstractString)
+    return Dict(
+        "label" => label,
+        "generated_at" => generated_at,
+        "julia_version" => string(VERSION),
+        "blas_threads" => BLAS.get_num_threads(),
+        "benchmark_config" => Dict(
+            "leaf_samples" => BENCH_SAMPLES,
+            "integrate_samples" => INTEGRATE_BENCH_SAMPLES,
+            "student_t_samples" => STUDENT_T_BENCH_SAMPLES,
+        ),
+    )
 end
 
+current_rss_kib() =
+    try
+        out = read(`ps -o rss= -p $(getpid())`, String)
+        return parse(Float64, strip(out))
+    catch
+        return nothing
+    end
+
 function measure_rss_delta_kib(bench::BenchmarkTools.Benchmark)
     # Warm once so retained-memory deltas reflect steady-state behavior rather
     # than first-call compilation/cache effects inside the benchmarkable body.
@@ -52,6 +72,7 @@ function measure_rss_delta_kib(bench::BenchmarkTools.Benchmark)
     GC.gc()
     GC.gc()
     rss_after = current_rss_kib()
+    (rss_before === nothing || rss_after === nothing) && return nothing
     return max(rss_after - rss_before, 0.0)
 end
 
@@ -61,8 +82,10 @@ function collect_rss_deltas(group::BenchmarkTools.BenchmarkGroup)
         subgroup = group[group_name]
         rows = Dict{String, Any}()
         for bench_name in sort(collect(keys(subgroup)))
-            rows[string(bench_name)] =
-                Dict("rss_delta_kib" => measure_rss_delta_kib(subgroup[bench_name]))
+            entry = Dict{String, Any}()
+            rss_delta_kib = measure_rss_delta_kib(subgroup[bench_name])
+            rss_delta_kib !== nothing && (entry["rss_delta_kib"] = rss_delta_kib)
+            rows[string(bench_name)] = entry
         end
         out[string(group_name)] = rows
     end
@@ -209,6 +232,9 @@ end
 
 println("QMC.jl Benchmarks")
 println("="^70)
+println(
+    "Benchmark config: BLAS threads=$(BLAS.get_num_threads()), integrate samples=$(INTEGRATE_BENCH_SAMPLES), StudentT samples=$(STUDENT_T_BENCH_SAMPLES)",
+)
 
 results = run_suite(SUITE)
 julia_memory = collect_rss_deltas(SUITE)
@@ -223,14 +249,16 @@ for group_name in sort(collect(keys(results)))
     group = results[group_name]
     for bench_name in sort(collect(keys(group)))
         med = median(group[bench_name])
-        rss_delta = julia_memory[string(group_name)][string(bench_name)]["rss_delta_kib"]
+        rss_entry = julia_memory[string(group_name)][string(bench_name)]
+        rss_delta = haskey(rss_entry, "rss_delta_kib") ? rss_entry["rss_delta_kib"] : nothing
+        rss_txt = rss_delta === nothing ? "n/a" : @sprintf("%.1f KiB", rss_delta)
         @printf(
-            "  %-45s  %10.3f ms  (%d allocs, %.1f KiB, rss Δ %.1f KiB)\n",
+            "  %-45s  %10.3f ms  (%d allocs, %.1f KiB, rss Δ %s)\n",
             bench_name,
             med.time / 1e6,
             med.allocs,
             med.memory / 1024,
-            rss_delta
+            rss_txt
         )
     end
 end
@@ -241,23 +269,18 @@ resdir = joinpath(@__DIR__, "results")
 mkpath(resdir)
 outfile = joinpath(resdir, "$(label).json")
 BenchmarkTools.save(outfile, results)
+julia_meta = benchmark_metadata(label)
 memfile = joinpath(resdir, "$(label)_memory.json")
 # JSON3.pretty(str) prints to stdout and returns `nothing`; pass an IO target so the
 # pretty-printed JSON is written to the file instead.
 open(memfile, "w") do io
-    JSON3.pretty(
-        io,
-        JSON3.write(Dict("julia_version" => string(VERSION), "results" => julia_memory)),
-    )
+    JSON3.pretty(io, JSON3.write(merge(julia_meta, Dict("results" => julia_memory))))
 end
 println("\nResults saved to benchmark/results/$(label).json")
 println("Julia memory sidecar saved to benchmark/results/$(label)_memory.json")
 
 solfile = joinpath(resdir, "$(label)_solutions.json")
 open(solfile, "w") do io
-    JSON3.pretty(
-        io,
-        JSON3.write(Dict("julia_version" => string(VERSION), "solutions" => julia_solutions)),
-    )
+    JSON3.pretty(io, JSON3.write(merge(julia_meta, Dict("solutions" => julia_solutions))))
 end
 println("Julia solution sidecar saved to benchmark/results/$(label)_solutions.json")
diff --git a/demos/README.md b/demos/README.md
index 48b1272..3e9af16 100644
--- a/demos/README.md
+++ b/demos/README.md
@@ -47,15 +47,31 @@ jupyter notebook demos/quickstart.ipynb
 
 When Jupyter opens, choose the `QMC` Julia kernel for the notebook. The Python environment only launches Jupyter; it does not run the notebook code.
 
-| Notebook | Topic |
-|---|---|
-| `quickstart.ipynb` | Keister integral with all 5 stopping criteria |
-| `pricing_options.ipynb` | European, Asian, Lookback, Digital options |
-| `digital_net_b2.ipynb` | Sobol' sequences, scrambling, replications |
-| `lattice.ipynb` | Lattice orderings, shifts, Kuo vectors |
-| `lebesgue_integration.ipynb` | Lebesgue measure integration |
-| `some_true_measures.ipynb` | All 10 true measures demonstrated |
-| `gbm_demo.ipynb` | Geometric Brownian Motion for finance |
-| `qmc.jl_intro.ipynb` | Package tour: building blocks explained |
-| `sample_scatter_plots.ipynb` | Point set visualizations and statistics |
-| `elliptic_pde.ipynb` | 1D elliptic PDE uncertainty quantification with MC and lattice QMC |
+Current notebook inventory (26 notebooks):
+
+- `acceptance_rejection.ipynb`
+- `asian_option_mlqmc.ipynb`
+- `control_variates.ipynb`
+- `digital_net_b2.ipynb`
+- `elliptic_pde.ipynb`
+- `financial_option_ml.ipynb`
+- `gbm_demo.ipynb`
+- `iris.ipynb`
+- `kronecker.ipynb`
+- `lattice.ipynb`
+- `lattice_random_generator.ipynb`
+- `lebesgue_integration.ipynb`
+- `linear_scrambled_halton.ipynb`
+- `nei_demo.ipynb`
+- `plot_proj_function.ipynb`
+- `pricing_options.ipynb`
+- `qei_demo.ipynb`
+- `qmc.jl_intro.ipynb`
+- `quickstart.ipynb`
+- `ray_tracing.ipynb`
+- `sample_scatter_plots.ipynb`
+- `sensitivity_indices.ipynb`
+- `some_true_measures.ipynb`
+- `umbridge.ipynb`
+- `vectorized_qmc.ipynb`
+- `vectorized_qmc_bayes.ipynb`
diff --git a/docs/src/api/discrete_distributions.md b/docs/src/api/discrete_distributions.md
index af0535e..6aab9cc 100644
--- a/docs/src/api/discrete_distributions.md
+++ b/docs/src/api/discrete_distributions.md
@@ -22,6 +22,21 @@ DigitalNetAnyBases
 Faure
 ```
 
+`DigitalNetB2` now covers the main QMCPy-style constructor surface used in
+this repository: direct custom matrices, LDData-format files or names, order
+aliases, widened `t` precision, `msb` handling, and higher-order interlacing
+through `alpha`. The bundled constructor path still uses 1024 raw Sobol'
+dimensions, so `alpha > 1` reduces the default effective-dimension cap to
+`floor(1024 / alpha)` unless you provide a larger LDData-style matrix source
+explicitly.
+
+`Lattice` likewise covers the main QMCPy-style constructor surface used here:
+direct custom generating vectors, the integer shortcut for random odd vectors,
+order aliases, and LDData-format file/name/URL inputs. Local LDData-style
+files are validated against their declared `d_limit` / `n_limit` metadata, and
+multilevel `spawn_dd` preserves explicit custom vectors when the stored source
+is long enough for the requested dimension.
+
 ## Sampling
 
 ```@docs
diff --git a/docs/src/api/stopping_criteria.md b/docs/src/api/stopping_criteria.md
index cae1416..a2f89e1 100644
--- a/docs/src/api/stopping_criteria.md
+++ b/docs/src/api/stopping_criteria.md
@@ -8,6 +8,13 @@ AbstractStoppingCriterion
 
 ## Algorithms
 
+`PFGPCI` is included below for API visibility, but it is currently a stub:
+constructing it works and `integrate` throws a descriptive "not yet
+implemented" error until a Julia GP backend is added.
+
+`CubQMCNetGSingle` remains available as an exported backward-compatible alias
+for `CubQMCNetG`.
+
 ```@docs
 CubMCCLT
 CubMCCLTVec
diff --git a/docs/src/api/utilities.md b/docs/src/api/utilities.md
index 29cd0f5..e646e07 100644
--- a/docs/src/api/utilities.md
+++ b/docs/src/api/utilities.md
@@ -28,6 +28,9 @@ to_float
 
 ## Fast Walsh-Hadamard Transform
 
+`gpu_fwht!` is currently a CPU-fallback placeholder. Its docstring is listed
+here for API completeness, but it does not yet provide true GPU acceleration.
+
 ```@docs
 fwht!
 fwht
diff --git a/docs/src/ci-testing.md b/docs/src/ci-testing.md
index 9476e7d..a1ba203 100644
--- a/docs/src/ci-testing.md
+++ b/docs/src/ci-testing.md
@@ -116,7 +116,7 @@ Tests are organized to mirror the Python QMCSoftware test suite:
 | File | Coverage |
 |------|----------|
 | `test_discrete_distributions.jl` | IIDStdUniform, Lattice, DigitalNetB2, Halton, Kronecker, scrambling, windowed sampling |
-| `test_true_measures.jl` | All true measures including AcceptanceRejection, DistributionsWrapper |
+| `test_true_measures.jl` | All true measures including AcceptanceRejection, AcceptanceRejectionReal, DistributionsWrapper |
 | `test_integrands.jl` | All integrands including FinancialOption variants, SensitivityIndices, BayesianLRCoeffs |
 | `test_kernels.jl` | Shift-invariant, digital-shift-invariant, Matern, Gaussian, combined kernels |
 | `test_stopping_criteria.jl` | All stopping criteria, rel\_tol, resume/checkpoint, IterationLog |
diff --git a/docs/src/components.md b/docs/src/components.md
index e1db606..da9c406 100644
--- a/docs/src/components.md
+++ b/docs/src/components.md
@@ -10,6 +10,18 @@ Generates low-discrepancy point sets in ``[0,1)^d``. Available generators:
 - **`Lattice`** — rank-1 integration lattice with optional random shift
 - **`DigitalNetB2`** — Sobol' digital net (base 2) with optional scrambling
 - **`Halton`** — Halton sequence using prime bases
+- **`Kronecker`** — Kronecker / irrational-rotation low-discrepancy sequence
+- **`DigitalNetAnyBases`** / **`Faure`** — more general digital-net constructions
+
+`DigitalNetB2` also supports QMCPy-style constructor options such as
+custom generating matrices, `order` aliases, widened `t` bit depth, `msb`
+handling for custom matrices, and higher-order interlacing via `alpha > 1`.
+
+`Lattice` also supports QMCPy-style constructor inputs such as direct custom
+generating vectors, the integer shortcut for random odd vectors, order
+aliases, and LDData-format files, names, or URLs. In multilevel workflows,
+`spawn_dd` now preserves explicit lattice-vector sources when they carry enough
+entries for the requested dimension.
 
 `Lattice`, `DigitalNetB2`, and `Halton` currently rely on the QMCToolsCL
 shared library, so Python is a current runtime dependency for these QMC generators.
@@ -28,6 +40,11 @@ Transforms ``[0,1)^d`` samples to the desired probability domain:
 - **`Kumaraswamy`** — Kumaraswamy distribution
 - **`JohnsonsSU`** — Johnson's SU distribution
 - **`BernoulliCont`** — continuous Bernoulli
+- **`AcceptanceRejection`** / **`AcceptanceRejectionReal`** — deterministic acceptance-rejection samplers on the cube / real space
+- **`DistributionsWrapper`** — bridge to Distributions.jl marginals
+- **`MaternGP`** — Matern Gaussian-process prior measure
+- **`UniformTriangle`** — uniform sampling on a reference triangle
+- **`ZeroInflatedExpUniform`** — mixed zero-inflated / exponential-uniform measure
 
 ## Integrand
 
@@ -36,9 +53,13 @@ The function to be integrated:
 - **`CustomFun`** — user-supplied function
 - **`Keister`** — ``\pi^{d/2} \cos(\|x\|)`` (standard QMC test function)
 - **`Genz`** — six Genz test functions (oscillatory, product peak, etc.)
-- **`FinancialOption`** — financial option pricing (European, Asian, lookback, digital)
+- **`AsianOption`** / **`FinancialOption`** — financial option pricing (European, Asian, lookback, digital, barrier)
+- **`FinancialOptionML`** — multilevel financial option interface
 - **`BoxIntegral`**, **`Linear0`**, **`Sin1D`**, **`Ishigami`**, **`Hartmann6D`** — standard test integrands
 - **`Multimodal2D`**, **`FourBranch2D`** — reliability analysis test functions
+- **`SensitivityIndices`** — Sobol sensitivity-index estimator
+- **`BayesianLRCoeffs`** — Bayesian logistic-regression coefficient posterior
+- **`UMBridgeWrapper`** — UMBridge.jl-backed external model integrand
 
 ### Multilevel Integrands
 
@@ -60,16 +81,21 @@ Adaptive algorithms that determine sample size:
 ### Single-Level Methods
 
 - **`CubMCCLT`** — IID Monte Carlo with CLT confidence interval
+- **`CubMCCLTVec`** — vectorized IID Monte Carlo with CLT-based doubling
+- **`CubMCG`** — guaranteed IID Monte Carlo using Berry-Esseen-style bounds
 - **`CubQMCLatticeG`** — replicated randomized lattice rule
 - **`CubQMCNetG`** — single randomized digital net in natural (radical-inverse) order
 - **`CubQMCNetGRep`** — replicated randomized digital net
 - **`CubQMCBayesLatticeG`** — Bayesian QMC for lattices (kernel-based error bound)
 - **`CubQMCBayesNetG`** — Bayesian QMC for digital nets (WHT-based error bound)
+- **`CubQMCRepStudentT`** — replicated QMC with Student's *t* confidence intervals
+- **`PFGPCI`** — probability-of-failure GP criterion; currently exported as a stub that errors on `integrate`
 
 ### Multilevel Methods
 
 - **`CubMLMC`** — multilevel Monte Carlo (Giles 2008) with optimal sample allocation
 - **`CubMLMCCont`** — continuation MLMC with progressively tighter tolerances
+- **`CubMLQMC`** — multilevel quasi-Monte Carlo with replicated low-discrepancy rules
 - **`CubMLQMCCont`** — continuation multilevel quasi-Monte Carlo using replicated lattice/digital net rules
 
 The multilevel methods use the telescoping sum ``E[Q_L] = E[Q_0] + \sum_{\ell=1}^{L} E[Q_\ell - Q_{\ell-1}]``
diff --git a/docs/src/demos.md b/docs/src/demos.md
index 5367b55..ba7abf1 100644
--- a/docs/src/demos.md
+++ b/docs/src/demos.md
@@ -40,20 +40,34 @@ println("Asian call price: $(result.solution)")
 
 ## Jupyter Notebooks
 
-Interactive Jupyter notebooks are available in the `demos/` directory:
+The `demos/` directory currently contains 26 notebooks. Current inventory:
 
-| Notebook | Description |
-|----------|-------------|
-| `quickstart.ipynb` | Getting started with QMC.jl |
-| `qmc.jl_intro.ipynb` | Introduction to the QMC framework |
-| `lebesgue_integration.ipynb` | Lebesgue integration examples |
-| `lattice.ipynb` | Lattice rule demonstrations |
-| `digital_net_b2.ipynb` | Digital net (Sobol') demonstrations |
-| `some_true_measures.ipynb` | True measure transforms |
-| `pricing_options.ipynb` | Financial option pricing (European, Asian, lookback) |
-| `gbm_demo.ipynb` | Geometric Brownian motion and volatility analysis |
-| `sample_scatter_plots.ipynb` | Scatter plots of QMC point sets |
-| `elliptic_pde.ipynb` | Elliptic PDE with multilevel MC/QMC methods |
+- `acceptance_rejection.ipynb`
+- `asian_option_mlqmc.ipynb`
+- `control_variates.ipynb`
+- `digital_net_b2.ipynb`
+- `elliptic_pde.ipynb`
+- `financial_option_ml.ipynb`
+- `gbm_demo.ipynb`
+- `iris.ipynb`
+- `kronecker.ipynb`
+- `lattice.ipynb`
+- `lattice_random_generator.ipynb`
+- `lebesgue_integration.ipynb`
+- `linear_scrambled_halton.ipynb`
+- `nei_demo.ipynb`
+- `plot_proj_function.ipynb`
+- `pricing_options.ipynb`
+- `qei_demo.ipynb`
+- `qmc.jl_intro.ipynb`
+- `quickstart.ipynb`
+- `ray_tracing.ipynb`
+- `sample_scatter_plots.ipynb`
+- `sensitivity_indices.ipynb`
+- `some_true_measures.ipynb`
+- `umbridge.ipynb`
+- `vectorized_qmc.ipynb`
+- `vectorized_qmc_bayes.ipynb`
 
 To run a notebook:
 ```bash
diff --git a/docs/src/index.md b/docs/src/index.md
index 50d904f..9792bff 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -1,11 +1,17 @@
 # QMC.jl: Quasi-Monte Carlo Software in Julia
 
-QMC.jl is a Julia translation of the original QMC package for Quasi-Monte Carlo (QMC) numerical integration.
+QMC.jl is a Julia port of QMCPy for Quasi-Monte Carlo (QMC) numerical
+integration.
 
 `IIDStdUniform` is Julia-native. The main low-discrepancy generators `Lattice`,
 `DigitalNetB2`, and `Halton` currently rely on the QMCToolsCL shared library,
 so Python is a current runtime dependency for full QMC functionality.
 
+Most core package families are implemented. A few advanced items remain
+partial: `PFGPCI` is currently a documented/exported stub that errors on
+`integrate`, and `gpu_fwht!` is presently a CPU-fallback placeholder rather
+than a real GPU backend.
+
 ## Overview
 
 QMC methods approximate multivariate integrals using four main components:
diff --git a/src/discrete_distribution/README.md b/src/discrete_distribution/README.md
index 66025d0..25dc485 100644
--- a/src/discrete_distribution/README.md
+++ b/src/discrete_distribution/README.md
@@ -5,8 +5,14 @@ This folder implements the low-discrepancy and IID point generators that produce
 ## Main generators
 
 - `iid_std_uniform.jl`: IID pseudo-random uniform samples.
-- `lattice.jl`: rank-1 lattice rules.
-- `digital_net_b2.jl`: Sobol'/digital net base-2 construction.
+- `lattice.jl`: rank-1 lattice rules, including custom generating vectors,
+  QMCPy-style order aliases, LDData-style file/name inputs, and preserved
+  custom-vector respawning for multilevel workflows.
+- `digital_net_b2.jl`: Sobol'/digital net base-2 construction, including
+  custom generating matrices, QMCPy-style constructor aliases, and
+  higher-order interlacing through `alpha`. The built-in constructor path uses
+  1024 raw Sobol' dimensions, and larger higher-order jobs can opt into
+  explicit LDData-style sources when `alpha > 1` would otherwise hit that cap.
 - `halton.jl`: Halton sequences.
 - `kronecker.jl`: Kronecker sequences.
 - `digital_net_any_bases.jl`: more general digital-net support.
diff --git a/src/discrete_distribution/digital_net_b2.jl b/src/discrete_distribution/digital_net_b2.jl
index a341435..2e00beb 100644
--- a/src/discrete_distribution/digital_net_b2.jl
+++ b/src/discrete_distribution/digital_net_b2.jl
@@ -1,15 +1,19 @@
 """
-    DigitalNetB2(dimension::Int; randomize="LMS_DS", seed=nothing, graycode=true,
-                 replications=nothing, generating_matrices=nothing)
+    DigitalNetB2(dimension::Int; randomize="LMS_DS", seed=nothing, graycode=nothing,
+                 order=nothing, replications=nothing, generating_matrices=nothing,
+                 t=nothing, alpha=1, msb=nothing)
 
 Digital net in base 2 (Sobol' sequence) with optional scrambling.
 
 Generates Sobol' points using the published Joe-Kuo direction numbers,
-supporting up to 1024 dimensions and 2^32 points.
+supporting up to 1024 raw Sobol' dimensions and 2^32 points with the bundled
+Joe-Kuo table. When `alpha > 1`, those raw dimensions are interlaced in groups
+of `alpha`, so the bundled effective-dimension limit becomes `floor(1024/alpha)`
+unless custom generating matrices are supplied explicitly.
 
 This generator currently relies on the QMCToolsCL shared library. Install
 `qmctoolscl` into a Python visible to Julia, or set `ENV["QMC_PYTHON"]`
-before `using QMC`.
+before the first `Lattice`, `DigitalNetB2`, or `Halton` use.
 
 # Randomization options
 - `"LMS_DS"`: linear matrix scramble followed by a digital shift (default).
@@ -19,17 +23,38 @@ before `using QMC`.
 - `"none"`: no randomization (deterministic Sobol' points).
 
 # Arguments
-- `dimension::Int`: number of dimensions (up to 1024).
+- `dimension::Int`: number of output dimensions. With the bundled Joe-Kuo
+  table, the raw Sobol' dimension cap is 1024, so higher-order
+  `alpha > 1` nets have a correspondingly lower effective-dimension limit
+  unless custom matrices are supplied.
 - `randomize::String`: randomization method.
 - `seed`: optional RNG seed for reproducibility.
-- `graycode::Bool=true`: use Gray-code ordering for efficiency.
+- `graycode::Union{Nothing,Bool}=nothing`: use Gray-code ordering for efficiency.
+- `order`: QMCPy-style ordering alias. Accepts `"GRAY"` / `"GRAY CODE"` or
+  `"RADICAL INVERSE"` / `"NATURAL"`. If omitted, Julia keeps the historical
+  default `graycode=true`.
 - `replications::Union{Nothing,Int}=nothing`: number of independent randomizations.
   If set, `gen_samples` returns an `R × n × d` array.
-- `generating_matrices`: optional custom direction-number matrix. May be
-  `nothing` (default Joe-Kuo table) or an integer matrix with at least
-  `dimension` rows and at most 32 columns. When the custom matrix has `m`
-  columns, `gen_samples` supports at most `2^m` points including any `n_start`
-  offset. Custom entries must already be represented in that `m`-bit precision.
+- `generating_matrices`: optional custom direction-number source. May be
+  `nothing` (default Joe-Kuo table), a direct integer matrix with at least
+  `dimension * alpha` rows and at most 32 columns, or a QMCPy/LDData-style
+  `dnet` text file/name/URL. For direct integer matrices, when the custom
+  matrix has `m` columns, `gen_samples` supports at most `2^m` points
+  including any `n_start` offset, and each entry must already fit within that
+  `m`-bit width. LDData-style text files may carry wider source precision
+  (up to 64 bits) and their own sample-limit metadata. A string may point to a
+  local file, a bare filename (for example `"joe_kuo.6.21201.txt"`), or a
+  GitHub/LDData URL. Non-bundled LDData files download on demand.
+- `t`: optional number of bits used after randomization / float conversion.
+  Must satisfy `t_source <= t <= 64`, where `t_source` is the bit width of the
+  generating matrices. If omitted, Julia preserves the historical `t = t_source`
+  for standard nets and uses a higher-order default when `alpha > 1`.
+- `alpha`: interlacing factor for higher-order digital nets. When `alpha > 1`,
+  Julia interlaces `alpha * dimension` underlying Sobol' dimensions into the
+  requested output dimension.
+- `msb`: optional flag for custom integer matrices. `true` means entries are
+  already in the expected most-significant-bit ordering; `false` bit-reverses
+  each entry across the generating-matrix bit width before use.
 
 # Examples
 ```julia
@@ -39,6 +64,10 @@ x = gen_samples(dn, 1024)  # 1024×3 scrambled Sobol' points
 dn_r = DigitalNetB2(3; seed=7, replications=16)
 x = gen_samples(dn_r, 1024)  # 16×1024×3 array
 
+# Higher-order interlacing
+dn_alpha = DigitalNetB2(3; seed=7, alpha=2)
+x = gen_samples(dn_alpha, 256)  # 256×3 interlaced Sobol' points
+
 # High-dimensional problems
 dn_hd = DigitalNetB2(52; seed=42)
 x = gen_samples(dn_hd, 4096)  # 4096×52 for e.g. Asian option with 52 timesteps
@@ -46,6 +75,9 @@ x = gen_samples(dn_hd, 4096)  # 4096×52 for e.g. Asian option with 52 timesteps
 """
 
 const _SOBOL_BITS = 32
+const _DEFAULT_DNET_SOURCE_NAME = "joe_kuo.6.21201.txt"
+const _DEFAULT_DNET_1024_SOURCE_NAME = "joe_kuo.6.1024.txt"
+const _DNET_LDDATA_RAW_BASE = "https://raw.githubusercontent.com/QMCSoftware/LDData/main/dnet/"
 
 # ──────────────────────────────────────────────────────────────────────────────
 # Load precomputed direction numbers from data file
@@ -54,7 +86,7 @@ const _SOBOL_BITS = 32
 # ──────────────────────────────────────────────────────────────────────────────
 
 """
-    _load_direction_numbers(ndim::Int) -> Matrix{UInt32}
+    _load_direction_numbers(ndim::Int) -> Matrix{UInt64}
 
 Load the precomputed Joe-Kuo direction numbers for dimensions 1..ndim.
 Returns a ndim × 32 matrix V where V[j,k] is the k-th direction number
@@ -62,16 +94,152 @@ for dimension j, already left-shifted so that the point value is V/2^32.
 """
 function _load_direction_numbers(ndim::Int)
     B = _SOBOL_BITS
-    V = Matrix{UInt32}(undef, ndim, B)
+    V = Matrix{UInt64}(undef, ndim, B)
     @inbounds for j in 1:ndim
         base = (j - 1) * B
         for k in 1:B
-            V[j, k] = _JK_DIRECTION_MATRIX[base + k]
+            V[j, k] = UInt64(_JK_DIRECTION_MATRIX[base + k])
         end
     end
     return V
 end
 
+function _canonical_digital_net_matrix_filename(name::AbstractString)
+    token = replace(strip(name), '\\' => '/')
+    token = split(token, '?'; limit=2)[1]
+    token = split(token, '#'; limit=2)[1]
+    isempty(token) && throw(ArgumentError("generating_matrices reference cannot be empty"))
+    parts = split(token, '/')
+    basename = parts[end]
+    isempty(basename) &&
+        throw(ArgumentError("generating_matrices reference \"$name\" has no filename"))
+    return endswith(lowercase(basename), ".txt") ? basename : string(basename, ".txt")
+end
+
+function _looks_like_lddata_digital_net_reference(name::AbstractString)
+    token = lowercase(replace(strip(name), '\\' => '/'))
+    isempty(token) && return false
+    if !(occursin('/', token) || occursin('\\', name))
+        return true
+    end
+    return (
+        occursin("github.com/qmcsoftware/lddata/", token) ||
+        occursin("raw.githubusercontent.com/qmcsoftware/lddata/", token) ||
+        startswith(token, "dnet/") ||
+        occursin("/dnet/", token)
+    )
+end
+
+function _parse_digital_net_matrix_entry(token::AbstractString, path::AbstractString)
+    value = tryparse(BigInt, token)
+    isnothing(value) &&
+        throw(ArgumentError("invalid generating-matrix entry in \"$path\": \"$token\""))
+    return value
+end
+
+function _read_digital_net_matrix_file(path::AbstractString, dimension::Int)
+    contents = String[]
+    for raw_line in eachline(path)
+        line = strip(split(raw_line, '#'; limit=2)[1])
+        isempty(line) && continue
+        push!(contents, line)
+    end
+    length(contents) >= 4 || throw(
+        ArgumentError(
+            "generating_matrices file \"$path\" must contain base, d_limit, n_limit, and bit metadata",
+        ),
+    )
+    base = tryparse(Int, contents[1])
+    base == 2 ||
+        throw(ArgumentError("DigitalNetB2 requires base=2 in \"$path\", got $(contents[1])"))
+    d_limit = tryparse(Int, contents[2])
+    isnothing(d_limit) && throw(
+        ArgumentError(
+            "invalid d_limit in generating_matrices file \"$path\": \"$(contents[2])\"",
+        ),
+    )
+    dimension <= d_limit || throw(
+        ArgumentError(
+            "dimension $dimension exceeds the file-supported maximum ($d_limit) in \"$path\"",
+        ),
+    )
+    n_limit = tryparse(Int, contents[3])
+    isnothing(n_limit) && throw(
+        ArgumentError(
+            "invalid n_limit in generating_matrices file \"$path\": \"$(contents[3])\"",
+        ),
+    )
+    bit_precision = tryparse(Int, contents[4])
+    isnothing(bit_precision) && throw(
+        ArgumentError(
+            "invalid bit precision in generating_matrices file \"$path\": \"$(contents[4])\"",
+        ),
+    )
+    0 < bit_precision <= 64 || throw(
+        ArgumentError(
+            "DigitalNetB2 currently supports at most 64-bit generating matrices, got $bit_precision in \"$path\"",
+        ),
+    )
+    length(contents) >= 4 + dimension || throw(
+        ArgumentError(
+            "generating_matrices file \"$path\" contains fewer than $dimension matrix rows",
+        ),
+    )
+    rows = split.(contents[5:(4 + dimension)])
+    mmax = length(rows[1])
+    mmax > 0 ||
+        throw(ArgumentError("generating_matrices file \"$path\" has an empty matrix row"))
+    V = Matrix{UInt64}(undef, dimension, mmax)
+    max_entry = (BigInt(1) << bit_precision) - 1
+    for j in 1:dimension
+        length(rows[j]) == mmax || throw(
+            ArgumentError(
+                "inconsistent generating-matrix row width in \"$path\": expected $mmax columns, got $(length(rows[j])) on row $j",
+            ),
+        )
+        for k in 1:mmax
+            value = _parse_digital_net_matrix_entry(rows[j][k], path)
+            0 < value <= max_entry || throw(
+                ArgumentError(
+                    "generating_matrices entry at ($j, $k) in \"$path\" must lie in 1:$max_entry, got $value",
+                ),
+            )
+            try
+                V[j, k] = UInt64(value)
+            catch err
+                if err isa InexactError
+                    throw(
+                        ArgumentError(
+                            "generating_matrices entry at ($j, $k) in \"$path\" cannot be represented as UInt64: $value",
+                        ),
+                    )
+                end
+                rethrow()
+            end
+        end
+    end
+    return V, n_limit, bit_precision
+end
+
+function _download_digital_net_matrix_from_lddata(filename::AbstractString, dimension::Int)
+    url = _DNET_LDDATA_RAW_BASE * filename
+    path, io = mktemp()
+    close(io)
+    try
+        _downloads_module().download(url, path)
+        return _read_digital_net_matrix_file(path, dimension)
+    catch err
+        msg = sprint(showerror, err)
+        throw(
+            ArgumentError(
+                "failed to fetch LDData digital-net generating matrices \"$filename\" from $url: $msg",
+            ),
+        )
+    finally
+        isfile(path) && rm(path; force=true)
+    end
+end
+
 # ──────────────────────────────────────────────────────────────────────────────
 # Struct
 # ──────────────────────────────────────────────────────────────────────────────
@@ -79,7 +247,8 @@ end
     DigitalNetB2 <: AbstractDiscreteDistribution
 
 Digital net in base 2 (Sobol' sequence) with optional scrambling.
-Supports up to 1024 dimensions using Joe-Kuo direction numbers.
+Supports up to 1024 raw Sobol' dimensions when using the bundled Joe-Kuo
+direction numbers.
 See the outer constructor for full documentation.
 
 QMCToolsCL is required for sample generation.
@@ -88,8 +257,11 @@ mutable struct DigitalNetB2{R <: AbstractRNG} <: AbstractDiscreteDistribution
     dimension::Int
     randomize::String
     graycode::Bool
+    t::Int
+    alpha::Int
+    source_bits::Int
     rng::R
-    direction_nums::Matrix{UInt32}   # ndim × BITS
+    direction_nums::Matrix{UInt64}   # raw ndim × mmax generating matrices
     n_limit::Union{Nothing, Int}
     mimics::String
     replications::Union{Nothing, Int}
@@ -106,9 +278,11 @@ function _normalize_digital_net_randomize(randomize)
         return "LMS"
     elseif token_lc == "ds"
         return "DS"
-    elseif token_lc == "nus"
+    elseif token_lc in ("nus", "owen")
         return "NUS"
-    elseif token_lc in ("none", "false")
+    elseif token_lc == "true"
+        return "LMS_DS"
+    elseif token_lc in ("none", "false", "no")
         return "none"
     end
     throw(
@@ -118,6 +292,22 @@ function _normalize_digital_net_randomize(randomize)
     )
 end
 
+function _normalize_digital_net_order(order)
+    token = order isa Symbol ? String(order) : order
+    token isa AbstractString ||
+        throw(ArgumentError("order must be a string or symbol, got $(typeof(order))"))
+    token_lc = replace(lowercase(strip(token)), "_" => " ", "-" => " ")
+    token_lc == "gray code" && return true
+    token_lc == "gray" && return true
+    token_lc == "natural" && return false
+    token_lc == "radical inverse" && return false
+    throw(
+        ArgumentError(
+            "order must be \"GRAY\" or \"RADICAL INVERSE\" (\"NATURAL\" is accepted as an alias), got \"$order\"",
+        ),
+    )
+end
+
 function _coerce_direction_matrix(values::AbstractMatrix{<:Integer}, dimension::Int)
     size(values, 1) >= dimension || throw(
         ArgumentError(
@@ -130,8 +320,8 @@ function _coerce_direction_matrix(values::AbstractMatrix{<:Integer}, dimension::
             "generating_matrices must have between 1 and $_SOBOL_BITS columns, got $mmax",
         ),
     )
+    V = Matrix{UInt64}(undef, dimension, mmax)
     max_entry = (BigInt(1) << mmax) - 1
-    V = Matrix{UInt32}(undef, dimension, mmax)
     for j in 1:dimension, k in 1:mmax
         value = values[j, k]
         value > 0 || throw(
@@ -145,30 +335,59 @@ function _coerce_direction_matrix(values::AbstractMatrix{<:Integer}, dimension::
             ),
         )
         try
-            V[j, k] = UInt32(value)
+            V[j, k] = UInt64(value)
         catch err
             if err isa InexactError
                 throw(
                     ArgumentError(
-                        "generating_matrices entry at ($j, $k) cannot be represented as UInt32: $value",
+                        "generating_matrices entry at ($j, $k) cannot be represented as UInt64: $value",
                     ),
                 )
             end
             rethrow()
         end
     end
-    return V, Int(1) << mmax
+    return V, Int(1) << mmax, mmax
+end
+
+_bitreverse_width(value::UInt64, width::Int) = bitreverse(value) >> (64 - width)
+
+function _coerce_direction_matrix(
+    values::AbstractMatrix{<:Integer},
+    dimension::Int,
+    msb::Union{Nothing, Bool},
+)
+    V, n_limit, t_source = _coerce_direction_matrix(values, dimension)
+    if msb === false
+        width = size(V, 2)
+        @inbounds for j in axes(V, 1), k in axes(V, 2)
+            V[j, k] = _bitreverse_width(V[j, k], width)
+        end
+    end
+    return V, n_limit, t_source
 end
 
-function _resolve_direction_numbers(dimension::Int, generating_matrices)
+function _resolve_direction_numbers(dimension::Int, generating_matrices, msb)
     if isnothing(generating_matrices)
-        return _load_direction_numbers(dimension), Int(1) << _SOBOL_BITS
+        return _load_direction_numbers(dimension), Int(1) << _SOBOL_BITS, _SOBOL_BITS
     elseif generating_matrices isa AbstractMatrix{<:Integer}
-        return _coerce_direction_matrix(generating_matrices, dimension)
+        return _coerce_direction_matrix(generating_matrices, dimension, msb)
+    elseif generating_matrices isa AbstractString
+        if isfile(generating_matrices)
+            return _read_digital_net_matrix_file(generating_matrices, dimension)
+        end
+        _looks_like_lddata_digital_net_reference(generating_matrices) ||
+            throw(ArgumentError("generating_matrices file \"$generating_matrices\" not found"))
+        filename = _canonical_digital_net_matrix_filename(generating_matrices)
+        if filename in (_DEFAULT_DNET_SOURCE_NAME, _DEFAULT_DNET_1024_SOURCE_NAME) &&
+           dimension <= _JK_DIRECTION_MATRIX_DIMS
+            return _load_direction_numbers(dimension), Int(1) << _SOBOL_BITS, _SOBOL_BITS
+        end
+        return _download_digital_net_matrix_from_lddata(filename, dimension)
     end
     throw(
         ArgumentError(
-            "generating_matrices must be nothing or an integer matrix with direction numbers",
+            "generating_matrices must be nothing, an integer matrix with direction numbers, a local LDData-format text file, or an LDData filename/URL",
         ),
     )
 end
@@ -177,16 +396,22 @@ function DigitalNetB2(
     dimension::Int;
     randomize="LMS_DS",
     seed=nothing,
-    graycode::Bool=true,
+    graycode::Union{Nothing, Bool}=nothing,
+    order=nothing,
     replications::Union{Nothing, Int}=nothing,
     generating_matrices=nothing,
+    t=nothing,
+    alpha::Int=1,
+    msb=nothing,
 )
     dimension > 0 || throw(ArgumentError("dimension must be positive, got $dimension"))
+    alpha > 0 || throw(ArgumentError("alpha must be positive, got $alpha"))
+    raw_dimension = Base.checked_mul(dimension, alpha)
     isnothing(generating_matrices) &&
-        dimension > _JK_DIRECTION_MATRIX_DIMS &&
+        raw_dimension > _JK_DIRECTION_MATRIX_DIMS &&
         throw(
             ArgumentError(
-                "dimension $dimension exceeds the maximum supported ($_JK_DIRECTION_MATRIX_DIMS)",
+                "dimension $dimension with alpha=$alpha requires $raw_dimension raw dimensions, exceeding the maximum supported ($_JK_DIRECTION_MATRIX_DIMS)",
             ),
         )
     if !isnothing(replications)
@@ -194,12 +419,33 @@ function DigitalNetB2(
     end
 
     rng = isnothing(seed) ? Random.default_rng() : MersenneTwister(seed)
-    V, n_limit = _resolve_direction_numbers(dimension, generating_matrices)
+    V, n_limit, source_bits =
+        _resolve_direction_numbers(raw_dimension, generating_matrices, msb)
+    graycode_from_order = isnothing(order) ? nothing : _normalize_digital_net_order(order)
+    if !isnothing(graycode_from_order) &&
+       !isnothing(graycode) &&
+       graycode != graycode_from_order
+        throw(ArgumentError("graycode=$(graycode) is inconsistent with order=\"$order\""))
+    end
+    graycode_bool =
+        isnothing(graycode_from_order) ? (isnothing(graycode) ? true : graycode) :
+        graycode_from_order
+    interlaced_cap = min(64, Base.checked_mul(alpha, source_bits))
+    t_bits = isnothing(t) ? (alpha == 1 ? source_bits : interlaced_cap) : t
+    t_bits isa Integer || throw(ArgumentError("t must be an integer, got $(typeof(t))"))
+    source_bits <= t_bits <= 64 ||
+        throw(ArgumentError("t must satisfy $source_bits <= t <= 64, got $t_bits"))
+    randomize_norm = _normalize_digital_net_randomize(randomize)
+    t_actual =
+        randomize_norm == "none" && alpha > 1 ? min(interlaced_cap, Int(t_bits)) : Int(t_bits)
 
     return DigitalNetB2(
         dimension,
-        _normalize_digital_net_randomize(randomize),
-        graycode,
+        randomize_norm,
+        graycode_bool,
+        t_actual,
+        alpha,
+        source_bits,
         rng,
         V,
         n_limit,
@@ -211,7 +457,8 @@ end
 # ──────────────────────────────────────────────────────────────────────────────
 # Matrix-space linear matrix scramble (LMS) over GF(2)
 #
-# Applies a random lower-triangular 32×32 matrix L (independent per dimension)
+# Applies a random lower-triangular `source_bits × source_bits` matrix `L`
+# (independent per dimension)
 # to each column of the generating matrix V, producing a scrambled copy in
 # which each direction number V[j,b] is replaced by L_j * V[j,b] over GF(2).
 #
@@ -221,31 +468,31 @@ end
 #
 # Returns a d × mmax Matrix{UInt64}.
 # ──────────────────────────────────────────────────────────────────────────────
-function _lms_direction_matrix(V::Matrix{UInt32}, d::Int, rng::AbstractRNG)
+function _lms_direction_matrix(
+    V::AbstractMatrix{<:Unsigned},
+    d::Int,
+    source_bits::Int,
+    rng::AbstractRNG,
+)
     mmax = size(V, 2)
     V_scr = Matrix{UInt64}(undef, d, mmax)
     for j in 1:d
         # Random lower-triangular matrix L over GF(2): L[k] is row k,
-        # with the diagonal bit at position (mmax-k) and random bits below it.
-        L = Vector{UInt32}(undef, mmax)
-        for k in 1:mmax
-            diag_bit = UInt32(1) << (mmax - k)
-            below_mask = diag_bit - UInt32(1)
-            L[k] = diag_bit | (rand(rng, UInt32) & below_mask)
+        # with the diagonal bit at position (source_bits-k) and random bits below it.
+        L = Vector{UInt64}(undef, source_bits)
+        for k in 1:source_bits
+            diag_bit = UInt64(1) << (source_bits - k)
+            below_mask = diag_bit - UInt64(1)
+            L[k] = diag_bit | (rand(rng, UInt64) & below_mask)
         end
         for b in 1:mmax
-            v = V[j, b]
-            newv = UInt32(0)
-            for k in 1:mmax
+            v = UInt64(V[j, b])
+            newv = UInt64(0)
+            for k in 1:source_bits
                 # Bit k of L*v = parity(L[k] & v)
                 t = L[k] & v
-                t = xor(t, t >> 16)
-                t = xor(t, t >> 8)
-                t = xor(t, t >> 4)
-                t = xor(t, t >> 2)
-                t = xor(t, t >> 1)
-                if t & UInt32(1) == UInt32(1)
-                    newv |= UInt32(1) << (mmax - k)
+                if isodd(count_ones(t))
+                    newv |= UInt64(1) << (source_bits - k)
                 end
             end
             V_scr[j, b] = UInt64(newv)
@@ -267,58 +514,178 @@ function _direction_matrix_to_C(V_mat::AbstractMatrix{T}, d::Int) where {T}
     return C_flat
 end
 
+function _rand_tbit_uint64s(rng::AbstractRNG, t::Int, n::Int)
+    t == 64 && return rand(rng, UInt64, n)
+    mask = (UInt64(1) << t) - UInt64(1)
+    return rand(rng, UInt64, n) .& mask
+end
+
+_left_shift_words(t::Int, bits::Int, R::Int) = fill(UInt64(t - bits), R)
+_interlaced_bits(alpha::Int, source_bits::Int, t::Int) = min(alpha * source_bits, t)
+
+# The fused `dnb2_gen_*_float` kernels fold the t-bit alignment left-shift
+# (`lshifts`) into the digital-shift pass, which they skip entirely when
+# `apply_shift == 0x00`. Whenever `t` exceeds the current bit width the points
+# still need that left-shift to land in [0, 1) at the correct scale, so the
+# pass must run even when there is no randomization shift. Enabling it with an
+# all-zero `shiftsb` performs a no-op digital shift while still applying the
+# alignment, matching the unfused gen → digital_shift → integer_to_float path.
+_apply_shift_flag(lshifts::AbstractVector{UInt64}) = any(!iszero, lshifts) ? 0x01 : 0x00
+
+function _interlace_direction_matrices(
+    C::Vector{UInt64},
+    R::Int,
+    d_alpha::Int,
+    mmax::Int,
+    d::Int,
+    tmax::Int,
+    tmax_alpha::Int,
+    alpha::Int,
+)
+    C_alpha = Vector{UInt64}(undef, R * d_alpha * mmax)
+    _c_dnb2_interlace!(R, d_alpha, mmax, d, tmax, tmax_alpha, alpha, C, C_alpha)
+    return C_alpha
+end
+
+_rowmaj_to_Rnxd_uint64(buf::Vector{UInt64}, R::Int, n::Int, d::Int) =
+    permutedims(reshape(buf, d, n, R), (3, 2, 1))
+_Rnxd_to_rowmaj_uint64(x::Array{UInt64, 3}) = vec(permutedims(x, (3, 2, 1)))
+
+function _interlace_point_buffer(
+    xb_buf::Vector{UInt64},
+    R::Int,
+    n::Int,
+    d_alpha::Int,
+    d::Int,
+    tmax::Int,
+    tmax_alpha::Int,
+    alpha::Int,
+)
+    xb_rnd = _rowmaj_to_Rnxd_uint64(xb_buf, R, n, d)
+    xb_rdn = permutedims(xb_rnd, (1, 3, 2))
+    C = _Rnxd_to_rowmaj_uint64(xb_rdn)
+    C_alpha = _interlace_direction_matrices(C, R, d_alpha, n, d, tmax, tmax_alpha, alpha)
+    xb_alpha_rdn = permutedims(reshape(C_alpha, n, d_alpha, R), (3, 2, 1))
+    xb_alpha_rnd = permutedims(xb_alpha_rdn, (1, 3, 2))
+    return _Rnxd_to_rowmaj_uint64(xb_alpha_rnd)
+end
+
 # ──────────────────────────────────────────────────────────────────────────────
 # Single-replication and multi-replication generation via qmctoolscl binary net
 # generation, optional digital shifts, and integer-to-float conversion.
 # ──────────────────────────────────────────────────────────────────────────────
 function _gen_single_replication(dd::DigitalNetB2, n::Int; n_start::Int=0)
     d = dd.dimension
+    d_raw = dd.alpha == 1 ? d : dd.alpha * d
     V = dd.direction_nums
     mmax = size(V, 2)
 
     if dd.randomize == "LMS_DS" || dd.randomize == "LMS"
-        V_scr = _lms_direction_matrix(V, d, dd.rng)
-        C_flat = _direction_matrix_to_C(V_scr, d)
-        lshifts = zeros(UInt64, 1)
+        curr_bits =
+            dd.alpha == 1 ? dd.source_bits : _interlaced_bits(dd.alpha, dd.source_bits, dd.t)
+        V_scr = _lms_direction_matrix(V, d_raw, dd.source_bits, dd.rng)
+        C_flat = if dd.alpha == 1
+            _direction_matrix_to_C(V_scr, d_raw)
+        else
+            _interlace_direction_matrices(
+                _direction_matrix_to_C(V_scr, d_raw),
+                1,
+                d,
+                mmax,
+                d_raw,
+                dd.source_bits,
+                curr_bits,
+                dd.alpha,
+            )
+        end
+        lshifts = _left_shift_words(dd.t, curr_bits, 1)
         if dd.randomize == "LMS_DS"
-            shiftsb = UInt64.(rand(dd.rng, UInt32, d))
+            shiftsb = _rand_tbit_uint64s(dd.rng, dd.t, d)
             apply_shift = 0x01
         else
             shiftsb = zeros(UInt64, d)
-            apply_shift = 0x00
+            apply_shift = _apply_shift_flag(lshifts)
         end
     elseif dd.randomize == "DS"
-        C_flat = _direction_matrix_to_C(V, d)
-        lshifts = zeros(UInt64, 1)
-        shiftsb = UInt64.(rand(dd.rng, UInt32, d))
+        curr_bits =
+            dd.alpha == 1 ? dd.source_bits : _interlaced_bits(dd.alpha, dd.source_bits, dd.t)
+        C_flat = if dd.alpha == 1
+            _direction_matrix_to_C(V, d_raw)
+        else
+            _interlace_direction_matrices(
+                _direction_matrix_to_C(V, d_raw),
+                1,
+                d,
+                mmax,
+                d_raw,
+                dd.source_bits,
+                curr_bits,
+                dd.alpha,
+            )
+        end
+        lshifts = _left_shift_words(dd.t, curr_bits, 1)
+        shiftsb = _rand_tbit_uint64s(dd.rng, dd.t, d)
         apply_shift = 0x01
     elseif dd.randomize == "NUS"
-        C_flat = _direction_matrix_to_C(V, d)
+        C_flat = _direction_matrix_to_C(V, d_raw)
         # Generate per-dimension seeds for Owen scrambling
-        dim_seeds = UInt64.(rand(dd.rng, UInt64, d))
+        dim_seeds = rand(dd.rng, UInt64, d_raw)
 
-        xb_buf = Vector{UInt64}(undef, n * d)
+        xb_buf = Vector{UInt64}(undef, n * d_raw)
         if dd.graycode
-            _c_dnb2_gen_gray!(1, n, d, n_start, mmax, C_flat, xb_buf)
+            _c_dnb2_gen_gray!(1, n, d_raw, n_start, mmax, C_flat, xb_buf)
         else
-            _c_dnb2_gen_natural!(1, n, d, n_start, mmax, C_flat, xb_buf)
+            _c_dnb2_gen_natural!(1, n, d_raw, n_start, mmax, C_flat, xb_buf)
         end
 
         # Apply Owen scrambling in-place
-        _owen_scramble_points!(xb_buf, n, d, mmax, dim_seeds)
+        _owen_scramble_points!(xb_buf, n, d_raw, dd.source_bits, dim_seeds)
 
-        tmaxes = fill(UInt64(mmax), 1)
+        point_bits = dd.source_bits
+        if dd.t > point_bits
+            xb_shifted = similar(xb_buf)
+            _c_dnb2_digital_shift!(
+                1,
+                n,
+                d_raw,
+                1,
+                _left_shift_words(dd.t, point_bits, 1),
+                xb_buf,
+                zeros(UInt64, d_raw),
+                xb_shifted,
+            )
+            xb_buf = xb_shifted
+            point_bits = dd.t
+        end
+        if dd.alpha > 1
+            xb_buf = _interlace_point_buffer(xb_buf, 1, n, d, d_raw, point_bits, dd.t, dd.alpha)
+        end
+        tmaxes = fill(UInt64(dd.t), 1)
         x_buf = Vector{Float64}(undef, n * d)
         _c_dnb2_integer_to_float!(1, n, d, tmaxes, xb_buf, x_buf)
         return _rowmaj_to_nxd(x_buf, n, d)
     else  # "none"
-        C_flat = _direction_matrix_to_C(V, d)
-        lshifts = zeros(UInt64, 1)
+        curr_bits = dd.alpha == 1 ? dd.source_bits : dd.t
+        C_flat = if dd.alpha == 1
+            _direction_matrix_to_C(V, d_raw)
+        else
+            _interlace_direction_matrices(
+                _direction_matrix_to_C(V, d_raw),
+                1,
+                d,
+                mmax,
+                d_raw,
+                dd.source_bits,
+                curr_bits,
+                dd.alpha,
+            )
+        end
+        lshifts = _left_shift_words(dd.t, curr_bits, 1)
         shiftsb = zeros(UInt64, d)
-        apply_shift = 0x00
+        apply_shift = _apply_shift_flag(lshifts)
     end
 
-    tmaxes = fill(UInt64(mmax), 1)
+    tmaxes = fill(UInt64(dd.t), 1)
     x_buf = Vector{Float64}(undef, n * d)
     if _HAS_DNB2_FUSED[]
         if dd.graycode
@@ -391,50 +758,97 @@ function gen_samples(dd::DigitalNetB2, n::Int; n_start::Int=0)
 
     R = dd.replications
     d = dd.dimension
+    d_raw = dd.alpha == 1 ? d : dd.alpha * d
     V = dd.direction_nums
     mmax = size(V, 2)
 
     if dd.randomize == "LMS_DS" || dd.randomize == "LMS"
-        # Build R independently scrambled generating matrices (r_x = R)
-        C_flat = Vector{UInt64}(undef, R * d * mmax)
+        curr_bits =
+            dd.alpha == 1 ? dd.source_bits : _interlaced_bits(dd.alpha, dd.source_bits, dd.t)
+        C_scr = Vector{UInt64}(undef, R * d_raw * mmax)
         for r in 1:R
-            V_scr = _lms_direction_matrix(V, d, dd.rng)
-            base = (r - 1) * d * mmax
-            for j in 1:d, b in 1:mmax
-                C_flat[base + (j - 1) * mmax + (b - 1) + 1] = UInt64(V_scr[j, b])
+            V_scr = _lms_direction_matrix(V, d_raw, dd.source_bits, dd.rng)
+            base = (r - 1) * d_raw * mmax
+            for j in 1:d_raw, b in 1:mmax
+                C_scr[base + (j - 1) * mmax + (b - 1) + 1] = UInt64(V_scr[j, b])
             end
         end
+        C_flat =
+            dd.alpha == 1 ? C_scr :
+            _interlace_direction_matrices(
+                C_scr,
+                R,
+                d,
+                mmax,
+                d_raw,
+                dd.source_bits,
+                curr_bits,
+                dd.alpha,
+            )
         r_x = R
-        lshifts = zeros(UInt64, R)
+        lshifts = _left_shift_words(dd.t, curr_bits, R)
         if dd.randomize == "LMS_DS"
-            shiftsb = UInt64.(rand(dd.rng, UInt32, R * d))
+            shiftsb = _rand_tbit_uint64s(dd.rng, dd.t, R * d)
             apply_shift = 0x01
         else
             shiftsb = zeros(UInt64, R * d)
-            apply_shift = 0x00
+            apply_shift = _apply_shift_flag(lshifts)
         end
     elseif dd.randomize == "DS"
-        C_flat = _direction_matrix_to_C(V, d)
+        curr_bits =
+            dd.alpha == 1 ? dd.source_bits : _interlaced_bits(dd.alpha, dd.source_bits, dd.t)
+        C_flat =
+            dd.alpha == 1 ? _direction_matrix_to_C(V, d_raw) :
+            _interlace_direction_matrices(
+                _direction_matrix_to_C(V, d_raw),
+                1,
+                d,
+                mmax,
+                d_raw,
+                dd.source_bits,
+                curr_bits,
+                dd.alpha,
+            )
         r_x = 1
-        lshifts = zeros(UInt64, 1)
-        shiftsb = UInt64.(rand(dd.rng, UInt32, R * d))
+        lshifts = _left_shift_words(dd.t, curr_bits, 1)
+        shiftsb = _rand_tbit_uint64s(dd.rng, dd.t, R * d)
         apply_shift = 0x01
     elseif dd.randomize == "NUS"
-        C_flat = _direction_matrix_to_C(V, d)
+        C_flat = _direction_matrix_to_C(V, d_raw)
         # Generate unscrambled binary points (single set)
-        xb_base = Vector{UInt64}(undef, n * d)
+        xb_base = Vector{UInt64}(undef, n * d_raw)
         if dd.graycode
-            _c_dnb2_gen_gray!(1, n, d, n_start, mmax, C_flat, xb_base)
+            _c_dnb2_gen_gray!(1, n, d_raw, n_start, mmax, C_flat, xb_base)
         else
-            _c_dnb2_gen_natural!(1, n, d, n_start, mmax, C_flat, xb_base)
+            _c_dnb2_gen_natural!(1, n, d_raw, n_start, mmax, C_flat, xb_base)
         end
         # Apply R independent Owen scramblings
         x_buf = Vector{Float64}(undef, R * n * d)
-        tmaxes_one = fill(UInt64(mmax), 1)
+        tmaxes_one = fill(UInt64(dd.t), 1)
         for r in 1:R
             xb_copy = copy(xb_base)
-            dim_seeds = UInt64.(rand(dd.rng, UInt64, d))
-            _owen_scramble_points!(xb_copy, n, d, mmax, dim_seeds)
+            dim_seeds = rand(dd.rng, UInt64, d_raw)
+            _owen_scramble_points!(xb_copy, n, d_raw, dd.source_bits, dim_seeds)
+            point_bits = dd.source_bits
+            if dd.t > point_bits
+                xb_shifted = similar(xb_copy)
+                _c_dnb2_digital_shift!(
+                    1,
+                    n,
+                    d_raw,
+                    1,
+                    _left_shift_words(dd.t, point_bits, 1),
+                    xb_copy,
+                    zeros(UInt64, d_raw),
+                    xb_shifted,
+                )
+                xb_copy = xb_shifted
+                point_bits = dd.t
+            end
+            if dd.alpha > 1
+                xb_copy =
+                    _interlace_point_buffer(xb_copy, 1, n, d, d_raw, point_bits, dd.t, dd.alpha)
+            end
             x_r = Vector{Float64}(undef, n * d)
             _c_dnb2_integer_to_float!(1, n, d, tmaxes_one, xb_copy, x_r)
             base = (r - 1) * n * d
@@ -442,14 +856,26 @@ function gen_samples(dd::DigitalNetB2, n::Int; n_start::Int=0)
         end
         return _rowmaj_to_Rnxd(x_buf, R, n, d)
     else  # "none"
-        C_flat = _direction_matrix_to_C(V, d)
+        curr_bits = dd.alpha == 1 ? dd.source_bits : dd.t
+        C_flat =
+            dd.alpha == 1 ? _direction_matrix_to_C(V, d_raw) :
+            _interlace_direction_matrices(
+                _direction_matrix_to_C(V, d_raw),
+                1,
+                d,
+                mmax,
+                d_raw,
+                dd.source_bits,
+                curr_bits,
+                dd.alpha,
+            )
         r_x = 1
-        lshifts = zeros(UInt64, 1)
+        lshifts = _left_shift_words(dd.t, curr_bits, 1)
         shiftsb = zeros(UInt64, R * d)
-        apply_shift = 0x00
+        apply_shift = _apply_shift_flag(lshifts)
     end
 
-    tmaxes = fill(UInt64(mmax), R)
+    tmaxes = fill(UInt64(dd.t), R)
     x_buf = Vector{Float64}(undef, R * n * d)
     if _HAS_DNB2_FUSED[]
         if dd.graycode
@@ -499,9 +925,10 @@ end
 
 function Base.show(io::IO, dd::DigitalNetB2)
     rep_str = isnothing(dd.replications) ? "" : ", R=$(dd.replications)"
+    alpha_str = dd.alpha == 1 ? "" : ", alpha=$(dd.alpha)"
     print(
         io,
-        "DigitalNetB2(d=$(dd.dimension), randomize=\"$(dd.randomize)\", graycode=$(dd.graycode)$rep_str)",
+        "DigitalNetB2(d=$(dd.dimension), randomize=\"$(dd.randomize)\", graycode=$(dd.graycode)$alpha_str$rep_str)",
     )
 end
 
diff --git a/src/discrete_distribution/lattice.jl b/src/discrete_distribution/lattice.jl
index 14b92d4..114a219 100644
--- a/src/discrete_distribution/lattice.jl
+++ b/src/discrete_distribution/lattice.jl
@@ -13,7 +13,7 @@ When `replications` is set to an integer R > 1, `gen_samples` returns an
 
 This generator currently relies on the QMCToolsCL shared library. Install
 `qmctoolscl` into a Python visible to Julia, or set `ENV["QMC_PYTHON"]`
-before `using QMC`.
+before the first `Lattice`, `DigitalNetB2`, or `Halton` use.
 
 # Arguments
 - `dimension`: number of dimensions (up to 9125).
@@ -30,8 +30,8 @@ before `using QMC`.
   - an integer `M`: generate a random odd vector `(1, v₂, …, v_d)` with
     `v_j ∈ {3, 5, …, 2^M - 1}`,
   - a local text-file path: load one integer per non-comment line, or a
-    QMCPy/LDData-style file whose first two integers are metadata followed by
-    the vector entries,
+    QMCPy/LDData-style file whose first two integers are `d_limit` / `n_limit`
+    metadata followed by exactly `d_limit` vector entries,
   - a QMCPy/LDData filename or GitHub URL (for example
     `"kuo.lattice-33002-1024-1048576.9125"` or the corresponding
     `QMCSoftware/LDData` `lattice/` URL). Non-bundled LDData files are
@@ -52,6 +52,7 @@ mutable struct Lattice{R <: AbstractRNG} <: AbstractDiscreteDistribution
     order::String
     replications::Union{Nothing, Int}
     gen_vector::Vector{UInt64}
+    source_gen_vector::Union{Nothing, Vector{UInt64}}
     shift::Matrix{Float64}   # R × d
     rng::R
     n_limit::Union{Nothing, Int}
@@ -60,6 +61,10 @@ end
 
 const _DEFAULT_LATTICE_VECTOR_NAME = "kuo.lattice-33002-1024-1048576.9125.txt"
 const _LATTICE_LDDATA_RAW_BASE = "https://raw.githubusercontent.com/QMCSoftware/LDData/main/lattice/"
+const _DOWNLOADS_PKGID =
+    Base.PkgId(Base.UUID("f43a241f-c20a-4ad4-852c-f6b1247861c6"), "Downloads")
+
+_downloads_module() = Base.require(_DOWNLOADS_PKGID)
 
 function _normalize_lattice_order(order::String)
     # Normalize order aliases to canonical tokens, matching QMCPy semantics:
@@ -108,9 +113,21 @@ end
 
 function _read_lattice_vector_file(path::AbstractString)
     values = UInt64[]
+    saw_d_limit = false
+    saw_n_limit = false
     open(path, "r") do io
         for raw_line in eachline(io)
-            line = strip(first(split(raw_line, '#'; limit=2)))
+            split_line = split(raw_line, '#'; limit=2)
+            line = strip(first(split_line))
+            if length(split_line) == 2
+                comment = replace(lowercase(strip(split_line[2])), " " => "_")
+                saw_d_limit |= occursin("d_limit", comment)
+                saw_n_limit |= occursin("n_limit", comment)
+            elseif startswith(strip(raw_line), '#')
+                comment = replace(lowercase(strip(raw_line[2:end])), " " => "_")
+                saw_d_limit |= occursin("d_limit", comment)
+                saw_n_limit |= occursin("n_limit", comment)
+            end
             isempty(line) && continue
             value = tryparse(UInt64, line)
             isnothing(value) && throw(
@@ -120,7 +137,37 @@ function _read_lattice_vector_file(path::AbstractString)
         end
     end
     isempty(values) && throw(ArgumentError("generating_vector file \"$path\" is empty"))
-    if length(values) >= 3 && Int(values[1]) == length(values) - 2
+    if saw_d_limit || saw_n_limit
+        saw_d_limit && saw_n_limit || throw(
+            ArgumentError(
+                "generating_vector file \"$path\" must specify both d_limit and n_limit metadata",
+            ),
+        )
+        length(values) >= 3 || throw(
+            ArgumentError(
+                "generating_vector file \"$path\" metadata must be followed by vector entries",
+            ),
+        )
+        values[1] <= typemax(Int) || throw(
+            ArgumentError("generating_vector file \"$path\" has d_limit too large for Int"),
+        )
+        values[2] <= typemax(Int) || throw(
+            ArgumentError("generating_vector file \"$path\" has n_limit too large for Int"),
+        )
+        d_limit = Int(values[1])
+        n_limit = Int(values[2])
+        d_limit > 0 ||
+            throw(ArgumentError("generating_vector file \"$path\" must have positive d_limit"))
+        n_limit > 0 ||
+            throw(ArgumentError("generating_vector file \"$path\" must have positive n_limit"))
+        length(values) == d_limit + 2 || throw(
+            ArgumentError(
+                "generating_vector file \"$path\" declares d_limit=$d_limit but contains $(length(values) - 2) vector entries",
+            ),
+        )
+        return values[3:end], n_limit
+    end
+    if length(values) >= 3 && values[1] <= typemax(Int) && Int(values[1]) == length(values) - 2
         values[2] <= typemax(Int) || throw(
             ArgumentError("generating_vector file \"$path\" has n_limit too large for Int"),
         )
@@ -160,7 +207,7 @@ function _download_lattice_vector_from_lddata(filename::AbstractString)
     path, io = mktemp()
     close(io)
     try
-        Downloads.download(url, path)
+        _downloads_module().download(url, path)
         return _read_lattice_vector_file(path)
     catch err
         msg = sprint(showerror, err)
@@ -195,15 +242,28 @@ function _resolve_lattice_generating_vector(dimension::Int, generating_vector, r
                 "dimension $dimension exceeds maximum supported ($_KUO_LATTICE_MAX_DIM)",
             ),
         )
-        return _KUO_LATTICE_GEN_VECTOR[1:dimension], 1 << 20
+        return _KUO_LATTICE_GEN_VECTOR[1:dimension], 1 << 20, nothing
     elseif generating_vector isa AbstractVector{<:Integer}
-        return _coerce_lattice_vector(generating_vector, dimension), nothing
+        length(generating_vector) >= dimension || throw(
+            ArgumentError(
+                "generating_vector must have at least $dimension entries, got $(length(generating_vector))",
+            ),
+        )
+        source = _coerce_lattice_vector(generating_vector, length(generating_vector))
+        return source[1:dimension], nothing, source
     elseif generating_vector isa Integer
-        return _random_lattice_vector(dimension, generating_vector, rng)
+        gv, n_limit = _random_lattice_vector(dimension, generating_vector, rng)
+        return gv, n_limit, copy(gv)
     elseif generating_vector isa AbstractString
         if isfile(generating_vector)
             values, n_limit = _read_lattice_vector_file(generating_vector)
-            return _coerce_lattice_vector(values, dimension), n_limit
+            source = _coerce_lattice_vector(values, length(values))
+            length(source) >= dimension || throw(
+                ArgumentError(
+                    "generating_vector must have at least $dimension entries, got $(length(source))",
+                ),
+            )
+            return source[1:dimension], n_limit, source
         end
         _looks_like_lddata_lattice_reference(generating_vector) ||
             throw(ArgumentError("generating_vector file \"$generating_vector\" not found"))
@@ -214,10 +274,16 @@ function _resolve_lattice_generating_vector(dimension::Int, generating_vector, r
                     "dimension $dimension exceeds maximum supported ($_KUO_LATTICE_MAX_DIM)",
                 ),
             )
-            return _KUO_LATTICE_GEN_VECTOR[1:dimension], 1 << 20
+            return _KUO_LATTICE_GEN_VECTOR[1:dimension], 1 << 20, nothing
         end
         values, n_limit = _download_lattice_vector_from_lddata(filename)
-        return _coerce_lattice_vector(values, dimension), n_limit
+        source = _coerce_lattice_vector(values, length(values))
+        length(source) >= dimension || throw(
+            ArgumentError(
+                "generating_vector must have at least $dimension entries, got $(length(source))",
+            ),
+        )
+        return source[1:dimension], n_limit, source
     else
         throw(
             ArgumentError(
@@ -242,7 +308,8 @@ function Lattice(
     R >= 1 || throw(ArgumentError("replications must be >= 1"))
 
     rng = isnothing(seed) ? Random.default_rng() : MersenneTwister(seed)
-    gv, n_limit = _resolve_lattice_generating_vector(dimension, generating_vector, rng)
+    gv, n_limit, source_gv =
+        _resolve_lattice_generating_vector(dimension, generating_vector, rng)
     shift = randomize ? rand(rng, R, dimension) : zeros(R, dimension)
 
     return Lattice(
@@ -251,6 +318,7 @@ function Lattice(
         order_lc,
         replications,
         gv,
+        source_gv,
         shift,
         rng,
         n_limit,
diff --git a/src/discrete_distribution/qmctoolscl_c.jl b/src/discrete_distribution/qmctoolscl_c.jl
index 49075eb..cf0325d 100644
--- a/src/discrete_distribution/qmctoolscl_c.jl
+++ b/src/discrete_distribution/qmctoolscl_c.jl
@@ -388,6 +388,49 @@ function _c_dnb2_digital_shift!(
     )
 end
 
+function _c_dnb2_interlace!(
+    R::Int,
+    d_alpha::Int,
+    mmax::Int,
+    d::Int,
+    tmax::Int,
+    tmax_alpha::Int,
+    alpha::Int,
+    C::Vector{UInt64},
+    C_alpha::Vector{UInt64},
+)
+    ccall(
+        (:dnb2_interlace, _qmctoolscl_lib_path()),
+        Cvoid,
+        (
+            UInt64,
+            UInt64,
+            UInt64,
+            UInt64,
+            UInt64,
+            UInt64,
+            UInt64,
+            UInt64,
+            UInt64,
+            UInt64,
+            Ptr{UInt64},
+            Ptr{UInt64},
+        ),
+        UInt64(R),
+        UInt64(d_alpha),
+        UInt64(mmax),
+        UInt64(R),
+        UInt64(d_alpha),
+        UInt64(mmax),
+        UInt64(d),
+        UInt64(tmax),
+        UInt64(tmax_alpha),
+        UInt64(alpha),
+        C,
+        C_alpha,
+    )
+end
+
 function _c_dnb2_integer_to_float!(
     R::Int,
     n::Int,
diff --git a/src/integrand/ml_integrand.jl b/src/integrand/ml_integrand.jl
index 6c81657..bdadf7d 100644
--- a/src/integrand/ml_integrand.jl
+++ b/src/integrand/ml_integrand.jl
@@ -63,11 +63,35 @@ spawn_dd(dd::IIDStdUniform, dimension::Int) = IIDStdUniform(dimension)
     spawn_dd(dd::Lattice, dimension::Int)
 
 Create a new Lattice sampler with the given dimension, preserving randomize
-and replications settings but with a fresh seed.
+and replications settings, as well as the original order and any explicit
+custom generating-vector source when enough entries are available.
 """
 function spawn_dd(dd::Lattice, dimension::Int)
     R = isnothing(dd.replications) ? nothing : dd.replications
-    return Lattice(dimension; randomize=dd.randomize, replications=R)
+    if !isnothing(dd.source_gen_vector)
+        length(dd.source_gen_vector) >= dimension || throw(
+            ArgumentError(
+                "cannot spawn Lattice dimension $dimension from a custom generating vector with only $(length(dd.source_gen_vector)) stored entries",
+            ),
+        )
+        rng = Random.default_rng()
+        reps = isnothing(R) ? 1 : R
+        shift = dd.randomize ? rand(rng, reps, dimension) : zeros(reps, dimension)
+        source = copy(dd.source_gen_vector)
+        return Lattice(
+            dimension,
+            dd.randomize,
+            dd.order,
+            R,
+            source[1:dimension],
+            source,
+            shift,
+            rng,
+            dd.n_limit,
+            dd.mimics,
+        )
+    end
+    return Lattice(dimension; randomize=dd.randomize, order=dd.order, replications=R)
 end
 
 """
@@ -78,16 +102,30 @@ randomize and replications settings but with a fresh seed.
 """
 function spawn_dd(dd::DigitalNetB2, dimension::Int)
     R = isnothing(dd.replications) ? nothing : dd.replications
-    if dimension <= size(dd.direction_nums, 1)
+    raw_dimension = Base.checked_mul(dimension, dd.alpha)
+    if raw_dimension <= size(dd.direction_nums, 1)
         return DigitalNetB2(
-            dimension;
-            randomize=dd.randomize,
-            graycode=dd.graycode,
-            replications=R,
-            generating_matrices=dd.direction_nums,
+            dimension,
+            dd.randomize,
+            dd.graycode,
+            dd.t,
+            dd.alpha,
+            dd.source_bits,
+            Random.default_rng(),
+            dd.direction_nums[1:raw_dimension, :],
+            dd.n_limit,
+            dd.mimics,
+            R,
         )
     end
-    return DigitalNetB2(dimension; randomize=dd.randomize, graycode=dd.graycode, replications=R)
+    return DigitalNetB2(
+        dimension;
+        randomize=dd.randomize,
+        graycode=dd.graycode,
+        replications=R,
+        t=dd.t,
+        alpha=dd.alpha,
+    )
 end
 
 """
diff --git a/src/stopping_criterion/cub_mc_clt.jl b/src/stopping_criterion/cub_mc_clt.jl
index c4bea57..fe2ad7e 100644
--- a/src/stopping_criterion/cub_mc_clt.jl
+++ b/src/stopping_criterion/cub_mc_clt.jl
@@ -1,8 +1,12 @@
 """
-    CubMCCLT(integrand; abs_tol=0.01, rel_tol=0.0, n_init=1024, n_max=2^30, alpha=0.01, inflate=1.2)
+    CubMCCLT(integrand; abs_tol=0.01, rel_tol=0.0, n_init=1024,
+             n_max=2^30, alpha=0.01, inflate=1.2, trace_iterations=false)
 
 IID Monte Carlo cubature with CLT-based confidence interval (two-stage method).
 
+Set `trace_iterations=true` to record an `IterationLog` in
+`result.data[:iteration_log]`.
+
 **Algorithm:**
 1. **Pilot stage:** Draw `n_init` IID samples, compute variance estimate σ̂.
 2. **Main stage:** Compute the number of *new* samples needed:
@@ -29,6 +33,7 @@ mutable struct CubMCCLT{I <: AbstractIntegrand} <: AbstractStoppingCriterion
     n_max::Int
     alpha::Float64
     inflate::Float64
+    trace_iterations::Bool
 end
 
 function CubMCCLT(
@@ -39,15 +44,27 @@ function CubMCCLT(
     n_max::Int=2^30,
     alpha::Float64=0.01,
     inflate::Float64=1.2,
+    trace_iterations::Bool=false,
 )
     n_max > 2 * n_init || throw(ArgumentError("n_max must be > 2 * n_init"))
     inflate >= 1.0 || throw(ArgumentError("inflate must be ≥ 1.0"))
     0.0 < alpha < 1.0 || throw(ArgumentError("alpha must be in (0, 1)"))
-    return CubMCCLT(integrand, abs_tol, rel_tol, n_init, n_max, alpha, inflate)
+    return CubMCCLT(
+        integrand,
+        abs_tol,
+        rel_tol,
+        n_init,
+        n_max,
+        alpha,
+        inflate,
+        trace_iterations,
+    )
 end
 
 function integrate(sc::CubMCCLT; resume::Union{Nothing, Dict{Symbol, Any}}=nothing)
+    t_start = time()
     z_star = quantile(Normal(), 1.0 - sc.alpha / 2.0)
+    log = IterationLog()
 
     # ── Stage 1: Pilot sample to estimate variance ──
     y0 = sample_and_evaluate(sc.integrand, sc.n_init)
@@ -56,6 +73,17 @@ function integrate(sc::CubMCCLT; resume::Union{Nothing, Dict{Symbol, Any}}=nothi
 
     # Determine tolerance using pilot estimate
     tol = max(sc.abs_tol, sc.rel_tol * abs(mu_hat0))
+    pilot_err = z_star * sc.inflate * sig_hat0 / sqrt(sc.n_init)
+    if sc.trace_iterations
+        push!(
+            log;
+            n=sc.n_init,
+            solution=mu_hat0,
+            error_bound=pilot_err,
+            tol=tol,
+            elapsed=time() - t_start,
+        )
+    end
 
     # Compute required main-stage sample size
     n_mu = ceil(Int, (z_star * sc.inflate * sig_hat0 / tol)^2)
@@ -84,6 +112,16 @@ function integrate(sc::CubMCCLT; resume::Union{Nothing, Dict{Symbol, Any}}=nothi
         @warn "CubMCCLT: did not converge within n_max=$(sc.n_max). " *
               "Error bound: $err, tolerance: $(max(sc.abs_tol, sc.rel_tol * abs(mu_hat)))"
     end
+    if sc.trace_iterations
+        push!(
+            log;
+            n=n_total,
+            solution=mu_hat,
+            error_bound=err,
+            tol=max(sc.abs_tol, sc.rel_tol * abs(mu_hat)),
+            elapsed=time() - t_start,
+        )
+    end
 
     data = Dict{Symbol, Any}(
         :n => n_total,
@@ -98,6 +136,9 @@ function integrate(sc::CubMCCLT; resume::Union{Nothing, Dict{Symbol, Any}}=nothi
         :converged => converged,
         :confidence_level => 1.0 - sc.alpha,
     )
+    if sc.trace_iterations
+        data[:iteration_log] = log
+    end
     return QMCResult(mu_hat, data)
 end
 
diff --git a/src/stopping_criterion/cub_mc_g.jl b/src/stopping_criterion/cub_mc_g.jl
index e5a7a09..040e758 100644
--- a/src/stopping_criterion/cub_mc_g.jl
+++ b/src/stopping_criterion/cub_mc_g.jl
@@ -1,6 +1,6 @@
 """
     CubMCG(integrand; abs_tol=0.01, rel_tol=0.0, n_init=1024,
-           n_max=2^30, alpha=0.01, inflate=1.2)
+           n_max=2^30, alpha=0.01, inflate=1.2, trace_iterations=false)
 
 Guaranteed IID Monte Carlo stopping criterion using Berry-Esseen inequalities.
 
@@ -18,6 +18,7 @@ functions with bounded kurtosis.
 - `n_max`: Maximum total number of samples.
 - `alpha`: Uncertainty level in (0, 1).
 - `inflate`: Inflation factor ≥ 1 for conservative variance estimation.
+- `trace_iterations`: record an `IterationLog` in `result.data[:iteration_log]`.
 
 # References
 1. Hickernell, Jiang, Liu, Owen. "Guaranteed conservative fixed width
@@ -33,6 +34,7 @@ mutable struct CubMCG{I <: AbstractIntegrand} <: AbstractStoppingCriterion
     inflate::Float64
     alpha_sigma::Float64
     kurtmax::Float64
+    trace_iterations::Bool
 end
 
 function CubMCG(
@@ -43,6 +45,7 @@ function CubMCG(
     n_max::Int=2^30,
     alpha::Float64=0.01,
     inflate::Float64=1.2,
+    trace_iterations::Bool=false,
 )
     abs_tol > 0 || throw(ArgumentError("abs_tol must be > 0"))
     rel_tol >= 0 || throw(ArgumentError("rel_tol must be ≥ 0"))
@@ -65,6 +68,7 @@ function CubMCG(
         inflate,
         alpha_sigma,
         kurtmax,
+        trace_iterations,
     )
 end
 
@@ -179,6 +183,7 @@ function integrate(sc::CubMCG; resume::Union{Nothing, Dict{Symbol, Any}}=nothing
     f = sc.integrand
     tm = f.true_measure
     dd = tm.dd
+    log = IterationLog()
 
     # Stage 1: pilot samples
     x0 = transform(tm, gen_samples(dd, sc.n_init))
@@ -186,10 +191,22 @@ function integrate(sc::CubMCG; resume::Union{Nothing, Dict{Symbol, Any}}=nothing
     mu0 = mean(y0)
     sig0 = std(y0; corrected=true)
     sigma_up = sc.inflate * sig0
+    tol0 = max(sc.abs_tol, sc.rel_tol * abs(mu0))
 
     if sc.rel_tol == 0.0
         # Fixed absolute tolerance
         alpha_mu = 1 - (1 - sc.alpha) / (1 - sc.alpha_sigma)
+        pilot_bound = sigma_up * _ncbinv(sc.n_init, alpha_mu, sc.kurtmax)
+        if sc.trace_iterations
+            push!(
+                log;
+                n=sc.n_init,
+                solution=mu0,
+                error_bound=pilot_bound,
+                tol=tol0,
+                elapsed=time() - t_start,
+            )
+        end
         toloversig = sc.abs_tol / max(sigma_up, 1e-300)
         n_mu, bound_hw = _nchebe(toloversig, alpha_mu, sc.kurtmax, sc.n_max, sigma_up)
 
@@ -202,6 +219,16 @@ function integrate(sc::CubMCG; resume::Union{Nothing, Dict{Symbol, Any}}=nothing
             solution = mu0
         end
         n_total = sc.n_init + max(n_mu, 0)
+        if sc.trace_iterations
+            push!(
+                log;
+                n=n_total,
+                solution=solution,
+                error_bound=bound_hw,
+                tol=sc.abs_tol,
+                elapsed=time() - t_start,
+            )
+        end
     else
         # Relative tolerance mode: iterative
         alphai = (sc.alpha - sc.alpha_sigma) / (2 * (1 - sc.alpha_sigma))
@@ -210,6 +237,16 @@ function integrate(sc::CubMCG; resume::Union{Nothing, Dict{Symbol, Any}}=nothing
         solution = mu0
         n_total = sc.n_init
         tau = 1.0
+        if sc.trace_iterations
+            push!(
+                log;
+                n=n_total,
+                solution=solution,
+                error_bound=bound_hw,
+                tol=tol0,
+                elapsed=time() - t_start,
+            )
+        end
 
         while true
             tol_eff = max(sc.abs_tol, sc.rel_tol * abs(solution))
@@ -233,6 +270,16 @@ function integrate(sc::CubMCG; resume::Union{Nothing, Dict{Symbol, Any}}=nothing
             y1 = evaluate(f, x1)
             solution = mean(y1)
             n_total += n_new
+            if sc.trace_iterations
+                push!(
+                    log;
+                    n=n_total,
+                    solution=solution,
+                    error_bound=bound_hw,
+                    tol=max(sc.abs_tol, sc.rel_tol * abs(solution)),
+                    elapsed=time() - t_start,
+                )
+            end
         end
     end
 
@@ -244,6 +291,9 @@ function integrate(sc::CubMCG; resume::Union{Nothing, Dict{Symbol, Any}}=nothing
         :bound_diff => 2 * bound_hw,
         :time_integrate => t_elapsed,
     )
+    if sc.trace_iterations
+        data[:iteration_log] = log
+    end
 
     return QMCResult(solution, data)
 end
diff --git a/src/stopping_criterion/cub_mlmc.jl b/src/stopping_criterion/cub_mlmc.jl
index 0cb08ae..f5934da 100644
--- a/src/stopping_criterion/cub_mlmc.jl
+++ b/src/stopping_criterion/cub_mlmc.jl
@@ -2,7 +2,8 @@
     CubMLMC(integrand::AbstractMLIntegrand; abs_tol=0.05, n_init=256,
             n_limit=10_000_000_000, alpha_ci=0.01,
             levels_min=2, levels_max=10,
-            alpha0=-1.0, beta0=-1.0, gamma0=-1.0)
+            alpha0=-1.0, beta0=-1.0, gamma0=-1.0,
+            trace_iterations=false)
 
 Multilevel Monte Carlo (MLMC) stopping criterion based on Giles (2008).
 
@@ -14,6 +15,9 @@ The rates α (weak error), β (variance decay), γ (cost growth) are
 estimated via linear regression if not provided (`alpha0`, `beta0`,
 `gamma0` ≤ 0).
 
+Set `trace_iterations=true` to record an `IterationLog` in
+`result.data[:iteration_log]`.
+
 # Arguments
 - `integrand`: A multilevel integrand (`AbstractMLIntegrand`).
 - `abs_tol`: Absolute error tolerance (converted to RMSE tolerance internally).
@@ -48,6 +52,7 @@ mutable struct CubMLMC{I <: AbstractMLIntegrand} <: AbstractStoppingCriterion
     alpha0::Float64
     beta0::Float64
     gamma0::Float64
+    trace_iterations::Bool
 end
 
 function CubMLMC(
@@ -62,6 +67,7 @@ function CubMLMC(
     alpha0::Float64=-1.0,
     beta0::Float64=-1.0,
     gamma0::Float64=-1.0,
+    trace_iterations::Bool=false,
 )
     levels_min >= 2 || throw(ArgumentError("levels_min must be ≥ 2"))
     levels_max >= levels_min || throw(ArgumentError("levels_max must be ≥ levels_min"))
@@ -84,6 +90,7 @@ function CubMLMC(
         alpha0,
         beta0,
         gamma0,
+        trace_iterations,
     )
 end
 
@@ -234,15 +241,57 @@ function _get_optimal_samples(sc::CubMLMC, state::_MLMCState)
     return ns
 end
 
+function _solution_mlmc(state::_MLMCState)
+    return sum(
+        state.sum_level[1, l + 1] / state.n_level[l + 1] for
+        l in 0:state.levels if state.n_level[l + 1] > 0
+    )
+end
+
+function _bias_estimate_mlmc(state::_MLMCState)
+    L = state.levels
+    L < 1 && return Inf
+    range_check = min(2, L - 1)
+    rem = 0.0
+    for k in 0:range_check
+        idx = L - k
+        rem = max(rem, state.mean_level[idx + 1] / 2.0^(k * state.alpha))
+    end
+    return rem / (2.0^state.alpha - 1)
+end
+
+function _varest_mlmc(state::_MLMCState)
+    return sum(
+        state.var_level[l + 1] / state.n_level[l + 1] for
+        l in 0:state.levels if state.n_level[l + 1] > 0
+    )
+end
+
+function _rmse_mlmc(sc::CubMLMC, state::_MLMCState)
+    bias = _bias_estimate_mlmc(state)
+    return sqrt(max(0.0, (1 - sc.theta) * _varest_mlmc(state) + sc.theta * bias^2))
+end
+
 # ── Main integrate method ──
 
 function integrate(sc::CubMLMC; resume::Union{Nothing, Dict{Symbol, Any}}=nothing)
     t_start = time()
     state = _init_mlmc_state(sc)
+    log = IterationLog()
 
     # Initial warm-up: n_init samples at each level
     diff_n = fill(sc.n_init, state.levels + 1)
     _update_samples!(sc, state, diff_n)
+    if sc.trace_iterations
+        push!(
+            log;
+            n=sum(state.n_level),
+            solution=_solution_mlmc(state),
+            error_bound=_rmse_mlmc(sc, state),
+            tol=sc.rmse_tol,
+            elapsed=time() - t_start,
+        )
+    end
 
     # Main loop
     while true
@@ -259,18 +308,21 @@ function integrate(sc::CubMLMC; resume::Union{Nothing, Dict{Symbol, Any}}=nothin
 
         # Take additional samples
         _update_samples!(sc, state, diff_n)
+        if sc.trace_iterations
+            push!(
+                log;
+                n=sum(state.n_level),
+                solution=_solution_mlmc(state),
+                error_bound=_rmse_mlmc(sc, state),
+                tol=sc.rmse_tol,
+                elapsed=time() - t_start,
+            )
+        end
 
         # Check if (almost) converged
         if all(diff_n .<= ceil.(Int, 0.01 .* state.n_level[1:(state.levels + 1)]))
             # Estimate remaining bias
-            L = state.levels
-            range_check = min(2, L - 1)
-            rem = 0.0
-            for k in 0:range_check
-                idx = L - k  # 0-based level
-                rem = max(rem, state.mean_level[idx + 1] / 2.0^(k * state.alpha))
-            end
-            rem /= (2.0^state.alpha - 1)
+            rem = _bias_estimate_mlmc(state)
 
             if rem > sqrt(sc.theta) * sc.rmse_tol
                 # Need more levels
@@ -312,6 +364,9 @@ function integrate(sc::CubMLMC; resume::Union{Nothing, Dict{Symbol, Any}}=nothin
         :rmse_tol => sc.rmse_tol,
         :time_integrate => t_elapsed,
     )
+    if sc.trace_iterations
+        data[:iteration_log] = log
+    end
 
     return QMCResult(solution, data)
 end
diff --git a/src/stopping_criterion/cub_mlmc_cont.jl b/src/stopping_criterion/cub_mlmc_cont.jl
index b7a7534..66ceb4b 100644
--- a/src/stopping_criterion/cub_mlmc_cont.jl
+++ b/src/stopping_criterion/cub_mlmc_cont.jl
@@ -2,7 +2,8 @@
     CubMLMCCont(integrand::AbstractMLIntegrand; abs_tol=0.05, n_init=256,
                 n_limit=10_000_000_000, alpha_ci=0.01,
                 levels_min=2, levels_max=10,
-                n_tols=10, inflate=100^(1/9), theta_init=0.5)
+                n_tols=10, inflate=100^(1/9), theta_init=0.5,
+                trace_iterations=false)
 
 Continuation multilevel Monte Carlo (MLMC) stopping criterion.
 
@@ -11,6 +12,9 @@ Runs the MLMC algorithm at a sequence of progressively tighter tolerances
 samples from coarser tolerances. This improves the estimation of the
 convergence rates α, β, γ before the final tolerance is reached.
 
+Set `trace_iterations=true` to record an `IterationLog` in
+`result.data[:iteration_log]`.
+
 # Arguments
 - `integrand`: A multilevel integrand (`AbstractMLIntegrand`).
 - `abs_tol`: Absolute error tolerance.
@@ -49,6 +53,7 @@ mutable struct CubMLMCCont{I <: AbstractMLIntegrand} <: AbstractStoppingCriterio
     alpha0::Float64
     beta0::Float64
     gamma0::Float64
+    trace_iterations::Bool
 end
 
 function CubMLMCCont(
@@ -63,6 +68,7 @@ function CubMLMCCont(
     n_tols::Int=10,
     inflate::Float64=100.0^(1/9),
     theta_init::Float64=0.5,
+    trace_iterations::Bool=false,
 )
     levels_min >= 2 || throw(ArgumentError("levels_min must be ≥ 2"))
     levels_max >= levels_min || throw(ArgumentError("levels_max must be ≥ levels_min"))
@@ -91,6 +97,7 @@ function CubMLMCCont(
         -1.0,
         -1.0,
         -1.0,
+        trace_iterations,
     )
 end
 
@@ -322,11 +329,22 @@ end
 function integrate(sc::CubMLMCCont; resume::Union{Nothing, Dict{Symbol, Any}}=nothing)
     t_start = time()
     state = _init_mlmc_cont_state(sc)
+    log = IterationLog()
 
     # Loop over progressively tighter tolerances
     for t in 0:(sc.n_tols - 1)
         step_tol = sc.inflate^(sc.n_tols - t - 1) * sc.target_tol
         _integrate_step!(sc, state, step_tol)
+        if sc.trace_iterations
+            push!(
+                log;
+                n=sum(state.n_level),
+                solution=_solution_mlmc(state),
+                error_bound=_rmse(sc, state, step_tol),
+                tol=step_tol,
+                elapsed=time() - t_start,
+            )
+        end
     end
 
     # Compute final solution
@@ -350,6 +368,9 @@ function integrate(sc::CubMLMCCont; resume::Union{Nothing, Dict{Symbol, Any}}=no
         :rmse_tol => sc.target_tol,
         :time_integrate => t_elapsed,
     )
+    if sc.trace_iterations
+        data[:iteration_log] = log
+    end
 
     return QMCResult(solution, data)
 end
diff --git a/src/stopping_criterion/cub_mlqmc.jl b/src/stopping_criterion/cub_mlqmc.jl
index c1e35dc..327e5ea 100644
--- a/src/stopping_criterion/cub_mlqmc.jl
+++ b/src/stopping_criterion/cub_mlqmc.jl
@@ -1,7 +1,7 @@
 """
     CubMLQMC(integrand::AbstractMLIntegrand; abs_tol=0.05, n_init=256,
              n_limit=10_000_000_000, alpha_ci=0.01,
-             levels_min=2, levels_max=10)
+             levels_min=2, levels_max=10, trace_iterations=false)
 
 Multilevel quasi-Monte Carlo (MLQMC) stopping criterion.
 
@@ -13,6 +13,9 @@ strategy — it targets the final tolerance directly.
 
 The discrete distribution must have `replications ≥ 4`.
 
+Set `trace_iterations=true` to record an `IterationLog` in
+`result.data[:iteration_log]`.
+
 # Arguments
 - `integrand`: A multilevel integrand (`AbstractMLIntegrand`).
 - `abs_tol`: Absolute error tolerance.
@@ -39,6 +42,7 @@ mutable struct CubMLQMC{I <: AbstractMLIntegrand} <: AbstractStoppingCriterion
     levels_min::Int
     levels_max::Int
     theta::Float64
+    trace_iterations::Bool
 end
 
 function Base.getproperty(sc::CubMLQMC, name::Symbol)
@@ -57,6 +61,7 @@ function CubMLQMC(
     alpha_ci::Float64=0.01,
     levels_min::Int=2,
     levels_max::Int=10,
+    trace_iterations::Bool=false,
 )
     levels_min >= 2 || throw(ArgumentError("levels_min must be ≥ 2"))
     levels_max >= levels_min || throw(ArgumentError("levels_max must be ≥ levels_min"))
@@ -75,7 +80,17 @@ function CubMLQMC(
         target_tol = rmse_tol
     end
 
-    return CubMLQMC(integrand, target_tol, n_init, n_limit, R, levels_min, levels_max, 0.5)
+    return CubMLQMC(
+        integrand,
+        target_tol,
+        n_init,
+        n_limit,
+        R,
+        levels_min,
+        levels_max,
+        0.5,
+        trace_iterations,
+    )
 end
 
 # Reuse _MLQMCState from cub_mlqmc_cont.jl
@@ -187,6 +202,7 @@ function integrate(sc::CubMLQMC; resume::Union{Nothing, Dict{Symbol, Any}}=nothi
     t_start = time()
     state = _init_mlqmc_direct_state(sc)
     target_tol = sc.target_tol
+    log = IterationLog()
 
     converged = false
     while !converged
@@ -222,11 +238,32 @@ function integrate(sc::CubMLQMC; resume::Union{Nothing, Dict{Symbol, Any}}=nothi
         end
 
         if converged
+            rmse = sqrt(max(0.0, (1 - sc.theta) * varest + sc.theta * state.bias_estimate^2))
+            if sc.trace_iterations
+                push!(
+                    log;
+                    n=(sc.replications * sum(state.n_level)),
+                    solution=sum(state.mean_level[1:state.levels]),
+                    error_bound=rmse,
+                    tol=target_tol,
+                    elapsed=time() - t_start,
+                )
+            end
             break
         end
 
         # Check overall convergence
         rmse = sqrt(max(0.0, (1 - sc.theta) * varest + sc.theta * state.bias_estimate^2))
+        if sc.trace_iterations
+            push!(
+                log;
+                n=(sc.replications * sum(state.n_level)),
+                solution=sum(state.mean_level[1:state.levels]),
+                error_bound=rmse,
+                tol=target_tol,
+                elapsed=time() - t_start,
+            )
+        end
         converged = rmse < target_tol
         if !converged
             if state.levels >= sc.levels_max
@@ -253,6 +290,9 @@ function integrate(sc::CubMLQMC; resume::Union{Nothing, Dict{Symbol, Any}}=nothi
         :rmse_tol => sc.target_tol,
         :time_integrate => t_elapsed,
     )
+    if sc.trace_iterations
+        data[:iteration_log] = log
+    end
 
     return QMCResult(solution, data)
 end
diff --git a/src/stopping_criterion/cub_mlqmc_cont.jl b/src/stopping_criterion/cub_mlqmc_cont.jl
index a3e8162..0dbb67f 100644
--- a/src/stopping_criterion/cub_mlqmc_cont.jl
+++ b/src/stopping_criterion/cub_mlqmc_cont.jl
@@ -2,7 +2,8 @@
     CubMLQMCCont(integrand::AbstractMLIntegrand; abs_tol=0.05, n_init=256,
                  n_limit=10_000_000_000, alpha_ci=0.01,
                  levels_min=2, levels_max=10,
-                 n_tols=10, inflate=100^(1/9), theta_init=0.5)
+                 n_tols=10, inflate=100^(1/9), theta_init=0.5,
+                 trace_iterations=false)
 
 Continuation multilevel quasi-Monte Carlo (MLQMC) stopping criterion.
 
@@ -14,6 +15,9 @@ tighter tolerances improves convergence rate estimation.
 
 The discrete distribution must have `replications ≥ 4`.
 
+Set `trace_iterations=true` to record an `IterationLog` in
+`result.data[:iteration_log]`.
+
 # Arguments
 - `integrand`: A multilevel integrand (`AbstractMLIntegrand`).
 - `abs_tol`: Absolute error tolerance.
@@ -50,6 +54,7 @@ mutable struct CubMLQMCCont{I <: AbstractMLIntegrand} <: AbstractStoppingCriteri
     inflate::Float64
     theta_init::Float64
     theta::Float64
+    trace_iterations::Bool
 end
 
 function CubMLQMCCont(
@@ -64,6 +69,7 @@ function CubMLQMCCont(
     n_tols::Int=10,
     inflate::Float64=100.0^(1/9),
     theta_init::Float64=0.5,
+    trace_iterations::Bool=false,
 )
     levels_min >= 2 || throw(ArgumentError("levels_min must be ≥ 2"))
     levels_max >= levels_min || throw(ArgumentError("levels_max must be ≥ levels_min"))
@@ -94,6 +100,7 @@ function CubMLQMCCont(
         inflate,
         theta_init,
         theta_init,
+        trace_iterations,
     )
 end
 
@@ -298,11 +305,22 @@ end
 function integrate(sc::CubMLQMCCont; resume::Union{Nothing, Dict{Symbol, Any}}=nothing)
     t_start = time()
     state = _init_mlqmc_state(sc)
+    log = IterationLog()
 
     # Continuation: loop over progressively tighter tolerances
     for t in 0:(sc.n_tols - 1)
         step_tol = sc.inflate^(sc.n_tols - t - 1) * sc.target_tol
         _integrate_step_qmc!(sc, state, step_tol)
+        if sc.trace_iterations
+            push!(
+                log;
+                n=(sc.replications * sum(state.n_level)),
+                solution=sum(state.mean_level[1:state.levels]),
+                error_bound=_rmse_qmc(sc, state),
+                tol=step_tol,
+                elapsed=time() - t_start,
+            )
+        end
     end
 
     # Final solution
@@ -321,6 +339,9 @@ function integrate(sc::CubMLQMCCont; resume::Union{Nothing, Dict{Symbol, Any}}=n
         :rmse_tol => sc.target_tol,
         :time_integrate => t_elapsed,
     )
+    if sc.trace_iterations
+        data[:iteration_log] = log
+    end
 
     return QMCResult(solution, data)
 end
diff --git a/src/stopping_criterion/cub_qmc_bayes_lattice_g.jl b/src/stopping_criterion/cub_qmc_bayes_lattice_g.jl
index 23c584d..dba0768 100644
--- a/src/stopping_criterion/cub_qmc_bayes_lattice_g.jl
+++ b/src/stopping_criterion/cub_qmc_bayes_lattice_g.jl
@@ -1,7 +1,8 @@
 """
     CubQMCBayesLatticeG(integrand; abs_tol=0.01, rel_tol=0.0,
                         n_init=2^8, n_max=2^22, order=2,
-                        ptransform=:C1SIN, errbd_type=:MLE, alpha=0.01)
+                        ptransform=:C1SIN, errbd_type=:MLE, alpha=0.01,
+                        trace_iterations=false)
 
 Bayesian QMC cubature for lattice rules with shift-invariant kernels.
 
@@ -9,6 +10,9 @@ Uses a Bernoulli-polynomial-based shift-invariant kernel diagonalized by FFT.
 The shape parameter θ is estimated via MLE or GCV. A periodization transform
 (default C1SIN) is applied to improve convergence for non-periodic integrands.
 
+Set `trace_iterations=true` to record an `IterationLog` in
+`result.data[:iteration_log]`.
+
 # Example
 ```julia
 dd = Lattice(3; randomize=true, seed=7)
@@ -28,6 +32,7 @@ mutable struct CubQMCBayesLatticeG{I <: AbstractIntegrand} <: AbstractStoppingCr
     ptransform::Symbol
     errbd_type::Symbol
     alpha::Float64
+    trace_iterations::Bool
 end
 
 function CubQMCBayesLatticeG(
@@ -40,6 +45,7 @@ function CubQMCBayesLatticeG(
     ptransform::Symbol=:C1SIN,
     errbd_type::Symbol=:MLE,
     alpha::Float64=0.01,
+    trace_iterations::Bool=false,
 )
     @assert ispow2(n_init) "n_init must be a power of 2"
     @assert ispow2(n_max) "n_max must be a power of 2"
@@ -56,6 +62,7 @@ function CubQMCBayesLatticeG(
         ptransform,
         errbd_type,
         alpha,
+        trace_iterations,
     )
 end
 
@@ -196,6 +203,7 @@ end
 # ── integrate ────────────────────────────────────────────────────────────────
 
 function integrate(sc::CubQMCBayesLatticeG; resume::Union{Nothing, Dict{Symbol, Any}}=nothing)
+    t_start = time()
     if resume !== nothing
         n_prev =
             haskey(resume, :n_per_rep) ? Int(resume[:n_per_rep]) :
@@ -207,6 +215,7 @@ function integrate(sc::CubQMCBayesLatticeG; resume::Union{Nothing, Dict{Symbol,
     mu_hat = 0.0;
     err = Inf;
     n_iter = 0
+    log = IterationLog()
 
     while n <= sc.n_max
         n_iter += 1
@@ -225,6 +234,9 @@ function integrate(sc::CubQMCBayesLatticeG; resume::Union{Nothing, Dict{Symbol,
             _bayes_lattice_stop(x_uniform, ftilde, n, sc.order, sc.errbd_type, sc.alpha)
 
         tol = max(sc.abs_tol, sc.rel_tol * abs(mu_hat))
+        if sc.trace_iterations
+            push!(log; n=n, solution=mu_hat, error_bound=err, tol=tol, elapsed=time() - t_start)
+        end
         err <= tol && break
 
         2n > sc.n_max &&
@@ -243,6 +255,9 @@ function integrate(sc::CubQMCBayesLatticeG; resume::Union{Nothing, Dict{Symbol,
         :ptransform => sc.ptransform,
         :errbd_type => sc.errbd_type,
     )
+    if sc.trace_iterations
+        data[:iteration_log] = log
+    end
     return QMCResult(mu_hat, data)
 end
 
diff --git a/src/stopping_criterion/cub_qmc_bayes_net_g.jl b/src/stopping_criterion/cub_qmc_bayes_net_g.jl
index ba75d91..661867d 100644
--- a/src/stopping_criterion/cub_qmc_bayes_net_g.jl
+++ b/src/stopping_criterion/cub_qmc_bayes_net_g.jl
@@ -1,11 +1,15 @@
 """
     CubQMCBayesNetG(integrand; abs_tol=0.01, rel_tol=0.0,
                     n_init=2^8, n_max=2^22, order=2,
-                    ptransform=:NONE, errbd_type=:MLE, alpha=0.01)
+                    ptransform=:NONE, errbd_type=:MLE, alpha=0.01,
+                    trace_iterations=false)
 
 Bayesian QMC cubature for digital net (Sobol') rules with digitally
 shift-invariant kernels, diagonalized by the Walsh-Hadamard Transform.
 
+Set `trace_iterations=true` to record an `IterationLog` in
+`result.data[:iteration_log]`.
+
 # Example
 ```julia
 dd = DigitalNetB2(3; randomize="LMS_DS", seed=7)
@@ -25,6 +29,7 @@ mutable struct CubQMCBayesNetG{I <: AbstractIntegrand} <: AbstractStoppingCriter
     ptransform::Symbol
     errbd_type::Symbol
     alpha::Float64
+    trace_iterations::Bool
 end
 
 function CubQMCBayesNetG(
@@ -37,6 +42,7 @@ function CubQMCBayesNetG(
     ptransform::Symbol=:NONE,
     errbd_type::Symbol=:MLE,
     alpha::Float64=0.01,
+    trace_iterations::Bool=false,
 )
     @assert ispow2(n_init) "n_init must be a power of 2"
     @assert ispow2(n_max) "n_max must be a power of 2"
@@ -53,6 +59,7 @@ function CubQMCBayesNetG(
         ptransform,
         errbd_type,
         alpha,
+        trace_iterations,
     )
 end
 
@@ -181,6 +188,7 @@ end
 # ── integrate ────────────────────────────────────────────────────────────────
 
 function integrate(sc::CubQMCBayesNetG; resume::Union{Nothing, Dict{Symbol, Any}}=nothing)
+    t_start = time()
     if resume !== nothing
         n_prev =
             haskey(resume, :n_per_rep) ? Int(resume[:n_per_rep]) :
@@ -192,6 +200,7 @@ function integrate(sc::CubQMCBayesNetG; resume::Union{Nothing, Dict{Symbol, Any}
     mu_hat = 0.0;
     err = Inf;
     n_iter = 0
+    log = IterationLog()
 
     while n <= sc.n_max
         n_iter += 1
@@ -208,6 +217,9 @@ function integrate(sc::CubQMCBayesNetG; resume::Union{Nothing, Dict{Symbol, Any}
         mu_hat, err = _bayes_net_stop(x_uniform, ftilde, n, sc.order, sc.errbd_type, sc.alpha)
 
         tol = max(sc.abs_tol, sc.rel_tol * abs(mu_hat))
+        if sc.trace_iterations
+            push!(log; n=n, solution=mu_hat, error_bound=err, tol=tol, elapsed=time() - t_start)
+        end
         err <= tol && break
 
         2n > sc.n_max &&
@@ -226,6 +238,9 @@ function integrate(sc::CubQMCBayesNetG; resume::Union{Nothing, Dict{Symbol, Any}
         :ptransform => sc.ptransform,
         :errbd_type => sc.errbd_type,
     )
+    if sc.trace_iterations
+        data[:iteration_log] = log
+    end
     return QMCResult(mu_hat, data)
 end
 
diff --git a/src/true_measure/student_t.jl b/src/true_measure/student_t.jl
index a5fa973..eb5d60e 100644
--- a/src/true_measure/student_t.jl
+++ b/src/true_measure/student_t.jl
@@ -6,13 +6,17 @@ Multivariate Student-t measure via sequential conditioning.
 Each dimension uses the inverse CDF of the marginal t-distribution,
 applied to the uniform samples.
 
+The exact `df=1.0` (Cauchy) and `df=2.0` cases use closed-form quantiles,
+with a dedicated standard-parameter fast path for `loc=0.0, scale=1.0`.
+
 # Arguments
 - `dd`: discrete distribution.
 - `df`: degrees of freedom (> 0).
 - `loc`: location (scalar or d-vector).
 - `scale`: scale (scalar or d-vector, positive).
 """
-struct StudentT{D <: AbstractDiscreteDistribution} <: AbstractTrueMeasure
+struct StudentT{D <: AbstractDiscreteDistribution, Standard, QuantileCase} <:
+       AbstractTrueMeasure
     dd::D
     dimension::Int
     df::Float64
@@ -28,35 +32,153 @@ function StudentT(dd::AbstractDiscreteDistribution; df::Float64=2.0, loc=0.0, sc
     all(scale_vec .> 0) || throw(ArgumentError("scale must be positive"))
     length(loc_vec) == d || throw(ArgumentError("loc length must match dimension"))
     length(scale_vec) == d || throw(ArgumentError("scale length must match dimension"))
-    return StudentT(dd, d, df, loc_vec, scale_vec)
+    is_standard = all(iszero, loc_vec) && all(==(1.0), scale_vec)
+    quantile_case = df == 1.0 ? :df1 : (df == 2.0 ? :df2 : :generic)
+    return StudentT{typeof(dd), is_standard, quantile_case}(dd, d, df, loc_vec, scale_vec)
 end
 
-function transform(tm::StudentT, x::AbstractMatrix)
+@inline _student_t_quantile_df1(u::Real) = tanpi(_open_unit_interval(u) - 0.5)
+
+@inline function _student_t_quantile_df2(u::Real)
+    p = _open_unit_interval(u)
+    return (2.0 * p - 1.0) / sqrt(2.0 * p * (1.0 - p))
+end
+
+function _transform_student_t_df1_standard!(y::Matrix{Float64}, x::AbstractMatrix)
     n, d = size(x)
-    y = Matrix{Float64}(undef, n, d)
-    if tm.df == 2.0
-        # Exact closed form for the Student-t quantile at ν = 2:
-        #   Q(p) = (2p − 1) / √(2p(1 − p)).
-        # This avoids the per-element inverse-incomplete-beta in
-        # Distributions.quantile and matches it to numerical precision.
-        @inbounds for j in 1:d
-            loc = tm.loc[j]
-            scale = tm.scale[j]
-            @simd for i in 1:n
-                p = _open_unit_interval(x[i, j])
-                y[i, j] = loc + scale * (2.0 * p - 1.0) / sqrt(2.0 * p * (1.0 - p))
-            end
+    @inbounds for j in 1:d
+        @simd for i in 1:n
+            y[i, j] = _student_t_quantile_df1(x[i, j])
         end
-    else
-        tdist = Distributions.TDist(tm.df)
-        @inbounds for j in 1:d
-            for i in 1:n
-                y[i, j] =
-                    tm.loc[j] + tm.scale[j] * quantile(tdist, _open_unit_interval(x[i, j]))
-            end
+    end
+    return y
+end
+
+function _transform_student_t_df1_affine!(y::Matrix{Float64}, tm::StudentT, x::AbstractMatrix)
+    n, d = size(x)
+    @inbounds for j in 1:d
+        loc = tm.loc[j]
+        scale = tm.scale[j]
+        @simd for i in 1:n
+            y[i, j] = loc + scale * _student_t_quantile_df1(x[i, j])
         end
     end
     return y
 end
 
+function _transform_student_t_df2_standard!(y::Matrix{Float64}, x::AbstractMatrix)
+    n, d = size(x)
+    @inbounds for j in 1:d
+        @simd for i in 1:n
+            y[i, j] = _student_t_quantile_df2(x[i, j])
+        end
+    end
+    return y
+end
+
+function _transform_student_t_df2_affine!(y::Matrix{Float64}, tm::StudentT, x::AbstractMatrix)
+    n, d = size(x)
+    @inbounds for j in 1:d
+        loc = tm.loc[j]
+        scale = tm.scale[j]
+        @simd for i in 1:n
+            y[i, j] = loc + scale * _student_t_quantile_df2(x[i, j])
+        end
+    end
+    return y
+end
+
+function _transform_student_t_quantile_standard!(
+    y::Matrix{Float64},
+    tm::StudentT,
+    x::AbstractMatrix,
+)
+    n, d = size(x)
+    tdist = Distributions.TDist(tm.df)
+    @inbounds for j in 1:d
+        for i in 1:n
+            y[i, j] = quantile(tdist, _open_unit_interval(x[i, j]))
+        end
+    end
+    return y
+end
+
+function _transform_student_t_quantile_affine!(
+    y::Matrix{Float64},
+    tm::StudentT,
+    x::AbstractMatrix,
+)
+    n, d = size(x)
+    tdist = Distributions.TDist(tm.df)
+    @inbounds for j in 1:d
+        loc = tm.loc[j]
+        scale = tm.scale[j]
+        for i in 1:n
+            y[i, j] = loc + scale * quantile(tdist, _open_unit_interval(x[i, j]))
+        end
+    end
+    return y
+end
+
+function transform(tm::StudentT, x::AbstractMatrix)
+    n, d = size(x)
+    y = Matrix{Float64}(undef, n, d)
+    _transform_student_t!(y, tm, x)
+    return y
+end
+
+function _transform_student_t!(
+    y::Matrix{Float64},
+    tm::StudentT{<:Any, true, :df1},
+    x::AbstractMatrix,
+)
+    # Exact Cauchy quantile for ν = 1:
+    #   Q(p) = tan(π(p − 1/2)).
+    return _transform_student_t_df1_standard!(y, x)
+end
+
+function _transform_student_t!(
+    y::Matrix{Float64},
+    tm::StudentT{<:Any, false, :df1},
+    x::AbstractMatrix,
+)
+    return _transform_student_t_df1_affine!(y, tm, x)
+end
+
+function _transform_student_t!(
+    y::Matrix{Float64},
+    tm::StudentT{<:Any, true, :df2},
+    x::AbstractMatrix,
+)
+    # Exact closed form for the Student-t quantile at ν = 2:
+    #   Q(p) = (2p − 1) / √(2p(1 − p)).
+    # This avoids the per-element inverse-incomplete-beta in
+    # Distributions.quantile and matches it to numerical precision.
+    return _transform_student_t_df2_standard!(y, x)
+end
+
+function _transform_student_t!(
+    y::Matrix{Float64},
+    tm::StudentT{<:Any, false, :df2},
+    x::AbstractMatrix,
+)
+    return _transform_student_t_df2_affine!(y, tm, x)
+end
+
+function _transform_student_t!(
+    y::Matrix{Float64},
+    tm::StudentT{<:Any, true, :generic},
+    x::AbstractMatrix,
+)
+    return _transform_student_t_quantile_standard!(y, tm, x)
+end
+
+function _transform_student_t!(
+    y::Matrix{Float64},
+    tm::StudentT{<:Any, false, :generic},
+    x::AbstractMatrix,
+)
+    return _transform_student_t_quantile_affine!(y, tm, x)
+end
+
 Base.show(io::IO, tm::StudentT) = print(io, "StudentT(d=$(tm.dimension), df=$(tm.df))")
diff --git a/test/README.md b/test/README.md
index 8d66cec..c9c6997 100644
--- a/test/README.md
+++ b/test/README.md
@@ -9,14 +9,14 @@ This directory contains the automated test suite for QMC.jl.
 
 ## Test files
 
-- `test_discrete_distributions.jl`
-- `test_true_measures.jl`
-- `test_integrands.jl`
-- `test_kernels.jl`
-- `test_stopping_criteria.jl`
-- `test_integration.jl`
-- `test_multilevel.jl`
-- `test_aqua.jl`
+- `test_discrete_distributions.jl`: tests discrete distribution constructions and sampling behavior.
+- `test_true_measures.jl`: validates true-measure definitions and related properties.
+- `test_integrands.jl`: checks integrand setup, transformations, and evaluations.
+- `test_kernels.jl`: verifies kernel implementations and expected numerical behavior.
+- `test_stopping_criteria.jl`: tests stopping rules and convergence-trigger logic.
+- `test_integration.jl`: exercises end-to-end integration workflows and outputs.
+- `test_multilevel.jl`: covers multilevel algorithms and level-coupling behavior.
+- `test_aqua.jl`: runs Aqua automatic quality checks.
 
 ## Local commands
 
diff --git a/test/runtests.jl b/test/runtests.jl
index 4b67798..dba1883 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -31,7 +31,14 @@ times = Dict{String, Float64}()
 
 @testset "QMC" begin
     for test_file in TEST_FILES
-        elapsed = @elapsed include(test_file)
+        elapsed = @elapsed begin
+            try
+                include(test_file)
+            catch err
+                println("  ✗ failed in $(test_file)")
+                rethrow(err)
+            end
+        end
         times[test_file] = elapsed
         println("  ⏱  $(test_file): $(fmt_duration(elapsed))")
     end
diff --git a/test/test_aqua.jl b/test/test_aqua.jl
index 4d3bf33..6edc1ad 100644
--- a/test/test_aqua.jl
+++ b/test/test_aqua.jl
@@ -5,7 +5,7 @@ using QMC
     Aqua.test_all(
         QMC;
         ambiguities=false,       # skip for now — many methods with Union types
-        stale_deps=(ignore=[:NBInclude],),  # used by test/run_notebooks.jl, not src/
+        stale_deps=(ignore=[:NBInclude, :Downloads],),  # NBInclude is test-only; Downloads is loaded dynamically for remote LDData fetches.
         deps_compat=true,
         piracies=false,          # we extend Base.push! for IterationLog
         project_extras=true,
diff --git a/test/test_discrete_distributions.jl b/test/test_discrete_distributions.jl
index d3faa15..1501c5a 100644
--- a/test/test_discrete_distributions.jl
+++ b/test/test_discrete_distributions.jl
@@ -93,8 +93,32 @@
             dd_direct =
                 Lattice(3; randomize=false, order="linear", generating_vector=[1, 3, 5, 7])
             @test gen_samples(dd_file, 8) ≈ gen_samples(dd_direct, 8)
+            @test dd_file.n_limit == 16
             @test_throws ArgumentError gen_samples(dd_file, 17)
         end
+        mktemp() do path, io
+            write(io, "# d_limit\n2\n# n_limit\n16\n1\n3\n5\n7\n")
+            close(io)
+            @test_throws ArgumentError Lattice(
+                3;
+                randomize=false,
+                order="linear",
+                generating_vector=path,
+            )
+        end
+
+        dd_default = Lattice(3; randomize=false, order="linear")
+        for ref in (
+            "kuo.lattice-33002-1024-1048576.9125",
+            "lattice/kuo.lattice-33002-1024-1048576.9125.txt",
+            "https://github.com/QMCSoftware/LDData/tree/main/lattice/kuo.lattice-33002-1024-1048576.9125.txt",
+            "https://raw.githubusercontent.com/QMCSoftware/LDData/main/lattice/kuo.lattice-33002-1024-1048576.9125.txt",
+        )
+            dd_named = Lattice(3; randomize=false, order="linear", generating_vector=ref)
+            @test dd_named.gen_vector == dd_default.gen_vector
+            @test dd_named.n_limit == dd_default.n_limit
+            @test gen_samples(dd_named, 16) ≈ gen_samples(dd_default, 16)
+        end
 
         dd_default = Lattice(3; randomize=false, order="linear")
         for ref in (
@@ -178,12 +202,181 @@
         @test DigitalNetB2(2; randomize="LMS DS").randomize == "LMS_DS"
         @test DigitalNetB2(2; randomize=:LMS_DS).randomize == "LMS_DS"
         @test DigitalNetB2(2; randomize="false").randomize == "none"
+        @test DigitalNetB2(2; randomize="OWEN").randomize == "NUS"
+
+        # Constructor parity: accept QMCPy-style order spellings.
+        @test DigitalNetB2(2; order="GRAY CODE").graycode
+        @test !DigitalNetB2(2; order="NATURAL").graycode
+        @test !DigitalNetB2(2; order=:RADICAL_INVERSE).graycode
+        let a = DigitalNetB2(3; randomize="none", seed=7, order="natural"),
+            b = DigitalNetB2(3; randomize="none", seed=7, graycode=false)
+
+            @test gen_samples(a, 16) ≈ gen_samples(b, 16)
+        end
+        @test_throws ArgumentError DigitalNetB2(2; order="bogus")
+        @test_throws ArgumentError DigitalNetB2(2; order="gray", graycode=false)
 
         # Custom direction matrices keep the existing generation path.
         V8 = Int.(dd.direction_nums[:, 1:8] .>> 24)
         dd_custom = DigitalNetB2(2; randomize="none", seed=1, generating_matrices=V8)
         @test gen_samples(dd_custom, 16) ≈ gen_samples(dd, 16)
-        @test dd_custom.direction_nums == UInt32.(V8)
+        @test dd_custom.direction_nums == UInt64.(V8)
+        V8_lsb = map(v -> Int(bitreverse(UInt32(v)) >> 24), V8)
+        dd_custom_lsb =
+            DigitalNetB2(2; randomize="none", seed=1, generating_matrices=V8_lsb, msb=false)
+        @test gen_samples(dd_custom_lsb, 16) ≈ gen_samples(dd_custom, 16)
+        Vsparse = [1 2 4 8 1 2 4 8; 1 2 4 8 1 2 4 8]
+        dd_sparse_lsb =
+            DigitalNetB2(2; randomize="none", seed=1, generating_matrices=Vsparse, msb=false)
+        @test dd_sparse_lsb.direction_nums == map(v -> bitreverse(UInt64(v)) >> 56, Vsparse)
+
+        # Advanced constructor parity: widened t-bits and msb handling.
+        @test DigitalNetB2(2; randomize="none", t=40).t == 40
+        @test gen_samples(DigitalNetB2(2; randomize="none", seed=1, t=40), 16) == x
+        xrep_t =
+            gen_samples(DigitalNetB2(2; randomize="none", seed=1, replications=3, t=40), 16)
+        @test size(xrep_t) == (3, 16, 2)
+        @test all(0.0 .<= xrep_t .< 1.0)
+        @test xrep_t[1, :, :] == xrep_t[2, :, :]
+        @test xrep_t[2, :, :] == xrep_t[3, :, :]
+
+        x_nus_t_a = gen_samples(DigitalNetB2(2; randomize="NUS", seed=9, t=40), 16)
+        x_nus_t_b = gen_samples(DigitalNetB2(2; randomize="NUS", seed=9, t=40), 16)
+        @test size(x_nus_t_a) == (16, 2)
+        @test all(0.0 .<= x_nus_t_a .< 1.0)
+        @test x_nus_t_a ≈ x_nus_t_b
+        x_nus_rep_t =
+            gen_samples(DigitalNetB2(2; randomize="NUS", seed=9, replications=2, t=40), 8)
+        @test size(x_nus_rep_t) == (2, 8, 2)
+        @test all(0.0 .<= x_nus_rep_t .< 1.0)
+
+        # Higher-order / interlaced digital nets.
+        dd_alpha = DigitalNetB2(3; randomize="none", seed=7, alpha=2, graycode=false)
+        @test dd_alpha.alpha == 2
+        @test dd_alpha.t == 64
+        @test gen_samples(dd_alpha, 4) == [
+            0.0 0.0 0.0
+            0.75 0.75 0.75
+            0.4375 0.9375 0.1875
+            0.6875 0.1875 0.9375
+        ]
+
+        x_alpha_lmsds_a = gen_samples(DigitalNetB2(3; randomize="LMS_DS", seed=7, alpha=2), 8)
+        x_alpha_lmsds_b = gen_samples(DigitalNetB2(3; randomize="LMS_DS", seed=7, alpha=2), 8)
+        @test size(x_alpha_lmsds_a) == (8, 3)
+        @test all(0.0 .<= x_alpha_lmsds_a .< 1.0)
+        @test x_alpha_lmsds_a ≈ x_alpha_lmsds_b
+
+        x_alpha_lms_a = gen_samples(DigitalNetB2(3; randomize="LMS", seed=7, alpha=2), 8)
+        x_alpha_lms_b = gen_samples(DigitalNetB2(3; randomize="LMS", seed=7, alpha=2), 8)
+        @test size(x_alpha_lms_a) == (8, 3)
+        @test all(0.0 .<= x_alpha_lms_a .< 1.0)
+        @test x_alpha_lms_a ≈ x_alpha_lms_b
+        x_alpha_lms_rep =
+            gen_samples(DigitalNetB2(3; randomize="LMS", seed=7, replications=2, alpha=2), 8)
+        @test size(x_alpha_lms_rep) == (2, 8, 3)
+        @test all(0.0 .<= x_alpha_lms_rep .< 1.0)
+
+        x_alpha_ds_a = gen_samples(DigitalNetB2(3; randomize="DS", seed=7, alpha=2), 8)
+        x_alpha_ds_b = gen_samples(DigitalNetB2(3; randomize="DS", seed=7, alpha=2), 8)
+        @test size(x_alpha_ds_a) == (8, 3)
+        @test all(0.0 .<= x_alpha_ds_a .< 1.0)
+        @test x_alpha_ds_a ≈ x_alpha_ds_b
+        x_alpha_ds_rep =
+            gen_samples(DigitalNetB2(3; randomize="DS", seed=7, replications=2, alpha=2), 8)
+        @test size(x_alpha_ds_rep) == (2, 8, 3)
+        @test all(0.0 .<= x_alpha_ds_rep .< 1.0)
+
+        x_alpha_nus_a = gen_samples(DigitalNetB2(3; randomize="NUS", seed=7, alpha=2), 8)
+        x_alpha_nus_b = gen_samples(DigitalNetB2(3; randomize="NUS", seed=7, alpha=2), 8)
+        @test size(x_alpha_nus_a) == (8, 3)
+        @test all(0.0 .<= x_alpha_nus_a .< 1.0)
+        @test x_alpha_nus_a ≈ x_alpha_nus_b
+        x_alpha_nus_rep =
+            gen_samples(DigitalNetB2(3; randomize="NUS", seed=7, replications=2, alpha=2), 8)
+        @test size(x_alpha_nus_rep) == (2, 8, 3)
+        @test all(0.0 .<= x_alpha_nus_rep .< 1.0)
+        if QMC._HAS_DNB2_FUSED[]
+            old_fused = QMC._HAS_DNB2_FUSED[]
+            try
+                QMC._HAS_DNB2_FUSED[] = true
+                x_alpha_fused =
+                    gen_samples(DigitalNetB2(3; randomize="LMS_DS", seed=7, alpha=2), 8)
+                QMC._HAS_DNB2_FUSED[] = false
+                x_alpha_unfused =
+                    gen_samples(DigitalNetB2(3; randomize="LMS_DS", seed=7, alpha=2), 8)
+                @test x_alpha_fused ≈ x_alpha_unfused
+            finally
+                QMC._HAS_DNB2_FUSED[] = old_fused
+            end
+        end
+
+        # LDData-style text sources: local file paths and QMCPy-compatible names.
+        mktemp() do path, io
+            write(io, "# base\n2\n# d_limit\n2\n# n_limit\n16\n# bit precision\n32\n")
+            for j in 1:2
+                write(io, join(string.(Int.(dd.direction_nums[j, :])), ' '))
+                write(io, "\n")
+            end
+            flush(io)
+
+            dd_file = DigitalNetB2(2; randomize="none", seed=1, generating_matrices=path)
+            @test gen_samples(dd_file, 16) ≈ gen_samples(dd, 16)
+            @test dd_file.n_limit == 16
+            @test_throws ArgumentError gen_samples(dd_file, 17)
+        end
+        mktemp() do path, io
+            write(io, "# base\n2\n# d_limit\n2\n# n_limit\n16\n# bit precision\n64\n")
+            for _ in 1:2
+                write(io, "1 2 4 8\n")
+            end
+            flush(io)
+
+            dd_file64 =
+                DigitalNetB2(1; randomize="none", seed=1, alpha=2, generating_matrices=path)
+            @test dd_file64.source_bits == 64
+            @test dd_file64.t == 64
+            @test size(gen_samples(dd_file64, 4)) == (4, 1)
+        end
+        mktemp() do path, io
+            write(io, "# base\n2\n# d_limit\n1024\n# n_limit\n2\n# bit precision\n32\n")
+            for _ in 1:1024
+                write(io, "1\n")
+            end
+            flush(io)
+
+            @test_throws ArgumentError DigitalNetB2(
+                513;
+                randomize="none",
+                alpha=2,
+                generating_matrices=path,
+            )
+        end
+        mktemp() do path, io
+            write(io, "# base\n2\n# d_limit\n1026\n# n_limit\n2\n# bit precision\n32\n")
+            for _ in 1:1026
+                write(io, "1\n")
+            end
+            flush(io)
+
+            dd_file_alpha_limit =
+                DigitalNetB2(513; randomize="none", alpha=2, generating_matrices=path)
+            @test dd_file_alpha_limit.source_bits == 32
+            @test dd_file_alpha_limit.n_limit == 2
+            @test size(gen_samples(dd_file_alpha_limit, 2)) == (2, 513)
+        end
+
+        dd_default = DigitalNetB2(3; randomize="none", seed=7)
+        for ref in (
+            "joe_kuo.6.21201.txt",
+            "joe_kuo.6.1024.txt",
+            "dnet/joe_kuo.6.21201.txt",
+            "https://github.com/QMCSoftware/LDData/tree/main/dnet/joe_kuo.6.21201.txt",
+            "https://raw.githubusercontent.com/QMCSoftware/LDData/main/dnet/joe_kuo.6.21201.txt",
+        )
+            dd_named = DigitalNetB2(3; randomize="none", seed=7, generating_matrices=ref)
+            @test gen_samples(dd_named, 16) ≈ gen_samples(dd_default, 16)
+        end
 
         # The custom matrix precision limits the allowable sample window.
         dd_short = DigitalNetB2(
@@ -195,16 +388,34 @@
         @test_throws ArgumentError gen_samples(dd_short, 9)
         @test_throws ArgumentError gen_samples(dd_short, 4; n_start=5)
 
+        # Higher-order bundled data uses raw Sobol' dimensions, so the default
+        # Joe-Kuo table tops out at floor(1024 / alpha) output dimensions.
+        dd_alpha_limit = DigitalNetB2(512; randomize="none", alpha=2)
+        @test size(gen_samples(dd_alpha_limit, 2)) == (2, 512)
+        @test_throws ArgumentError DigitalNetB2(513; randomize="none", alpha=2)
+
+        # Explicit custom matrices can still exceed the bundled effective-dimension
+        # limit as long as enough raw rows are supplied.
+        dd_alpha_custom =
+            DigitalNetB2(513; randomize="none", alpha=2, generating_matrices=ones(Int, 1026, 1))
+        @test size(gen_samples(dd_alpha_custom, 2)) == (2, 513)
+
         @test_throws ArgumentError DigitalNetB2(2; generating_matrices=ones(Int, 1, 4))
         @test_throws ArgumentError DigitalNetB2(2; generating_matrices=ones(Int, 2, 33))
         @test_throws ArgumentError DigitalNetB2(2; generating_matrices=zeros(Int, 2, 4))
+        @test_throws ArgumentError DigitalNetB2(2; generating_matrices=fill(1 << 20, 2, 8))
+        @test_throws ArgumentError DigitalNetB2(2; randomize="none", t=31)
+        @test_throws ArgumentError DigitalNetB2(2; randomize="none", t=65)
+        @test_throws ArgumentError DigitalNetB2(2; generating_matrices=V8, t=7)
         @test_throws ArgumentError DigitalNetB2(
             2;
-            generating_matrices=Int.(dd.direction_nums[:, 1:8]),
+            generating_matrices=fill(BigInt(typemax(UInt64)) + 1, 2, 4),
         )
         @test_throws ArgumentError DigitalNetB2(
-            2;
-            generating_matrices=fill(BigInt(typemax(UInt32)) + 1, 2, 4),
+            3;
+            randomize="none",
+            alpha=2,
+            generating_matrices=V8,
         )
     end
 
diff --git a/test/test_integrands.jl b/test/test_integrands.jl
index e8281bf..0fc238a 100644
--- a/test/test_integrands.jl
+++ b/test/test_integrands.jl
@@ -348,8 +348,8 @@
         xc = transform(tmc, gen_samples(ddc, 2^16))
         closed, total = compute_sensitivity_indices(sic, xc)
         ref = ishigami_exact()
-        @test all(abs.(closed .- ref.closed) .<= 0.05)
-        @test all(abs.(total .- ref.total) .<= 0.05)
+        @test all(abs.(closed .- ref.closed) .<= 0.06)
+        @test all(abs.(total .- ref.total) .<= 0.06)
     end
 
     @testset "BayesianLRCoeffs" begin
diff --git a/test/test_kernels.jl b/test/test_kernels.jl
index 86750c6..3b8af01 100644
--- a/test/test_kernels.jl
+++ b/test/test_kernels.jl
@@ -39,7 +39,7 @@
         K = kernel_matrix(k, x)
         @test size(K) == (3, 3)
         @test K ≈ K'  # symmetric
-        @test all(eigvals(K) .>= -1e-10)  # PSD
+        @test all(eigvals(K) .>= -1e-8)  # PSD up to numerical roundoff
     end
 
     @testset "KernelMatern32" begin
@@ -71,7 +71,7 @@
         K = kernel_matrix(k, x)
         @test size(K) == (3, 3)
         @test K ≈ K'
-        @test all(eigvals(K) .>= -1e-10)
+        @test all(eigvals(K) .>= -1e-8)
         @test_throws ArgumentError KernelRationalQuadratic(alpha=0.0)
         @test_throws ArgumentError KernelRationalQuadratic(lengthscale=-1.0)
     end
diff --git a/test/test_multilevel.jl b/test/test_multilevel.jl
index 5dd7c69..276f0cd 100644
--- a/test/test_multilevel.jl
+++ b/test/test_multilevel.jl
@@ -93,15 +93,52 @@ end
         dd2 = spawn_dd(dd, 8)
         @test dd2.dimension == 8
 
-        dd_lat = Lattice(4; replications=8)
+        dd_lat = Lattice(4; replications=8, order="gray")
         dd_lat2 = spawn_dd(dd_lat, 16)
         @test dd_lat2.dimension == 16
+        @test dd_lat2.order == "gray"
         @test dd_lat2.replications == 8
 
+        lat_source = [1, 3, 5, 7, 9, 11]
+        dd_lat_custom =
+            Lattice(4; randomize=false, order="linear", generating_vector=lat_source)
+        dd_lat_custom2 = spawn_dd(dd_lat_custom, 6)
+        @test dd_lat_custom2.dimension == 6
+        @test dd_lat_custom2.order == "linear"
+        @test dd_lat_custom2.gen_vector == UInt64.(lat_source)
+        @test gen_samples(dd_lat_custom2, 8) == gen_samples(
+            Lattice(6; randomize=false, order="linear", generating_vector=lat_source),
+            8,
+        )
+        @test_throws ArgumentError spawn_dd(dd_lat_custom, 7)
+
+        mktemp() do path, io
+            write(io, "# d_limit\n6\n# n_limit\n16\n1\n3\n5\n7\n9\n11\n")
+            close(io)
+            dd_lat_file = Lattice(4; randomize=false, order="linear", generating_vector=path)
+            dd_lat_file2 = spawn_dd(dd_lat_file, 6)
+            @test dd_lat_file2.dimension == 6
+            @test dd_lat_file2.order == "linear"
+            @test dd_lat_file2.n_limit == 16
+            @test dd_lat_file2.gen_vector == UInt64.(lat_source)
+            @test gen_samples(dd_lat_file2, 8) == gen_samples(
+                Lattice(6; randomize=false, order="linear", generating_vector=lat_source),
+                8,
+            )
+        end
+
         dd_dn = DigitalNetB2(4; replications=8)
         dd_dn2 = spawn_dd(dd_dn, 16)
         @test dd_dn2.dimension == 16
 
+        dd_dn_alpha = DigitalNetB2(4; randomize="LMS_DS", seed=7, replications=2, alpha=2)
+        dd_dn_alpha2 = spawn_dd(dd_dn_alpha, 8)
+        @test dd_dn_alpha2.dimension == 8
+        @test dd_dn_alpha2.alpha == 2
+        @test dd_dn_alpha2.t == dd_dn_alpha.t
+        @test dd_dn_alpha2.replications == 2
+        @test size(gen_samples(dd_dn_alpha2, 8)) == (2, 8, 8)
+
         dd_dn_custom = DigitalNetB2(
             4;
             randomize="none",
@@ -109,6 +146,21 @@ end
         )
         dd_dn_custom2 = spawn_dd(dd_dn_custom, 2)
         @test dd_dn_custom2.direction_nums == dd_dn_custom.direction_nums[1:2, :]
+
+        mktemp() do path, io
+            write(io, "2\n4\n16\n64\n")
+            for _ in 1:4
+                write(io, "1 2 4 8\n")
+            end
+            flush(io)
+
+            dd_dn_file = DigitalNetB2(1; randomize="none", alpha=2, generating_matrices=path)
+            dd_dn_file2 = spawn_dd(dd_dn_file, 1)
+            @test dd_dn_file2.source_bits == dd_dn_file.source_bits
+            @test dd_dn_file2.t == dd_dn_file.t
+            @test dd_dn_file2.alpha == dd_dn_file.alpha
+            @test gen_samples(dd_dn_file2, 4) == gen_samples(dd_dn_file, 4)
+        end
     end
 
     @testset "spawn_tm" begin
@@ -139,6 +191,12 @@ end
     @test !isnan(result.solution)
     @test result.data[:n_total] > 0
     @test result.data[:levels] >= 3  # levels_min + 1
+
+    traced = integrate(
+        CubMLMC(f; abs_tol=0.5, n_init=256, levels_min=2, levels_max=6, trace_iterations=true),
+    )
+    @test haskey(traced.data, :iteration_log)
+    @test length(traced.data[:iteration_log]) >= 1
 end
 
 @testset "CubMLMCCont" begin
@@ -150,6 +208,38 @@ end
     @test result.solution isa Float64
     @test !isnan(result.solution)
     @test result.data[:n_total] > 0
+
+    traced = integrate(
+        CubMLMCCont(
+            f;
+            abs_tol=0.5,
+            n_init=256,
+            levels_min=2,
+            levels_max=6,
+            n_tols=5,
+            trace_iterations=true,
+        ),
+    )
+    @test haskey(traced.data, :iteration_log)
+    @test length(traced.data[:iteration_log]) >= 1
+end
+
+@testset "CubMLQMC" begin
+    dd = Lattice(32; replications=8)
+    f = TestMLIntegrand(dd; d_coarsest=1, volatility=0.5, start_price=30.0, strike_price=35.0)
+    sc = CubMLQMC(f; abs_tol=0.5, n_init=64, levels_min=2, levels_max=6)
+    result = integrate(sc)
+
+    @test result.solution isa Float64
+    @test !isnan(result.solution)
+    @test result.data[:n_total] > 0
+    @test result.data[:replications] == 8
+
+    traced = integrate(
+        CubMLQMC(f; abs_tol=0.5, n_init=64, levels_min=2, levels_max=6, trace_iterations=true),
+    )
+    @test haskey(traced.data, :iteration_log)
+    @test length(traced.data[:iteration_log]) >= 1
 end
 
 @testset "CubMLQMCCont" begin
@@ -162,4 +252,18 @@ end
     @test !isnan(result.solution)
     @test result.data[:n_total] > 0
     @test result.data[:replications] == 8
+
+    traced = integrate(
+        CubMLQMCCont(
+            f;
+            abs_tol=0.5,
+            n_init=64,
+            levels_min=2,
+            levels_max=6,
+            n_tols=5,
+            trace_iterations=true,
+        ),
+    )
+    @test haskey(traced.data, :iteration_log)
+    @test length(traced.data[:iteration_log]) >= 1
 end
diff --git a/test/test_stopping_criteria.jl b/test/test_stopping_criteria.jl
index 97d49f5..47073db 100644
--- a/test/test_stopping_criteria.jl
+++ b/test/test_stopping_criteria.jl
@@ -79,6 +79,11 @@
         @test result.data[:n] >= 2^8
         @test result.data[:n_total] == result.data[:n]
         @test result.data[:n_per_rep] == result.data[:n]
+        traced = integrate(
+            CubQMCBayesLatticeG(f; abs_tol=0.1, n_init=2^8, n_max=2^12, trace_iterations=true),
+        )
+        @test haskey(traced.data, :iteration_log)
+        @test length(traced.data[:iteration_log]) >= 1
     end
 
     @testset "CubQMCBayesNetG (smoke)" begin
@@ -94,6 +99,11 @@
         @test result.data[:n] >= 2^8
         @test result.data[:n_total] == result.data[:n]
         @test result.data[:n_per_rep] == result.data[:n]
+        traced = integrate(
+            CubQMCBayesNetG(f; abs_tol=0.1, n_init=2^8, n_max=2^12, trace_iterations=true),
+        )
+        @test haskey(traced.data, :iteration_log)
+        @test length(traced.data[:iteration_log]) >= 1
     end
 
     @testset "CubQMCRepStudentT (smoke)" begin
@@ -138,6 +148,9 @@
         sc = CubMCG(f; abs_tol=0.1)
         @test sc.abs_tol == 0.1
         @test sc.kurtmax > 0
+        r = integrate(CubMCG(f; abs_tol=0.1, trace_iterations=true))
+        @test haskey(r.data, :iteration_log)
+        @test length(r.data[:iteration_log]) >= 1
     end
 
     @testset "CubMCCLTVec" begin
@@ -177,6 +190,11 @@
         r = integrate(sc)
         @test haskey(r.data, :iteration_log)
         @test length(r.data[:iteration_log]) >= 1
+
+        sc_clt = CubMCCLT(f; abs_tol=0.05, trace_iterations=true)
+        r_clt = integrate(sc_clt)
+        @test haskey(r_clt.data, :iteration_log)
+        @test length(r_clt.data[:iteration_log]) == 2
     end
 
     @testset "QMCResult show" begin
diff --git a/test/test_true_measures.jl b/test/test_true_measures.jl
index a670bf2..065d610 100644
--- a/test/test_true_measures.jl
+++ b/test/test_true_measures.jl
@@ -100,6 +100,45 @@
         @test abs(mean(xt)) < 0.15
     end
 
+    @testset "StudentT df=2 closed form" begin
+        dd = DigitalNetB2(2; randomize="none", seed=71)
+        x = gen_samples(dd, 16)
+
+        tm_std = StudentT(dd)
+        xt_std = transform(tm_std, x)
+        expected_std = Matrix{Float64}(undef, size(x))
+        @inbounds for j in axes(x, 2), i in axes(x, 1)
+            p = QMC._open_unit_interval(x[i, j])
+            expected_std[i, j] = (2.0 * p - 1.0) / sqrt(2.0 * p * (1.0 - p))
+        end
+        @test xt_std ≈ expected_std
+
+        loc = [1.0, -2.0]
+        scale = [0.5, 2.0]
+        tm_affine = StudentT(dd; df=2.0, loc=loc, scale=scale)
+        xt_affine = transform(tm_affine, x)
+        @test xt_affine ≈ expected_std .* transpose(scale) .+ transpose(loc)
+    end
+
+    @testset "StudentT df=1 closed form" begin
+        dd = DigitalNetB2(2; randomize="none", seed=72)
+        x = gen_samples(dd, 16)
+
+        tm_std = StudentT(dd; df=1.0)
+        xt_std = transform(tm_std, x)
+        expected_std = Matrix{Float64}(undef, size(x))
+        @inbounds for j in axes(x, 2), i in axes(x, 1)
+            expected_std[i, j] = tanpi(QMC._open_unit_interval(x[i, j]) - 0.5)
+        end
+        @test xt_std ≈ expected_std
+
+        loc = [0.25, -1.5]
+        scale = [1.5, 0.25]
+        tm_affine = StudentT(dd; df=1.0, loc=loc, scale=scale)
+        xt_affine = transform(tm_affine, x)
+        @test xt_affine ≈ expected_std .* transpose(scale) .+ transpose(loc)
+    end
+
     @testset "Triangular" begin
         dd = IIDStdUniform(2; seed=80)
         tm = Triangular(dd; lower=0.0, upper=1.0, mode=0.5)