diff --git a/.gitignore b/.gitignore
index 16c4bc5..19135d1 100644
--- a/.gitignore
+++ b/.gitignore
@@ -2,6 +2,7 @@
 *.jl.cov
 *.jl.*.cov
 *.jl.mem
+lcov.info
 
 # Ignore the root package lockfile and the formatter tool lockfile.
 # Keep docs/Manifest.toml trackable for reproducible documentation builds.
@@ -41,4 +42,4 @@ benchmark/results/
 
 # Developer-specific files
 sc_*
-lcov.info
+
diff --git a/CHANGELOG.md b/CHANGELOG.md
deleted file mode 100644
index c2a4b95..0000000
--- a/CHANGELOG.md
+++ /dev/null
@@ -1,82 +0,0 @@
-# Changelog
-
-All notable changes to QMC.jl will be documented in this file.
-
-The format follows [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
-and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
-
-## [0.1.0] - 2026-05-29
-
-Initial release — a full Julia port of QMCPy v2.3.
-
-### Discrete Distributions
-
-- `IIDStdUniform` — i.i.d. uniform sampling (pure Julia)
-- `Lattice` — rank-1 lattice with Kuo generating vectors, shift randomization, radical-inverse/Gray-code ordering, replications
-- `DigitalNetB2` — Sobol digital nets with LMS, DS, NUS, Owen scrambling, replications
-- `Halton` — Halton sequence with Owen scrambling
-- `Kronecker` — Kronecker/golden-ratio low-discrepancy sequence
-- `DigitalNetAnyBases` / `Faure` — generalized digital nets
-- Windowed sampling (`n_start` parameter) for all generators
-
-### True Measures
-
-- `Uniform`, `Gaussian`, `BrownianMotion`, `Lebesgue`, `GeometricBrownianMotion`
-- `StudentT`, `Triangular`, `Kumaraswamy`, `JohnsonsSU`, `BernoulliCont`
-- `AcceptanceRejection` — rejection sampling with custom PDF
-- `DistributionsWrapper` — bridge to Distributions.jl marginals
-- `MaternGP`, `UniformTriangle`, `ZeroInflatedExpUniform`
-
-### Integrands
-
-- `Keister`, `Genz`, `BoxIntegral`, `Linear0`, `Sin1D`
-- `Ishigami`, `Hartmann6D`, `Multimodal2D`, `FourBranch2D`
-- `AsianOption`, `FinancialOption` (European, Asian, lookback, digital, barrier)
-- `FinancialOptionML` — multilevel financial option
-- `SensitivityIndices` — Sobol sensitivity indices via pick-freeze
-- `BayesianLRCoeffs` — Bayesian logistic regression posterior
-- `UMBridgeWrapper` — UMBridge.jl model integration
-- `CustomFun` — user-defined integrand
-
-### Stopping Criteria
-
-- `CubMCCLT` — MC with CLT confidence interval
-- `CubMCCLTVec` — vectorized MC with CLT, resume/checkpoint support
-- `CubMCG` — MC with guaranteed error (kurtosis-based)
-- `CubQMCLatticeG` / `CubQMCNetG` — QMC with guaranteed error
-- `CubQMCBayesLatticeG` / `CubQMCBayesNetG` — Bayesian QMC
-- `CubQMCRepStudentT` — replicated QMC with Student-t CI
-- `CubMLMC` / `CubMLMCCont` — multilevel MC
-- `CubMLQMC` / `CubMLQMCCont` — multilevel QMC
-- `PFGPCI` — stub (requires AbstractGPs.jl backend)
-- `rel_tol` support across all criteria
-- Resume/checkpoint support for iterative criteria
-- `IterationLog` diagnostics framework
-
-### Kernels
-
-- `KernelShiftInvar`, `KernelDigShiftInvar`
-- `KernelMatern12`, `KernelMatern32`, `KernelMatern52`, `KernelGaussian`
-- `KernelMultiTask`, `SumKernel`, `ProductKernel`
-
-### Utilities
-
-- Fast Walsh-Hadamard transform (`fwht!`, `ifwht!`)
-- BRO-FFT transform
-- Periodization transforms (Baker, C1, C2, C3, C2SIN)
-- Bernoulli polynomials
-- LatNetBuilder linker
-- `mlmc_test` multilevel diagnostic utility
-- GPU backend stub
-
-### Demos
-
-- 26 Jupyter notebooks covering: quickstart, sampling, lattices, digital nets, Halton, Kronecker, financial options, multilevel QMC, sensitivity analysis, Bayesian optimization, ray tracing, UMBridge, and more
-
-### Infrastructure
-
-- Documenter.jl documentation with full API reference
-- CI/CD: fast Linux CI, cross-platform full sweep, nightly compatibility, docs build
-- JuliaFormatter configuration
-- Benchmark suite with BenchmarkTools.jl
-- Concurrency cancellation on all CI workflows
diff --git a/Makefile b/Makefile
index c76e49e..3e8c929 100644
--- a/Makefile
+++ b/Makefile
@@ -1,4 +1,8 @@
-.PHONY: test coverage doc format format-check lint clean bench bench-compare bench-compare-py bench-compare-py-label bench-all-label bench-compare-labels check-qmcpy-python
+.PHONY: test coverage doc format format-check lint clean bench bench-compare bench-compare-py bench-compare-py-label bench-all bench-compare-labels bench-coverage bench-compare-coverage bench-compare-py-coverage bench-all-coverage local-ci check-qmcpy-python
+
+# ============================================================================
+# Configuration and helpers
+# ============================================================================
 
 FORMATTER_PROJECT=devtools/formatter
 DOC_DEPOT ?= $(if $(TMPDIR),$(TMPDIR),/tmp/)qmcju-doc-depot
@@ -10,6 +14,18 @@ QMCPY_PYTHON_AUTO := $(shell \
 		fi; \
 	done)
 PYTHON ?= $(if $(QMCPY_PYTHON_AUTO),$(QMCPY_PYTHON_AUTO),python)
+BENCH_COVERAGE ?= 0
+JULIA_BENCH_COVERAGE_FLAG := $(if $(filter 1,$(BENCH_COVERAGE)),--code-coverage=user,)
+
+define RUN_TIMED
+	@start=$$(date +%s); \
+	$(1); \
+	exit_code=$$?; \
+	end=$$(date +%s); \
+	elapsed=$$((end - start)); \
+	printf '\n[%s] total runtime: %ss\n' "$(2)" "$$elapsed"; \
+	exit $$exit_code
+endef
 
 # Allow `make bench gaus` / `make bench-compare gaus` / `make bench-compare-py gaus`
 # as shorthand for `LABEL=gaus`. `make` treats `gaus` as an extra goal, so consume
@@ -33,15 +49,33 @@ ifneq ($(filter bench-compare-labels,$(firstword $(MAKECMDGOALS))),)
     LABEL_B ?= $(word 2,$(EXTRA_LABEL_COMPARE_GOALS))
     OUT_LABEL ?= $(word 3,$(EXTRA_LABEL_COMPARE_GOALS))
     .PHONY: $(EXTRA_LABEL_COMPARE_GOALS)
-    $(EXTRA_LABEL_COMPARE_GOALS):
+	    $(EXTRA_LABEL_COMPARE_GOALS):
 	@:
   endif
 endif
 
+# ============================================================================
+# Project maintenance and setup
+# ============================================================================
+
 # Update packages and resolve dependencies
 update:
 	julia --project=. -e 'using Pkg; Pkg.update(); Pkg.resolve'
 
+# Instantiate project dependencies (includes Plots and all other deps). Download what Manifest.toml says.
+setup:
+	julia --project=. -e 'using Pkg; Pkg.instantiate()'
+
+# Clean build artifacts
+clean:
+	rm -rf docs/build
+	rm -rf *.jl.cov *.jl.*.cov *.jl.mem lcov.info
+	find src benchmark test -name '*.cov' -delete
+
+# ============================================================================
+# Testing and coverage
+# ============================================================================
+
 # Run all tests
 test: 
 	julia --project=. -e 'using Pkg; Pkg.instantiate(); Pkg.test()'
@@ -52,16 +86,23 @@ coverage:
 	rm -f lcov.info
 	julia --project=. -e 'using Pkg; Pkg.instantiate(); Pkg.test(coverage=true)'
 	julia --project=. devtools/process_coverage.jl
+	find src test -name '*.cov' -delete
 
 # Run specific test file
 test-%:
 	julia --project=. -e 'include("test/$*.jl")'
 
+# ============================================================================
+# Documentation
+# ============================================================================
+
 # Build documentation
 doc:
-	rm -rf docs/build
-	JULIA_DEPOT_PATH="$(DOC_DEPOT):$(HOME)/.julia" julia --project=docs -e 'using Pkg; Pkg.instantiate(); Pkg.resolve()'
-	JULIA_DEPOT_PATH="$(DOC_DEPOT):$(HOME)/.julia" julia --project=docs docs/make.jl
+	$(call RUN_TIMED,rm -rf docs/build && JULIA_DEPOT_PATH="$(DOC_DEPOT):$(HOME)/.julia" julia --project=docs -e 'using Pkg; Pkg.instantiate(); Pkg.resolve()' && JULIA_DEPOT_PATH="$(DOC_DEPOT):$(HOME)/.julia" julia --project=docs docs/make.jl,doc)
+
+# ============================================================================
+# Formatting
+# ============================================================================
 
 # Format code with JuliaFormatter (uses the repo .JuliaFormatter.toml for all paths)
 format:
@@ -71,14 +112,9 @@ format:
 format-check:
 	julia --project=$(FORMATTER_PROJECT) -e 'using Pkg; Pkg.instantiate(); using JuliaFormatter; @assert format(["src/", "test/", "benchmark/"], overwrite=false)'
 
-# Clean build artifacts
-clean:
-	rm -rf docs/build
-	rm -rf *.jl.cov *.jl.*.cov *.jl.mem lcov.info
-
-# Instantiate project dependencies (includes Plots and all other deps). Download what Manifest.toml says.
-setup:
-	julia --project=. -e 'using Pkg; Pkg.instantiate()'
+# ============================================================================
+# Smoke checks and notebooks
+# ============================================================================
 
 # Run a quick smoke test
 smoke:
@@ -92,11 +128,15 @@ notebook:
 notebook-%:
 	julia --project=. test/run_notebooks.jl $*
 
+# ============================================================================
+# Benchmarking
+# ============================================================================
+
 # Run the benchmark suite (uses its own environment in benchmark/, set up on first run)
 # Override output label with: make bench LABEL=gaus or make bench gaus
 # Saves results to benchmark/results/latest.json (default) or <LABEL>.json
 bench:
-	julia benchmark/runbenchmarks.jl $(LABEL)
+	$(call RUN_TIMED,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 || \
@@ -112,7 +152,7 @@ check-qmcpy-python:
 REV ?= HEAD
 LABEL ?=
 bench-compare:
-	julia benchmark/compare.jl $(REV) $(LABEL)
+	$(call RUN_TIMED,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.
@@ -124,24 +164,46 @@ bench-compare:
 JL_LABEL ?= $(if $(LABEL),$(LABEL),latest)
 PY_LABEL ?= $(JL_LABEL)
 bench-compare-py: bench check-qmcpy-python
-	$(PYTHON) benchmark/benchmark_qmcpy.py $(PY_LABEL)
-	julia benchmark/compare_py.jl $(JL_LABEL) $(PY_LABEL) $(LABEL)
+	$(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)
 
 # Explicit labeled Julia-vs-QMCPy comparison flow.
 # Usage: make bench-compare-py-label LABEL=base
 bench-compare-py-label: check-qmcpy-python
-	julia benchmark/runbenchmarks.jl $(LABEL)
-	$(PYTHON) benchmark/benchmark_qmcpy.py $(LABEL)
-	julia benchmark/compare_py.jl $(LABEL) $(LABEL) $(LABEL)
+	$(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)
 
 # Run the labeled Julia-only comparison and Julia-vs-QMCPy comparison in one task.
-# Usage: make bench-all-label LABEL=base
-bench-all-label:
-	$(MAKE) bench-compare LABEL=$(LABEL)
-	$(MAKE) bench-compare-py LABEL=$(LABEL)
+# 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)
+
+# 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)
+
+# 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)
+
+# 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)
+
+# 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)
 
 # 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]
 #    or: make bench-compare-labels a b [report]
 bench-compare-labels:
 	julia benchmark/compare_labels.jl $(LABEL_A) $(LABEL_B) $(OUT_LABEL)
+
+# ============================================================================
+# Combination of above targets
+# ============================================================================
+
+# Run the local 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
diff --git a/Project.toml b/Project.toml
index 8233f2e..543317e 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,6 +7,7 @@ authors = ["Sou-Cheng T. Choi", "Fred J. Hickernell", "QMC Contributors"]
 [deps]
 NBInclude = "0db19996-df87-5ea3-a455-e3a50d440464"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -18,6 +19,7 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 [compat]
 Aqua = "0.8"
 Distributions = "0.25"
+Downloads = "1"
 FFTW = "1"
 Libdl = "1.10"
 LinearAlgebra = "1.10"
diff --git a/README.md b/README.md
index 37a9c3c..8a38328 100644
--- a/README.md
+++ b/README.md
@@ -32,7 +32,7 @@ QMC.jl provides four building blocks that snap together:
 | **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` |
 | **Integrand** | `Keister`, `Genz`, `AsianOption`, `FinancialOption`, `FinancialOptionML`, `BoxIntegral`, `Linear0`, `Sin1D`, `Ishigami`, `Hartmann6D`, `Multimodal2D`, `FourBranch2D`, `SensitivityIndices`, `BayesianLRCoeffs`, `UMBridgeWrapper`, `CustomFun` |
-| **Stopping Criterion** | `CubMCCLT`, `CubMCCLTVec`, `CubMCG`, `CubQMCLatticeG`, `CubQMCNetG`, `CubQMCBayesLatticeG`, `CubQMCBayesNetG`, `CubQMCRepStudentT`, `CubMLMC`, `CubMLMCCont`, `CubMLQMC`, `CubMLQMCCont`, `PFGPCI` |
+| **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` |
 
 Additional utilities include periodization transforms, iteration diagnostics (`IterationLog`), resume/checkpoint support, Walsh-Hadamard and BRO-FFT transforms, and a LatNetBuilder linker.
diff --git a/benchmark/README.md b/benchmark/README.md
index 8e0be7b..fe9109b 100644
--- a/benchmark/README.md
+++ b/benchmark/README.md
@@ -1,11 +1,8 @@
 # Benchmarking QMC.jl
 
-This directory contains the standalone benchmark suite, the git-revision
-comparison tooling, and the Julia-vs-QMCPy comparison tooling.
+This directory contains the standalone benchmark suite, the git-revision comparison tooling, and the Julia-vs-QMCPy comparison tooling.
 
-The benchmark scripts use the local `benchmark/Project.toml` environment, so
-benchmark-only dependencies do not need to be added to the main package
-environment.
+The benchmark scripts use the local `benchmark/Project.toml` environment, so benchmark-only dependencies do not need to be added to the main package environment.
 
 ## Files
 
@@ -17,6 +14,21 @@ environment.
 - `compare_py.jl`: compares Julia results against QMCPy results
 - `results/`: generated benchmark result files and Markdown reports
 
+## Results Directory
+
+`benchmark/results/` 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
+
+These JSON and Markdown outputs are machine-specific benchmark artifacts and are intentionally not tracked in git.
+
 ## Quick Start
 
 From the repository root:
@@ -104,8 +116,7 @@ So:
 - ratio `> 1` means `base` is better
 - ratio `< 1` means `latest` is better
 
-If one label wins both time and memory, it wins overall. If time and memory trade
-off, the script uses a combined score to break the tie and reports that explicitly.
+If one label wins both time and memory, it wins overall. If time and memory trade off, the script uses a combined score to break the tie and reports that explicitly.
 
 ## Julia vs QMCPy Comparison
 
@@ -140,6 +151,14 @@ benchmark/results/qmcpy_base.json
 benchmark/results/compare_python_base.md
 ```
 
+### Coverage Caveat
+
+Do not interpret coverage-enabled benchmark runs as representative performance comparisons.
+
+When the Julia benchmark flow is run with `--code-coverage=user`, every native Julia line is instrumented. That instrumentation can inhibit important loop optimizations, including bounds-check elision and SIMD, and adds per-iteration overhead to Julia-native loops. In practice this can disproportionately slow rows whose cost is dominated by native Julia iteration (`O(d*n)`), while BLAS-bound rows may remain comparatively stable.
+
+If a Julia-vs-QMCPy report shows Python improving sharply against Julia, first check whether the run was produced through a `bench-*-coverage` target or with `BENCH_COVERAGE=1`. Those runs are useful for correctness and coverage collection, not for clean performance claims.
+
 The Julia runner now also writes a sidecar memory file per label:
 
 - `<label>_memory.json`: per-row Julia `rss_delta_kib` from one warmed call
@@ -149,8 +168,7 @@ The Python benchmark JSON records two approximate memory metrics per row:
 - `tracemalloc_peak_kib`: Python-managed peak memory during one warmed call
 - `rss_delta_kib`: retained process RSS delta after one warmed call
 
-Together with Julia's `KiB` allocation metric from `BenchmarkTools`, this gives
-three distinct memory signals:
+Together with Julia's `KiB` allocation metric from `BenchmarkTools`, this gives three distinct memory signals:
 
 - Julia `alloc KiB`: bytes allocated during the call
 - Julia `rss_delta_kib`: retained process RSS delta after the call
@@ -169,30 +187,22 @@ When reading the weighted ratios in `compare_python*.md`:
 - ratio `< 1` means the Python metric is smaller than the Julia total
 - ratio `> 1` means the Python metric is larger than the Julia total
 
-If one memory ratio is `< 1` and the other is `> 1`, that is normal rather than
-contradictory. It usually means Python allocated more temporary memory during the
-call but released most of it afterward, or conversely retained more process
-memory even though its traced Python-level peak was modest. Use:
+If one memory ratio is `< 1` and the other is `> 1`, that is normal rather than contradictory. It usually means Python allocated more temporary memory during the call but released most of it afterward, or conversely retained more process memory even though its traced Python-level peak was modest. Use:
 
 - time ratio as the cleanest cross-language comparison
 - `tracemalloc` as a temporary-allocation signal
 - `Julia RSS delta` / `Python RSS delta` as retained-footprint signals
 
-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`.
+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 Python harness also mirrors most of the newer Julia-only benchmark rows:
 
 - large-`d` `Gaussian(diag)` and `Gaussian(dense)` transforms
 - `StudentT` and `JohnsonsSU` transforms
 - `BoxIntegral` and `Linear0` evaluate rows, including large-`d` variants
-- `Genz(gaussian_peak)` and `Genz(continuous)` via direct NumPy formulas with the
-  same default parameters Julia uses, since QMCPy exposes only oscillatory and
-  corner-peak Genz variants
+- `Genz(gaussian_peak)` and `Genz(continuous)` via direct NumPy formulas with the same default parameters Julia uses, since QMCPy exposes only oscillatory and corner-peak Genz variants
 
-If `compare_python*.md` still shows `n/a` rows after regenerating `qmcpy_<label>.json`,
-those rows do not currently have a meaningful Python counterpart in the harness.
+If `compare_python*.md` still shows `n/a` rows after regenerating `qmcpy_<label>.json`, those rows do not currently have a meaningful Python counterpart in the harness.
 
 If the Python interpreter should be overridden:
 
@@ -200,10 +210,7 @@ If the Python interpreter should be overridden:
 make bench-compare-py base PYTHON=python3
 ```
 
-By default, the Makefile now tries to auto-detect a Python that can
-`import qmcpy`, checking common candidates such as `python`, `python3`, and
-common Miniconda locations. If that guess is wrong on your machine, override it
-explicitly with `PYTHON=/path/to/python`.
+By default, the Makefile now tries to auto-detect a Python that can `import qmcpy`, checking common candidates such as `python`, `python3`, and common Miniconda locations. If that guess is wrong on your machine, override it explicitly with `PYTHON=/path/to/python`.
 
 If the Julia and QMCPy input labels should differ:
 
@@ -235,11 +242,10 @@ julia benchmark/compare_py.jl base base base
 
 ## Combined Labeled Comparison Target
 
-To run both the labeled Julia-only comparison and the labeled Julia-vs-QMCPy
-comparison in one task:
+To run both the labeled Julia-only comparison and the labeled Julia-vs-QMCPy comparison in one task:
 
 ```bash
-make bench-all-label LABEL=base
+make bench-all LABEL=base
 ```
 
 This runs:
@@ -284,11 +290,7 @@ Common generated files:
 
 ## Notes
 
-- `bench-compare` uses `PkgBenchmark` and may benchmark git revisions in a
-  temporary worktree.
+- `bench-compare` uses `PkgBenchmark` and may benchmark git revisions in a temporary worktree.
 - `bench-compare-py` requires a Python environment where `qmcpy` is installed.
-- Python-vs-Julia benchmark reports now include Python `tracemalloc` peak and, when
-  available, retained RSS delta. Those memory metrics are approximate and should be
-  read as supporting evidence, not as exact equivalents of Julia allocation `KiB`.
-- Current benchmark reports are written in a sanitized form and avoid embedding
-  host-specific system details or absolute local paths.
+- Python-vs-Julia benchmark reports now include Python `tracemalloc` peak and, when available, retained RSS delta. Those memory metrics are approximate and should be read as supporting evidence, not as exact equivalents of Julia allocation `KiB`.
+- Current benchmark reports are written in a sanitized form and avoid embedding host-specific system details or absolute local paths.
diff --git a/benchmark/benchmark_qmcpy.py b/benchmark/benchmark_qmcpy.py
index 8679bd9..aed43da 100644
--- a/benchmark/benchmark_qmcpy.py
+++ b/benchmark/benchmark_qmcpy.py
@@ -138,7 +138,7 @@ def main():
     label = sys.argv[1] if len(sys.argv) > 1 else "latest"
     results = {"gen_samples": {}, "transform": {}, "evaluate": {}, "integrate": {}}
 
-    print("QMCPy Benchmarks (qmcpy %s)" % getattr(qp, "__version__", "?"))
+    print("\n\nQMCPy Benchmarks (qmcpy %s)" % getattr(qp, "__version__", "?"))
     print("=" * 70)
 
     # 1. Discrete distribution sampling --------------------------------------
@@ -289,6 +289,9 @@ def _scalar(x):
 
     # (name, make_sc, warmup). make_sc builds a fresh stopping criterion; it is
     # reused both for timing and for the one-shot accuracy measurement below.
+    # Julia's default `CubQMCNetG` now matches QMCPy's single-net `CubQMCNetG`
+    # directly. The Julia-only replicated companion `CubQMCNetGRep` is omitted
+    # here and will therefore show `n/a` in cross-language reports.
     integrate_cases = [
         ("CubMCCLT Keister",
          lambda: qp.CubMCCLT(qp.Keister(qp.IIDStdUniform(3, seed=SEED)), abs_tol=0.01),
@@ -311,6 +314,24 @@ def _scalar(x):
                                 volatility=0.2, start_price=100, strike_price=100,
                                 interest_rate=0.05, t_final=1), abs_tol=0.5),
          False),
+        ("CubMCCLT EuropeanOption",
+         lambda: qp.CubMCCLT(
+             qp.FinancialOption(qp.IIDStdUniform(50, seed=SEED), option="EUROPEAN",
+                                volatility=0.2, start_price=100, strike_price=100,
+                                interest_rate=0.05, t_final=1), abs_tol=0.5),
+         False),
+        ("CubQMCLatticeG EuropeanOption",
+         lambda: qp.CubQMCLatticeG(
+             qp.FinancialOption(qp.Lattice(50, seed=SEED), option="EUROPEAN",
+                                volatility=0.2, start_price=100, strike_price=100,
+                                interest_rate=0.05, t_final=1), abs_tol=0.5),
+         False),
+        ("CubQMCNetG AsianOption",
+         lambda: qp.CubQMCNetG(
+             qp.FinancialOption(qp.DigitalNetB2(50, seed=SEED), option="ASIAN",
+                                volatility=0.2, start_price=100, strike_price=100,
+                                interest_rate=0.05, t_final=1), abs_tol=0.5),
+         False),
     ]
 
     for name, make_sc, warm in integrate_cases:
diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl
index 285b52e..73fae96 100644
--- a/benchmark/benchmarks.jl
+++ b/benchmark/benchmarks.jl
@@ -190,12 +190,19 @@ function _int_cubqmclatticeg_keister()
 end
 
 function _int_cubqmcnetg_keister()
-    dd = DigitalNetB2(3; seed=42, randomize="LMS_DS")
+    dd = DigitalNetB2(3; seed=42, randomize="LMS_DS", graycode=false)
     tm = Gaussian(dd; covariance=0.5)   # Keister requires N(0, I/2); matches QMCPy
     f = Keister(tm)
     CubQMCNetG(f; abs_tol=0.01)
 end
 
+function _int_cubqmcnetgrep_keister()
+    dd = DigitalNetB2(3; seed=42, randomize="LMS_DS")
+    tm = Gaussian(dd; covariance=0.5)   # Keister requires N(0, I/2); matches QMCPy
+    f = Keister(tm)
+    CubQMCNetGRep(f; abs_tol=0.01)
+end
+
 function _int_cubmcclt_asian()
     dd = IIDStdUniform(50; seed=42)
     tm = GeometricBrownianMotion(
@@ -210,6 +217,19 @@ function _int_cubmcclt_asian()
 end
 
 function _int_cubqmcnetg_european()
+    dd = DigitalNetB2(50; seed=42, randomize="LMS_DS", graycode=false)
+    tm = GeometricBrownianMotion(
+        dd;
+        volatility=0.2,
+        start_price=100.0,
+        interest_rate=0.05,
+        t_final=1.0,
+    )
+    f = FinancialOption(tm; option_type=:european, strike_price=100.0)
+    CubQMCNetG(f; abs_tol=0.5)
+end
+
+function _int_cubqmcnetgrep_european()
     dd = DigitalNetB2(50; seed=42, randomize="LMS_DS")
     tm = GeometricBrownianMotion(
         dd;
@@ -219,15 +239,73 @@ function _int_cubqmcnetg_european()
         t_final=1.0,
     )
     f = FinancialOption(tm; option_type=:european, strike_price=100.0)
+    CubQMCNetGRep(f; abs_tol=0.5)
+end
+
+function _int_cubmcclt_european()
+    dd = IIDStdUniform(50; seed=42)
+    tm = GeometricBrownianMotion(
+        dd;
+        volatility=0.2,
+        start_price=100.0,
+        interest_rate=0.05,
+        t_final=1.0,
+    )
+    f = FinancialOption(tm; option_type=:european, strike_price=100.0)
+    CubMCCLT(f; abs_tol=0.5)
+end
+
+function _int_cubqmclatticeg_european()
+    dd = Lattice(50; seed=42, randomize=true)
+    tm = GeometricBrownianMotion(
+        dd;
+        volatility=0.2,
+        start_price=100.0,
+        interest_rate=0.05,
+        t_final=1.0,
+    )
+    f = FinancialOption(tm; option_type=:european, strike_price=100.0)
+    CubQMCLatticeG(f; abs_tol=0.5)
+end
+
+function _int_cubqmcnetg_asian()
+    dd = DigitalNetB2(50; seed=42, randomize="LMS_DS", graycode=false)
+    tm = GeometricBrownianMotion(
+        dd;
+        volatility=0.2,
+        start_price=100.0,
+        interest_rate=0.05,
+        t_final=1.0,
+    )
+    f = FinancialOption(tm; option_type=:asian, strike_price=100.0)
     CubQMCNetG(f; abs_tol=0.5)
 end
 
+function _int_cubqmcnetgrep_asian()
+    dd = DigitalNetB2(50; seed=42, randomize="LMS_DS")
+    tm = GeometricBrownianMotion(
+        dd;
+        volatility=0.2,
+        start_price=100.0,
+        interest_rate=0.05,
+        t_final=1.0,
+    )
+    f = FinancialOption(tm; option_type=:asian, strike_price=100.0)
+    CubQMCNetGRep(f; abs_tol=0.5)
+end
+
 const INTEGRATE_CASES = Pair{String, Function}[
     "CubMCCLT Keister" => _int_cubmcclt_keister,
     "CubQMCLatticeG Keister" => _int_cubqmclatticeg_keister,
     "CubQMCNetG Keister" => _int_cubqmcnetg_keister,
+    "CubQMCNetGRep Keister" => _int_cubqmcnetgrep_keister,
     "CubMCCLT AsianOption" => _int_cubmcclt_asian,
     "CubQMCNetG EuropeanOption" => _int_cubqmcnetg_european,
+    "CubQMCNetGRep EuropeanOption" => _int_cubqmcnetgrep_european,
+    "CubMCCLT EuropeanOption" => _int_cubmcclt_european,
+    "CubQMCLatticeG EuropeanOption" => _int_cubqmclatticeg_european,
+    "CubQMCNetG AsianOption" => _int_cubqmcnetg_asian,
+    "CubQMCNetGRep AsianOption" => _int_cubqmcnetgrep_asian,
 ]
 
 SUITE["integrate"] = BenchmarkGroup()
diff --git a/benchmark/compare.jl b/benchmark/compare.jl
index 1111ed4..39ad742 100644
--- a/benchmark/compare.jl
+++ b/benchmark/compare.jl
@@ -28,6 +28,8 @@ let deps = keys(Pkg.project().dependencies)
     "BenchmarkTools" in deps || Pkg.add("BenchmarkTools")
     "PkgBenchmark" in deps || Pkg.add("PkgBenchmark")
 end
+# Keep the benchmark manifest in sync when the local path package changes deps.
+Pkg.resolve()
 Pkg.instantiate()
 
 using PkgBenchmark
@@ -36,6 +38,7 @@ using Printf
 
 const PKG = dirname(@__DIR__)
 const RESDIR = joinpath(@__DIR__, "results")
+const BENCH_COVERAGE = get(ENV, "BENCH_COVERAGE", "0") == "1"
 mkpath(RESDIR)
 
 comparison_outfile(label::AbstractString) =
@@ -86,6 +89,22 @@ function summary_metrics(rows)
     )
 end
 
+coverage_juliacmd() =
+    BENCH_COVERAGE ? `$(Base.julia_cmd()) --code-coverage=user` : Base.julia_cmd()
+
+function copy_coverage_files(src::AbstractString, dest::AbstractString)
+    for (dir, _, files) in walkdir(src)
+        for file in files
+            endswith(file, ".cov") || continue
+            srcfile = joinpath(dir, file)
+            rel = relpath(srcfile, src)
+            destfile = joinpath(dest, rel)
+            mkpath(dirname(destfile))
+            cp(srcfile, destfile; force=true)
+        end
+    end
+end
+
 "Write a sanitized single-run benchmark summary without host or path metadata."
 function write_single_run_md(outfile, result; label="local")
     groupdata = result.benchmarkgroup
@@ -180,7 +199,23 @@ function write_comparison_md(
 end
 
 "Benchmark the current working tree in place (uncommitted changes included)."
-bench_worktree() = benchmarkpkg(PKG; verbose=false)
+function bench_worktree()
+    snap = snapshot_package_tree(PKG)
+    try
+        return with_benchmark_env(
+            joinpath(snap, "benchmark"),
+            snap,
+            () -> benchmarkpkg(
+                snap,
+                BenchmarkConfig(; juliacmd=coverage_juliacmd());
+                verbose=false,
+            ),
+        )
+    finally
+        BENCH_COVERAGE && copy_coverage_files(snap, PKG)
+        rm(snap; force=true, recursive=true)
+    end
+end
 
 "Run `f()` with `project_dir` active and `pkgdir` developed as `QMC`."
 function with_benchmark_env(project_dir::AbstractString, pkgdir::AbstractString, f::Function)
@@ -197,10 +232,19 @@ function with_benchmark_env(project_dir::AbstractString, pkgdir::AbstractString,
     end
 end
 
+"Copy a package tree into a temporary directory without its `.git` metadata."
+function snapshot_package_tree(src::AbstractString)
+    dest = mktempdir()
+    for entry in readdir(src; join=true)
+        basename(entry) == ".git" && continue
+        cp(entry, joinpath(dest, basename(entry)); force=true)
+    end
+    return dest
+end
+
 "Benchmark a committed `rev` in a throwaway git worktree, leaving PKG untouched.
-The current benchmark/benchmarks.jl is copied into the worktree first, so the
-revision's own (possibly old or broken) benchmarks.jl is never used — both sides
-run today's suite against their respective package source."
+The worktree is snapshotted into a git-free temporary directory before running
+PkgBenchmark so the benchmark runner never needs to restore a git sha."
 function bench_revision(rev::AbstractString)
     if !success(`git -C $PKG rev-parse --verify --quiet $rev`)
         error(
@@ -213,39 +257,30 @@ function bench_revision(rev::AbstractString)
     wt = joinpath(parent, "wt")   # must not pre-exist; `git worktree add` creates it
     run(`git -C $PKG worktree add --quiet --detach $wt $rev`)
     try
-        # Run the CURRENT benchmark suite against the revision's package source, so
-        # an old/broken benchmarks.jl committed at `rev` doesn't break the run and
-        # both sides measure the same suite.
+        # Run the CURRENT benchmark suite against the revision's package source,
+        # so an old/broken benchmarks.jl committed at `rev` doesn't break the run
+        # and both sides measure the same suite.
         cp(
             joinpath(PKG, "benchmark", "benchmarks.jl"),
             joinpath(wt, "benchmark", "benchmarks.jl");
             force=true,
         )
-        # Commit the synced suite inside the throwaway worktree so its working tree
-        # is clean. PkgBenchmark checks out the commit it benchmarks and then tries
-        # to restore the original sha; on a *dirty* detached worktree that restore
-        # fails and emits "Failed to return back to original sha …". A clean tree
-        # avoids the stash/checkout dance. This only moves the worktree's detached
-        # HEAD (the main repo and `rev` are untouched), and the worktree is removed
-        # below regardless.
+
+        snap = snapshot_package_tree(wt)
         try
-            run(`git -C $wt add -A`)
-            run(
-                pipeline(
-                    `git -C $wt -c user.email=bench@localhost -c user.name=bench commit -q --allow-empty -m bench-sync-suite`;
-                    stdout=devnull,
-                    stderr=devnull,
+            return with_benchmark_env(
+                joinpath(snap, "benchmark"),
+                snap,
+                () -> benchmarkpkg(
+                    snap,
+                    BenchmarkConfig(; juliacmd=coverage_juliacmd());
+                    verbose=false,
                 ),
             )
-        catch
-            # If committing fails (e.g. git identity unavailable), benchmarking
-            # still works; PkgBenchmark may just re-emit its restore warning.
+        finally
+            BENCH_COVERAGE && copy_coverage_files(snap, PKG)
+            rm(snap; force=true, recursive=true)
         end
-        return with_benchmark_env(
-            joinpath(wt, "benchmark"),
-            wt,
-            () -> benchmarkpkg(wt; verbose=false),
-        )
     finally
         try
             run(`git -C $PKG worktree remove --force $wt`)
diff --git a/benchmark/compare_labels.jl b/benchmark/compare_labels.jl
index 755ee19..8b17d22 100644
--- a/benchmark/compare_labels.jl
+++ b/benchmark/compare_labels.jl
@@ -24,6 +24,8 @@ Pkg.activate(@__DIR__)
 let deps = keys(Pkg.project().dependencies)
     "BenchmarkTools" in deps || Pkg.add("BenchmarkTools")
 end
+# Keep the benchmark manifest in sync when the local path package changes deps.
+Pkg.resolve()
 Pkg.instantiate()
 
 using BenchmarkTools: load, median
diff --git a/benchmark/compare_py.jl b/benchmark/compare_py.jl
index 2afefaa..cf03d03 100644
--- a/benchmark/compare_py.jl
+++ b/benchmark/compare_py.jl
@@ -29,6 +29,8 @@ let deps = keys(Pkg.project().dependencies)
     "BenchmarkTools" in deps || Pkg.add("BenchmarkTools")
     "JSON3" in deps || Pkg.add("JSON3")
 end
+# Keep the benchmark manifest in sync when the local path package changes deps.
+Pkg.resolve()
 Pkg.instantiate()
 
 using BenchmarkTools
diff --git a/benchmark/runbenchmarks.jl b/benchmark/runbenchmarks.jl
index 753af8d..572e370 100644
--- a/benchmark/runbenchmarks.jl
+++ b/benchmark/runbenchmarks.jl
@@ -19,6 +19,8 @@ let deps = keys(Pkg.project().dependencies)
     "QMC" in deps || Pkg.develop(; path=dirname(@__DIR__))
     "BenchmarkTools" in deps || Pkg.add("BenchmarkTools")
 end
+# Keep the benchmark manifest in sync when the local path package changes deps.
+Pkg.resolve()
 Pkg.instantiate()
 
 using BenchmarkTools
@@ -127,14 +129,14 @@ end
 """
     collect_integrate_solutions() -> Dict
 
-Run each integrate accuracy case once (outside the timing loop) and record its
+Run each integrate case once (outside the timing loop) and record its
 solution value and the tolerances it was configured with. Used by the
 Julia-vs-Python accuracy check in `compare_py.jl` and by the standalone exact-
-value smoke summary. Relies on `INTEGRATE_ACCURACY_CASES` from benchmarks.jl.
+value smoke summary. Relies on `INTEGRATE_CASES` from benchmarks.jl.
 """
 function collect_integrate_solutions()
     out = Dict{String, Any}()
-    for (name, make_sc) in INTEGRATE_ACCURACY_CASES
+    for (name, make_sc) in INTEGRATE_CASES
         try
             sc = make_sc()
             res = integrate(sc)
diff --git a/cite_qmcju.bib b/cite_qmcju.bib
deleted file mode 100644
index 27d6deb..0000000
--- a/cite_qmcju.bib
+++ /dev/null
@@ -1,7 +0,0 @@
-@misc{qmcju2026,
-  author = {Sou-Cheng T. Choi and Fred J. Hickernell and Aleksei G. Sorokin and contributors},
-  title  = {{QMC.jl}: A Quasi-{M}onte {C}arlo {J}ulia Library},
-  year   = {2026},
-  note   = {Translated from the {QMC} {P}ython version 2.3 package},
-  url    = {https://github.com/QMCSoftware/QMC.jl}
-}
diff --git a/demos/asian_option_mlqmc.ipynb b/demos/asian_option_mlqmc.ipynb
index 038bd46..7754830 100644
--- a/demos/asian_option_mlqmc.ipynb
+++ b/demos/asian_option_mlqmc.ipynb
@@ -36,7 +36,7 @@
     "let\n",
     "    for level in 0:4\n",
     "        d_steps = 2 * 2^level\n",
-    "        dd_level = DigitalNetB2(d_steps; seed=7)\n",
+    "        dd_level = DigitalNetB2(d_steps; seed=7, randomize=\"LMS_DS\", graycode=false)\n",
     "        tm_level = GeometricBrownianMotion(dd_level;\n",
     "            volatility=0.2, start_price=100.0,\n",
     "            interest_rate=0.05, t_final=1.0)\n",
diff --git a/demos/digital_net_b2.ipynb b/demos/digital_net_b2.ipynb
index c26eac0..8ad1c4d 100644
--- a/demos/digital_net_b2.ipynb
+++ b/demos/digital_net_b2.ipynb
@@ -15,7 +15,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": null,
    "id": "9eb4cc7c",
    "metadata": {
     "execution": {
@@ -42,7 +42,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": null,
    "id": "bcf2bf7b",
    "metadata": {
     "execution": {
@@ -98,7 +98,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": null,
    "id": "8ab53c6f",
    "metadata": {
     "execution": {
@@ -147,7 +147,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": null,
    "id": "9369097a",
    "metadata": {
     "execution": {
@@ -216,7 +216,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": null,
    "id": "6fbb3613",
    "metadata": {
     "execution": {
@@ -267,7 +267,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": null,
    "id": "314f54f1",
    "metadata": {
     "execution": {
@@ -329,7 +329,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": null,
    "id": "119e54b7",
    "metadata": {
     "execution": {
@@ -392,7 +392,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 8,
+   "execution_count": null,
    "id": "b8968ec7",
    "metadata": {
     "execution": {
@@ -444,7 +444,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 9,
+   "execution_count": null,
    "id": "79c3db13",
    "metadata": {
     "execution": {
@@ -476,14 +476,14 @@
     "println(\"Integration: ∫₀¹ ∫₀¹ (x₁ + x₂) dx = 1.0\")\n",
     "println(\"=\"^60)\n",
     "\n",
-    "dd = DigitalNetB2(2; randomize=\"LMS_DS\", seed=7)\n",
+    "dd = DigitalNetB2(2; randomize=\"LMS_DS\", graycode=false, seed=7)\n",
     "tm = Uniform(dd)\n",
     "f = CustomFun(tm, x -> sum(x, dims=2)[:])\n",
-    "sc = CubQMCNetG(f; abs_tol=0.001, n_init=2^8, n_reps=16)\n",
+    "sc = CubQMCNetG(f; abs_tol=0.001, n_init=2^8)\n",
     "result = integrate(sc)\n",
     "@printf(\"  Solution:  %.6f  (exact = 1.0)\\n\", result.solution)\n",
     "@printf(\"  Error:     %.2e\\n\", abs(result.solution - 1.0))\n",
-    "@printf(\"  Samples:   %d per replication\\n\", result.data[:n])\n",
+    "@printf(\"  Samples:   %d\\n\", result.data[:n])\n",
     "\n",
     "println()\n",
     "println(\"=\"^60)\n",
diff --git a/demos/iris.ipynb b/demos/iris.ipynb
index 46f669e..a487c06 100644
--- a/demos/iris.ipynb
+++ b/demos/iris.ipynb
@@ -203,7 +203,7 @@
    "metadata": {},
    "source": [
     "if HAVE_DECISIONTREE\n",
-    "    dd = DigitalNetB2(d_hp; seed=7)\n",
+    "    dd = DigitalNetB2(d_hp; seed=7, randomize=\"LMS_DS\", graycode=false)\n",
     "    tm = Uniform(dd; lower_bound=hp_lb, upper_bound=hp_ub)\n",
     "    cf = CustomFun(tm, x -> get_dt_accuracy(x))\n",
     "\n",
diff --git a/demos/lebesgue_integration.ipynb b/demos/lebesgue_integration.ipynb
index 4756502..1bcc06e 100644
--- a/demos/lebesgue_integration.ipynb
+++ b/demos/lebesgue_integration.ipynb
@@ -186,12 +186,12 @@
     "@printf(\"  IID MC:       %.5f  (error ≈ %.2e, n = %d)\\n\",\n",
     "        result.solution, abs(result.solution - true_value_2), result.data[:n])\n",
     "\n",
-    "dd = DigitalNetB2(3; seed=7)\n",
+    "dd = DigitalNetB2(3; seed=7, randomize=\"LMS_DS\", graycode=false)\n",
     "tm = Uniform(dd)\n",
     "f = CustomFun(tm, x -> exp.(x[:, 1] .* x[:, 2] .* x[:, 3]))\n",
-    "sc = CubQMCNetG(f; abs_tol=0.001, n_init=2^8, n_reps=16)\n",
+    "sc = CubQMCNetG(f; abs_tol=0.001, n_init=2^8)\n",
     "result = integrate(sc)\n",
-    "@printf(\"  Sobol' QMC:   %.5f  (error ≈ %.2e, n/rep = %d)\\n\",\n",
+    "@printf(\"  Sobol' QMC:   %.5f  (error ≈ %.2e, n = %d)\\n\",\n",
     "        result.solution, abs(result.solution - true_value_2), result.data[:n])\n",
     "println()"
    ]
diff --git a/demos/pricing_options.ipynb b/demos/pricing_options.ipynb
index a3c93fd..922d639 100644
--- a/demos/pricing_options.ipynb
+++ b/demos/pricing_options.ipynb
@@ -342,14 +342,14 @@
     }
    ],
    "source": [
-    "dd = DigitalNetB2(d; randomize=\"LMS_DS\", seed=7)\n",
+    "dd = DigitalNetB2(d; randomize=\"LMS_DS\", graycode=false, seed=7)\n",
     "tm = BrownianMotion(dd)\n",
     "f_as_call_dn = FinancialOption(tm;\n",
     "    option_type=:asian, call_put=:call, mean_type=:arithmetic,\n",
     "    volatility=σ, start_price=S0, strike_price=K, interest_rate=r)\n",
-    "sc = CubQMCNetG(f_as_call_dn; abs_tol=0.05, n_init=2^8, n_reps=16)\n",
+    "sc = CubQMCNetG(f_as_call_dn; abs_tol=0.05, n_init=2^8)\n",
     "result = integrate(sc)\n",
-    "@printf(\"  Asian Call (Sobol' QMC):  %.4f  (n/rep = %d)\\n\", result.solution, result.data[:n])"
+    "@printf(\"  Asian Call (Sobol' QMC):  %.4f  (n = %d)\\n\", result.solution, result.data[:n])"
    ]
   },
   {
diff --git a/demos/qmc.jl_intro.ipynb b/demos/qmc.jl_intro.ipynb
index e8941a3..627e554 100644
--- a/demos/qmc.jl_intro.ipynb
+++ b/demos/qmc.jl_intro.ipynb
@@ -555,12 +555,12 @@
     }
    ],
    "source": [
-    "dd = DigitalNetB2(d; seed=7)\n",
+    "dd = DigitalNetB2(d; seed=7, randomize=\"LMS_DS\", graycode=false)\n",
     "tm = Uniform(dd)\n",
     "f = Genz(tm; kind=:continuous, a=a, u=u)\n",
-    "sc = CubQMCNetG(f; abs_tol=0.001, n_init=2^8, n_reps=16)\n",
+    "sc = CubQMCNetG(f; abs_tol=0.001, n_init=2^8)\n",
     "result = integrate(sc)\n",
-    "@printf(\"  Sobol' QMC:    %.6f (exact = %.6f, n/rep = %d)\\n\",\n",
+    "@printf(\"  Sobol' QMC:    %.6f (exact = %.6f, n = %d)\\n\",\n",
     "        result.solution, exact, result.data[:n])\n",
     "println()\n",
     "\n",
diff --git a/demos/quickstart.ipynb b/demos/quickstart.ipynb
index 82d3294..24bca62 100644
--- a/demos/quickstart.ipynb
+++ b/demos/quickstart.ipynb
@@ -189,10 +189,10 @@
     "println(\"Method 3: Digital Net (Sobol') QMC (CubQMCNetG)\")\n",
     "println(\"=\"^60)\n",
     "\n",
-    "dd = DigitalNetB2(d; randomize=\"LMS_DS\", seed=7)\n",
+    "dd = DigitalNetB2(d; randomize=\"LMS_DS\", graycode=false, seed=7)\n",
     "tm = Gaussian(dd; covariance=0.5)\n",
     "f = Keister(tm)\n",
-    "sc = CubQMCNetG(f; abs_tol=0.001, n_init=2^10, n_reps=16)\n",
+    "sc = CubQMCNetG(f; abs_tol=0.001, n_init=2^10)\n",
     "result = integrate(sc)\n",
     "\n",
     "@printf(\"  QMC solution:  %.6f\\n\", result.solution)\n",
diff --git a/demos/umbridge.ipynb b/demos/umbridge.ipynb
index 291fa54..e97ba68 100644
--- a/demos/umbridge.ipynb
+++ b/demos/umbridge.ipynb
@@ -87,7 +87,7 @@
    "metadata": {},
    "source": [
     "# Discrete distribution and true measure\n",
-    "sampler = DigitalNetB2(3; seed=7)\n",
+    "sampler = DigitalNetB2(3; seed=7, randomize=\"LMS_DS\", graycode=false)\n",
     "distribution = Uniform(sampler; lower_bound=1.0, upper_bound=1.05)\n",
     "\n",
     "println(\"Sampler: \", sampler)\n",
diff --git a/demos/vectorized_qmc_bayes.ipynb b/demos/vectorized_qmc_bayes.ipynb
index 4a8c811..7dfd77d 100644
--- a/demos/vectorized_qmc_bayes.ipynb
+++ b/demos/vectorized_qmc_bayes.ipynb
@@ -186,7 +186,7 @@
     "            CubQMCBayesLatticeG(Keister(tm_cmp); abs_tol=tol)\n",
     "        end),\n",
     "        (\"CubQMCNetG\",          () -> begin\n",
-    "            dd_cmp = DigitalNetB2(d; seed=7)\n",
+    "            dd_cmp = DigitalNetB2(d; seed=7, randomize=\"LMS_DS\", graycode=false)\n",
     "            tm_cmp = Gaussian(dd_cmp; mean=0.0, covariance=0.5)\n",
     "            CubQMCNetG(Keister(tm_cmp); abs_tol=tol)\n",
     "        end),\n",
diff --git a/devtools/process_coverage.jl b/devtools/process_coverage.jl
index 3d90a8b..ae10e35 100644
--- a/devtools/process_coverage.jl
+++ b/devtools/process_coverage.jl
@@ -1,6 +1,6 @@
 #!/usr/bin/env julia
 
-const COVERAGE_ROOTS = ("src",)
+coverage_roots() = isempty(ARGS) ? ("src",) : Tuple(ARGS)
 
 function cov_to_source_path(path::AbstractString)
     return replace(path, r"\.(?:\d+\.)?cov$" => "")
@@ -78,7 +78,7 @@ function write_lcov(path::AbstractString, merged::Dict{String, Vector{Int}})
     end
 end
 
-function print_summary(merged::Dict{String, Vector{Int}})
+function print_summary(merged::Dict{String, Vector{Int}}, roots)
     rows = Tuple{String, Int, Int, Float64}[]
     total_covered = 0
     total_exec = 0
@@ -89,7 +89,7 @@ function print_summary(merged::Dict{String, Vector{Int}})
         push!(rows, (source, covered, exec, pct))
     end
     total_pct = total_exec == 0 ? 100.0 : 100 * total_covered / total_exec
-    println("Coverage summary (src):")
+    println("Coverage summary ($(join(roots, ", "))):")
     for (source, covered, exec, pct) in rows
         pct_str = string(round(pct; digits=2), "%")
         println("  $(rpad(source, 56)) $(lpad(pct_str, 8))  ($covered/$exec)")
@@ -98,6 +98,7 @@ function print_summary(merged::Dict{String, Vector{Int}})
     println("  total: $(round(total_pct; digits=2))% ($total_covered/$total_exec executable lines)")
 end
 
-merged, _ = collect_coverage(COVERAGE_ROOTS)
+roots = coverage_roots()
+merged, _ = collect_coverage(roots)
 write_lcov("lcov.info", merged)
-print_summary(merged)
+print_summary(merged, roots)
diff --git a/docs/Manifest.toml b/docs/Manifest.toml
index a64b3f6..ee1f724 100644
--- a/docs/Manifest.toml
+++ b/docs/Manifest.toml
@@ -406,7 +406,7 @@ uuid = "43287f4e-b6f4-7ad1-bb20-aadabca52c3d"
 version = "1.4.0"
 
 [[deps.QMC]]
-deps = ["Distributions", "FFTW", "Libdl", "LinearAlgebra", "NBInclude", "Printf", "Random", "SpecialFunctions", "Statistics"]
+deps = ["Distributions", "Downloads", "FFTW", "Libdl", "LinearAlgebra", "NBInclude", "Printf", "Random", "SpecialFunctions", "Statistics"]
 path = ".."
 uuid = "b8e7c4a1-3f5d-4e9a-b2c6-1a8d9f0e7c3b"
 version = "0.1.0"
diff --git a/docs/src/api/stopping_criteria.md b/docs/src/api/stopping_criteria.md
index 45e2610..cae1416 100644
--- a/docs/src/api/stopping_criteria.md
+++ b/docs/src/api/stopping_criteria.md
@@ -14,6 +14,7 @@ CubMCCLTVec
 CubMCG
 CubQMCLatticeG
 CubQMCNetG
+CubQMCNetGRep
 CubQMCBayesLatticeG
 CubQMCBayesNetG
 CubQMCRepStudentT
diff --git a/docs/src/components.md b/docs/src/components.md
index 447901f..e1db606 100644
--- a/docs/src/components.md
+++ b/docs/src/components.md
@@ -61,7 +61,8 @@ Adaptive algorithms that determine sample size:
 
 - **`CubMCCLT`** — IID Monte Carlo with CLT confidence interval
 - **`CubQMCLatticeG`** — replicated randomized lattice rule
-- **`CubQMCNetG`** — replicated randomized digital net
+- **`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)
 
diff --git a/src/QMC.jl b/src/QMC.jl
index 2f185dd..d3331ea 100644
--- a/src/QMC.jl
+++ b/src/QMC.jl
@@ -1,6 +1,7 @@
 module QMC
 
-using LinearAlgebra, Libdl, Random, Statistics, FFTW, SpecialFunctions, Distributions, Printf
+using LinearAlgebra,
+    Libdl, Random, Statistics, FFTW, SpecialFunctions, Distributions, Printf, Downloads
 
 # Utilities
 include("util/bernoulli.jl")
@@ -77,6 +78,7 @@ include("kernel/multitask.jl")
 include("stopping_criterion/cub_mc_clt.jl")
 include("stopping_criterion/cub_qmc_lattice_g.jl")
 include("stopping_criterion/cub_qmc_net_g.jl")
+include("stopping_criterion/cub_qmc_net_g_rep.jl")
 include("stopping_criterion/cub_qmc_bayes_lattice_g.jl")
 include("stopping_criterion/cub_qmc_bayes_net_g.jl")
 include("stopping_criterion/cub_mlmc.jl")
@@ -122,7 +124,7 @@ export SumKernel, ProductKernel, kernel_eval, kernel_matrix
 export KernelMultiTask
 
 # Exports — stopping criteria
-export CubMCCLT, CubQMCLatticeG, CubQMCNetG
+export CubMCCLT, CubQMCLatticeG, CubQMCNetG, CubQMCNetGRep, CubQMCNetGSingle
 export CubQMCBayesLatticeG, CubQMCBayesNetG
 export CubMCG
 export CubMCCLTVec
diff --git a/src/abstract_types.jl b/src/abstract_types.jl
index 4c2c018..c6f5ec4 100644
--- a/src/abstract_types.jl
+++ b/src/abstract_types.jl
@@ -17,11 +17,14 @@ abstract type AbstractStoppingCriterion end
 Update the error tolerance(s) on a stopping criterion in place, mirroring
 QMCPy's `set_tolerance`. Only the keyword(s) you pass are changed; the others are
 left untouched. The single-level `Cub*` criteria expose `abs_tol`/`rel_tol`; the
-multilevel criteria expose `rmse_tol`. Throws if the criterion has no field for a
-requested tolerance. Returns `sc`.
+multilevel criteria expose `rmse_tol` or `target_tol`. Throws if the criterion
+has no field for a requested tolerance. Returns `sc`.
 
 ```julia
-sc = CubQMCNetG(Keister(Gaussian(DigitalNetB2(3); covariance=0.5)); abs_tol=0.05)
+sc = CubQMCNetG(
+    Keister(Gaussian(DigitalNetB2(3; randomize="LMS_DS", graycode=false, seed=7); covariance=0.5));
+    abs_tol=0.05,
+)
 set_tolerance!(sc; abs_tol=0.01, rel_tol=0.0)
 ```
 """
@@ -44,10 +47,14 @@ function set_tolerance!(
         sc.rel_tol = Float64(rel_tol)
     end
     if rmse_tol !== nothing
-        hasfield(typeof(sc), :rmse_tol) ||
-            throw(ArgumentError("$(typeof(sc)) has no rmse_tol field"))
         rmse_tol >= 0 || throw(ArgumentError("rmse_tol must be >= 0"))
-        sc.rmse_tol = Float64(rmse_tol)
+        if hasfield(typeof(sc), :rmse_tol)
+            sc.rmse_tol = Float64(rmse_tol)
+        elseif hasfield(typeof(sc), :target_tol)
+            sc.target_tol = Float64(rmse_tol)
+        else
+            throw(ArgumentError("$(typeof(sc)) has no rmse_tol/target_tol field"))
+        end
     end
     return sc
 end
@@ -60,6 +67,8 @@ abstract type AbstractKernel end
 const _OPEN01_LOW = eps(Float64)
 const _OPEN01_HIGH = 1.0 - eps(Float64)
 @inline _open_unit_interval(u::Float64) = clamp(u, _OPEN01_LOW, _OPEN01_HIGH)
+@inline _open_unit_interval(u::T) where {T <: AbstractFloat} = clamp(u, eps(T), one(T) - eps(T))
+@inline _open_unit_interval(u::Real) = _open_unit_interval(Float64(u))
 
 # Common interface functions
 """
diff --git a/src/discrete_distribution/kronecker.jl b/src/discrete_distribution/kronecker.jl
index 1f03015..6ff9f4e 100644
--- a/src/discrete_distribution/kronecker.jl
+++ b/src/discrete_distribution/kronecker.jl
@@ -90,9 +90,13 @@ function gen_samples(kr::Kronecker, n::Int; n_start::Int=0)
     d = kr.dimension
     x = Matrix{Float64}(undef, n, d)
 
-    @inbounds for i in 1:n
-        for j in 1:d
-            x[i, j] = mod((n_start + i) * kr.alpha[j] + kr.shift[j], 1.0)
+    # Julia matrices are column-major, so fill one column at a time rather than
+    # striding by `n` across rows in the inner loop.
+    @inbounds for j in 1:d
+        αj = kr.alpha[j]
+        shiftj = kr.shift[j]
+        for i in 1:n
+            x[i, j] = mod((n_start + i) * αj + shiftj, 1.0)
         end
     end
 
diff --git a/src/discrete_distribution/lattice.jl b/src/discrete_distribution/lattice.jl
index 6636f2a..14b92d4 100644
--- a/src/discrete_distribution/lattice.jl
+++ b/src/discrete_distribution/lattice.jl
@@ -31,7 +31,11 @@ before `using QMC`.
     `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.
+    the 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
+    downloaded on demand.
 
 # Examples
 ```julia
@@ -55,6 +59,7 @@ mutable struct Lattice{R <: AbstractRNG} <: AbstractDiscreteDistribution
 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/"
 
 function _normalize_lattice_order(order::String)
     # Normalize order aliases to canonical tokens, matching QMCPy semantics:
@@ -124,6 +129,51 @@ function _read_lattice_vector_file(path::AbstractString)
     return values, nothing
 end
 
+function _canonical_lattice_vector_filename(name::AbstractString)
+    token = replace(strip(name), '\\' => '/')
+    token = split(token, '?'; limit=2)[1]
+    token = split(token, '#'; limit=2)[1]
+    isempty(token) && throw(ArgumentError("generating_vector reference cannot be empty"))
+    parts = split(token, '/')
+    basename = parts[end]
+    isempty(basename) &&
+        throw(ArgumentError("generating_vector reference \"$name\" has no filename"))
+    return endswith(lowercase(basename), ".txt") ? basename : string(basename, ".txt")
+end
+
+function _looks_like_lddata_lattice_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, "lattice/") ||
+        occursin("/lattice/", token)
+    )
+end
+
+function _download_lattice_vector_from_lddata(filename::AbstractString)
+    url = _LATTICE_LDDATA_RAW_BASE * filename
+    path, io = mktemp()
+    close(io)
+    try
+        Downloads.download(url, path)
+        return _read_lattice_vector_file(path)
+    catch err
+        msg = sprint(showerror, err)
+        throw(
+            ArgumentError(
+                "failed to fetch LDData lattice generating vector \"$filename\" from $url: $msg",
+            ),
+        )
+    finally
+        isfile(path) && rm(path; force=true)
+    end
+end
+
 function _random_lattice_vector(dimension::Int, m::Integer, rng::AbstractRNG)
     1 < m < 27 ||
         throw(ArgumentError("integer generating_vector must satisfy 1 < M < 27, got $m"))
@@ -139,7 +189,7 @@ function _random_lattice_vector(dimension::Int, m::Integer, rng::AbstractRNG)
 end
 
 function _resolve_lattice_generating_vector(dimension::Int, generating_vector, rng::AbstractRNG)
-    if isnothing(generating_vector) || generating_vector == _DEFAULT_LATTICE_VECTOR_NAME
+    if isnothing(generating_vector)
         dimension <= _KUO_LATTICE_MAX_DIM || throw(
             ArgumentError(
                 "dimension $dimension exceeds maximum supported ($_KUO_LATTICE_MAX_DIM)",
@@ -151,14 +201,27 @@ function _resolve_lattice_generating_vector(dimension::Int, generating_vector, r
     elseif generating_vector isa Integer
         return _random_lattice_vector(dimension, generating_vector, rng)
     elseif generating_vector isa AbstractString
-        isfile(generating_vector) ||
+        if isfile(generating_vector)
+            values, n_limit = _read_lattice_vector_file(generating_vector)
+            return _coerce_lattice_vector(values, dimension), n_limit
+        end
+        _looks_like_lddata_lattice_reference(generating_vector) ||
             throw(ArgumentError("generating_vector file \"$generating_vector\" not found"))
-        values, n_limit = _read_lattice_vector_file(generating_vector)
+        filename = _canonical_lattice_vector_filename(generating_vector)
+        if filename == _DEFAULT_LATTICE_VECTOR_NAME
+            dimension <= _KUO_LATTICE_MAX_DIM || throw(
+                ArgumentError(
+                    "dimension $dimension exceeds maximum supported ($_KUO_LATTICE_MAX_DIM)",
+                ),
+            )
+            return _KUO_LATTICE_GEN_VECTOR[1:dimension], 1 << 20
+        end
+        values, n_limit = _download_lattice_vector_from_lddata(filename)
         return _coerce_lattice_vector(values, dimension), n_limit
     else
         throw(
             ArgumentError(
-                "generating_vector must be nothing, an integer, a vector of integers, or a local file path",
+                "generating_vector must be nothing, an integer, a vector of integers, a local file path, or an LDData lattice filename/URL",
             ),
         )
     end
diff --git a/src/integrand/asian_option.jl b/src/integrand/asian_option.jl
index b1ecbb5..eae6264 100644
--- a/src/integrand/asian_option.jl
+++ b/src/integrand/asian_option.jl
@@ -72,32 +72,40 @@ function evaluate(f::AsianOption, x::AbstractMatrix)
     r = f.interest_rate
     T = f._time_vector[end]
     tv = f._time_vector
-    y = Vector{Float64}(undef, n)
-
-    for i in 1:n
-        # Compute stock prices at each monitoring date
-        # x[i,:] contains Brownian motion values W(t_1), ..., W(t_d)
-        stock_prices = Vector{Float64}(undef, d)
-        for j in 1:d
-            stock_prices[j] = S0 * exp((r - 0.5 * σ^2) * tv[j] + σ * x[i, j])
+    drift = r - 0.5 * σ^2
+    avg = zeros(Float64, n)
+    # Julia matrices are column-major, so the natural `i`-outer / `j`-inner
+    # stock-price loop would stride by `n` across `x[i, j]`. Accumulating the
+    # average with `j` outer / `i` inner reads each Brownian-motion column
+    # contiguously while preserving the per-row monitoring-date order.
+    if f.mean_type == :arithmetic
+        @inbounds for j in 1:d
+            tj = tv[j]
+            @simd for i in 1:n
+                avg[i] += S0 * exp(drift * tj + σ * x[i, j])
+            end
         end
-
-        # Compute average
-        if f.mean_type == :arithmetic
-            avg_price = mean(stock_prices)
-        else  # :geometric
-            avg_price = exp(mean(log.(stock_prices)))
+        @. avg = avg / d
+    else  # geometric
+        logS0 = log(S0)
+        @inbounds for j in 1:d
+            tj = tv[j]
+            @simd for i in 1:n
+                avg[i] += logS0 + drift * tj + σ * x[i, j]
+            end
         end
+        @. avg = exp(avg / d)
+    end
 
-        # Compute payoff
-        if f.call_put == :call
-            payoff = max(avg_price - K, 0.0)
-        else  # :put
-            payoff = max(K - avg_price, 0.0)
+    discount = exp(-r * T)
+    y = Vector{Float64}(undef, n)
+    @inbounds for i in 1:n
+        payoff = if f.call_put == :call
+            max(avg[i] - K, 0.0)
+        else
+            max(K - avg[i], 0.0)
         end
-
-        # Discount
-        y[i] = exp(-r * T) * payoff
+        y[i] = discount * payoff
     end
 
     return y
diff --git a/src/integrand/bayesian_lr_coeffs.jl b/src/integrand/bayesian_lr_coeffs.jl
index d540730..7767525 100644
--- a/src/integrand/bayesian_lr_coeffs.jl
+++ b/src/integrand/bayesian_lr_coeffs.jl
@@ -65,11 +65,15 @@ where eta_i = beta_0 + X_i' * beta_{1:p}
 function _log_likelihood(f::BayesianLRCoeffs, beta::AbstractVector)
     ll = 0.0
     offset = f.dimension == f.n_features + 1 ? 1 : 0
+    # Julia stores matrices column-major, so `X * β` uses a cache-friendly
+    # column reduction / BLAS GEMV instead of striding across rows in a nested
+    # `for i; for j; X[i, j]` loop.
+    eta_vec = f.feature_array * @view(beta[(offset + 1):end])
+    if offset == 1
+        eta_vec .+= beta[1]
+    end
     @inbounds for i in 1:f.n_obs
-        eta = offset == 1 ? beta[1] : 0.0
-        for j in 1:f.n_features
-            eta += f.feature_array[i, j] * beta[j + offset]
-        end
+        eta = eta_vec[i]
         # Numerically stable log(1 + exp(eta))
         if eta > 20
             ll += f.response_vector[i] * eta - eta
diff --git a/src/integrand/box_integral.jl b/src/integrand/box_integral.jl
index 1ba738c..5de88f5 100644
--- a/src/integrand/box_integral.jl
+++ b/src/integrand/box_integral.jl
@@ -22,13 +22,20 @@ BoxIntegral(tm::AbstractTrueMeasure; s::Float64=2.0) = BoxIntegral(tm, tm.dimens
 function evaluate(f::BoxIntegral, x::AbstractMatrix)
     n, d = size(x)
     s_half = f.s / 2.0
+    # Julia matrices are column-major, so the natural row loop
+    # `for i in 1:n; for j in 1:d; x[i, j]^2 end end` makes the inner loop walk
+    # memory with stride `n`, which is cache-unfriendly at large `d`. Making `j`
+    # the outer loop reads each column contiguously and accumulates into a
+    # length-`n` vector; the per-row summation order remains unchanged.
+    r2 = zeros(Float64, n)
+    @inbounds for j in 1:d
+        @simd for i in 1:n
+            r2[i] += x[i, j]^2
+        end
+    end
     y = Vector{Float64}(undef, n)
     @inbounds for i in 1:n
-        r2 = 0.0
-        for j in 1:d
-            r2 += x[i, j]^2
-        end
-        y[i] = r2 == 0.0 && s_half < 0.0 ? 0.0 : r2^s_half
+        y[i] = r2[i] == 0.0 && s_half < 0.0 ? 0.0 : r2[i]^s_half
     end
     return y
 end
diff --git a/src/integrand/financial_option.jl b/src/integrand/financial_option.jl
index 3e05e52..6363c72 100644
--- a/src/integrand/financial_option.jl
+++ b/src/integrand/financial_option.jl
@@ -210,11 +210,128 @@ function evaluate(f::FinancialOption, x::AbstractMatrix)
         end
     end
 
+    # Fast path: Asian option with GBM transform already applied. The transform
+    # gives the full price path in x, so the payoff is the mean across the d
+    # monitoring dates. Doing it as a column reduction (j outer, i inner) reads x
+    # contiguously and skips the per-row _stock_prices/_payoff calls and the
+    # per-row log allocation; per-row accumulation order (j = 1:d) is unchanged,
+    # so results match the generic path.
+    if f.option_type == :asian && f._use_gbm_transform
+        n, d = size(x)
+        K = f.strike_price
+        avg = zeros(Float64, n)
+        if f.mean_type == :arithmetic
+            @inbounds for j in 1:d
+                @simd for i in 1:n
+                    avg[i] += x[i, j]
+                end
+            end
+            @. avg = avg / d
+        else  # geometric: exp(mean(log S))
+            @inbounds for j in 1:d
+                @simd for i in 1:n
+                    avg[i] += log(x[i, j])
+                end
+            end
+            @. avg = exp(avg / d)
+        end
+        return if f.call_put == :call
+            @. discount * max(avg - K, 0.0)
+        else
+            @. discount * max(K - avg, 0.0)
+        end
+    end
+
+    # Fast path: digital option with GBM transform already applied.
+    # Only the terminal price S(T) = x[:, end] matters.
+    if f.option_type == :digital && f._use_gbm_transform
+        S_T = @view x[:, end]
+        K = f.strike_price
+        return if f.call_put == :call
+            @. discount * ifelse(S_T > K, 1.0, 0.0)
+        else
+            @. discount * ifelse(S_T < K, 1.0, 0.0)
+        end
+    end
+
+    # Fast path: lookback option with GBM transform already applied.
+    # Reduce the path extrema in column-major order instead of rebuilding each
+    # row's stock-price vector.
+    if f.option_type == :lookback && f._use_gbm_transform
+        n, d = size(x)
+        extrema = copy(@view x[:, 1])
+        if f.call_put == :call
+            @inbounds for j in 2:d
+                @simd for i in 1:n
+                    extrema[i] = max(extrema[i], x[i, j])
+                end
+            end
+            K = f.strike_price
+            return @. discount * max(extrema - K, 0.0)
+        else
+            @inbounds for j in 2:d
+                @simd for i in 1:n
+                    extrema[i] = min(extrema[i], x[i, j])
+                end
+            end
+            K = f.strike_price
+            return @. discount * max(K - extrema, 0.0)
+        end
+    end
+
+    # Fast path: barrier option with GBM transform already applied.
+    # Track the terminal price and the path crossing state by reducing columns.
+    if f.option_type == :barrier && f._use_gbm_transform
+        n, d = size(x)
+        B = f.barrier_price
+        K = f.strike_price
+        ST = @view x[:, end]
+        crossed = falses(n)
+        if f.start_price < B
+            @inbounds for j in 1:d
+                @simd for i in 1:n
+                    crossed[i] |= x[i, j] >= B
+                end
+            end
+        else
+            @inbounds for j in 1:d
+                @simd for i in 1:n
+                    crossed[i] |= x[i, j] <= B
+                end
+            end
+        end
+        active = f.barrier_in_out == :in ? crossed : .!crossed
+        return if f.call_put == :call
+            @. discount * ifelse(active, max(ST - K, 0.0), 0.0)
+        else
+            @. discount * ifelse(active, max(K - ST, 0.0), 0.0)
+        end
+    end
+
     n = size(x, 1)
     y = Vector{Float64}(undef, n)
-    @inbounds for i in 1:n
-        S = _stock_prices(f, @view(x[i, :]))
-        y[i] = discount * _payoff(f, S)
+    if f._use_gbm_transform
+        # Prices already in x; pass row views straight to the payoff — no per-row
+        # price-vector allocation.
+        @inbounds for i in 1:n
+            y[i] = discount * _payoff(f, @view(x[i, :]))
+        end
+    else
+        # Non-GBM: build the price path into a single reused buffer instead of
+        # allocating a length-d vector per row inside _stock_prices. Also hoists
+        # the GBM/non-GBM branch out of the loop.
+        d = f.dimension
+        σ = f.volatility
+        S0 = f.start_price
+        r = f.interest_rate
+        tv = f._time_vector
+        Sbuf = Vector{Float64}(undef, d)
+        @inbounds for i in 1:n
+            for j in 1:d
+                Sbuf[j] = S0 * exp((r - 0.5 * σ^2) * tv[j] + σ * x[i, j])
+            end
+            y[i] = discount * _payoff(f, Sbuf)
+        end
     end
     return y
 end
diff --git a/src/integrand/genz.jl b/src/integrand/genz.jl
index ba5421e..10c285e 100644
--- a/src/integrand/genz.jl
+++ b/src/integrand/genz.jl
@@ -59,53 +59,79 @@ function evaluate(f::Genz, x::AbstractMatrix)
         @. y = cos(2π * u[1] + s)
 
     elseif f.kind == :product_peak
-        for i in 1:n
-            p = 1.0
-            for j in 1:d
-                p *= 1.0 / (a[j]^(-2) + (x[i, j] - u[j])^2)
+        # Julia matrices are column-major, so keeping `i` outer and `j` inner
+        # would make the inner loop stride by `n`. Swapping to `j` outer / `i`
+        # inner reads each column contiguously while preserving the per-row
+        # product order.
+        fill!(y, 1.0)
+        @inbounds for j in 1:d
+            aj_m2 = a[j]^(-2)
+            uj = u[j]
+            @simd for i in 1:n
+                y[i] *= 1.0 / (aj_m2 + (x[i, j] - uj)^2)
             end
-            y[i] = p
         end
 
     elseif f.kind == :corner_peak
-        for i in 1:n
-            s = 1.0
-            for j in 1:d
-                s += a[j] * x[i, j]
+        # Same column-major cache argument as above: `j` outer / `i` inner keeps
+        # reads contiguous instead of stride-`n` across columns.
+        s = ones(n)
+        @inbounds for j in 1:d
+            aj = a[j]
+            @simd for i in 1:n
+                s[i] += aj * x[i, j]
             end
-            y[i] = s^(-(d + 1))
         end
+        @inbounds @. y = s^(-(d + 1))
 
     elseif f.kind == :gaussian_peak
-        for i in 1:n
-            s = 0.0
-            for j in 1:d
-                s += a[j]^2 * (x[i, j] - u[j])^2
+        # Column reduction avoids the stride-`n` access pattern of a row-major
+        # loop nest on Julia's column-major arrays.
+        s = zeros(n)
+        @inbounds for j in 1:d
+            aj2 = a[j]^2
+            uj = u[j]
+            @simd for i in 1:n
+                s[i] += aj2 * (x[i, j] - uj)^2
             end
-            y[i] = exp(-s)
         end
+        @. y = exp(-s)
 
     elseif f.kind == :continuous
-        for i in 1:n
-            s = 0.0
-            for j in 1:d
-                s += a[j] * abs(x[i, j] - u[j])
+        # Same column-major tip: loop over columns outside, rows inside, so each
+        # column is read contiguously and accumulated into a length-`n` vector.
+        s = zeros(n)
+        @inbounds for j in 1:d
+            aj = a[j]
+            uj = u[j]
+            @simd for i in 1:n
+                s[i] += aj * abs(x[i, j] - uj)
             end
-            y[i] = exp(-s)
         end
+        @. y = exp(-s)
 
     elseif f.kind == :discontinuous
-        for i in 1:n
-            discont = false
-            s = 0.0
-            for j in 1:d
-                if x[i, j] > u[j]
-                    discont = true
-                    break
+        # Preserve the early-stop semantics with a per-row activity mask while
+        # iterating in column-major order; once a row crosses the discontinuity,
+        # later coordinates are ignored exactly as in the row loop.
+        active = trues(n)
+        s = zeros(n)
+        @inbounds for j in 1:d
+            aj = a[j]
+            uj = u[j]
+            for i in 1:n
+                if active[i]
+                    xij = x[i, j]
+                    if xij > uj
+                        active[i] = false
+                    else
+                        s[i] += aj * xij
+                    end
                 end
-                s += a[j] * x[i, j]
             end
-            y[i] = discont ? 0.0 : exp(s)
+        end
+        @inbounds for i in 1:n
+            y[i] = active[i] ? exp(s[i]) : 0.0
         end
     end
 
diff --git a/src/integrand/hartmann6d.jl b/src/integrand/hartmann6d.jl
index 46b2e8f..e27fc3e 100644
--- a/src/integrand/hartmann6d.jl
+++ b/src/integrand/hartmann6d.jl
@@ -37,17 +37,37 @@ const _HARTMANN_P = [
 
 function evaluate(f::Hartmann6D, x::AbstractMatrix)
     n = size(x, 1)
+    s1 = zeros(Float64, n)
+    s2 = zeros(Float64, n)
+    s3 = zeros(Float64, n)
+    s4 = zeros(Float64, n)
+    # Julia matrices are column-major, so accumulate one coordinate at a time
+    # rather than striding by `n` across a per-row inner loop.
+    @inbounds for j in 1:6
+        a1 = _HARTMANN_A[1, j]
+        a2 = _HARTMANN_A[2, j]
+        a3 = _HARTMANN_A[3, j]
+        a4 = _HARTMANN_A[4, j]
+        p1 = _HARTMANN_P[1, j]
+        p2 = _HARTMANN_P[2, j]
+        p3 = _HARTMANN_P[3, j]
+        p4 = _HARTMANN_P[4, j]
+        @simd for i in 1:n
+            xij = x[i, j]
+            s1[i] += a1 * (xij - p1)^2
+            s2[i] += a2 * (xij - p2)^2
+            s3[i] += a3 * (xij - p3)^2
+            s4[i] += a4 * (xij - p4)^2
+        end
+    end
     y = Vector{Float64}(undef, n)
     @inbounds for i in 1:n
-        val = 0.0
-        for k in 1:4
-            inner = 0.0
-            for j in 1:6
-                inner += _HARTMANN_A[k, j] * (x[i, j] - _HARTMANN_P[k, j])^2
-            end
-            val -= _HARTMANN_ALPHA[k] * exp(-inner)
-        end
-        y[i] = val
+        y[i] = -(
+            _HARTMANN_ALPHA[1] * exp(-s1[i]) +
+            _HARTMANN_ALPHA[2] * exp(-s2[i]) +
+            _HARTMANN_ALPHA[3] * exp(-s3[i]) +
+            _HARTMANN_ALPHA[4] * exp(-s4[i])
+        )
     end
     return y
 end
diff --git a/src/integrand/linear0.jl b/src/integrand/linear0.jl
index 5d77b1b..70d7f82 100644
--- a/src/integrand/linear0.jl
+++ b/src/integrand/linear0.jl
@@ -17,13 +17,17 @@ Linear0(tm::AbstractTrueMeasure) = Linear0(tm, tm.dimension)
 
 function evaluate(f::Linear0, x::AbstractMatrix)
     n = size(x, 1)
-    y = Vector{Float64}(undef, n)
-    @inbounds for i in 1:n
-        s = 0.0
-        for j in 1:f.dimension
-            s += x[i, j]
+    y = zeros(Float64, n)
+    # Julia stores matrices in column-major order, so a loop nest like
+    # `for i in 1:n; for j in 1:d; x[i, j] end end` strides by `n` in the inner
+    # loop and starts thrashing cache at large `d`. Swapping the loop nest to
+    # `j` outer / `i` inner reads each column contiguously while accumulating
+    # into a length-`n` vector. For each row `i`, the additions still happen in
+    # order `j = 1:d`, so the result matches the original formulation.
+    @inbounds for j in 1:f.dimension
+        @simd for i in 1:n
+            y[i] += x[i, j]
         end
-        y[i] = s
     end
     return y
 end
diff --git a/src/integrand/sensitivity_indices.jl b/src/integrand/sensitivity_indices.jl
index 5bf7b26..5db39b8 100644
--- a/src/integrand/sensitivity_indices.jl
+++ b/src/integrand/sensitivity_indices.jl
@@ -100,13 +100,16 @@ function evaluate(si::SensitivityIndices, x::AbstractMatrix)
 
     # Output: n × 3 × k
     y = Array{Float64}(undef, n, 3, k)
+    V = copy(Z)
 
     for j in 1:k
-        # Build V_j: X on subset, Z on complement
-        V = copy(Z)
+        # Build V_j: X on subset, Z on complement. Reuse one work matrix and
+        # update/reset whole columns so the mutation stays column-major.
         for col in 1:d
             if si.indices[j, col]
                 V[:, col] .= X[:, col]
+            else
+                V[:, col] .= Z[:, col]
             end
         end
         f_V = evaluate(si.base_integrand, V)
diff --git a/src/integrand/umbridge_wrapper.jl b/src/integrand/umbridge_wrapper.jl
index 490f693..02df057 100644
--- a/src/integrand/umbridge_wrapper.jl
+++ b/src/integrand/umbridge_wrapper.jl
@@ -31,7 +31,7 @@ docker run -d -p 4243:4243 linusseelinger/benchmark-muq-beam-propagation:latest
 using QMC
 import QMC: Uniform
 
-dd = DigitalNetB2(3; seed=7)
+dd = DigitalNetB2(3; seed=7, randomize="LMS_DS", graycode=false)
 tm = Uniform(dd; lower_bound=1.0, upper_bound=1.05)
 f  = UMBridgeWrapper(tm;
          url="http://localhost:4243",
diff --git a/src/kernel/dig_shift_invariant.jl b/src/kernel/dig_shift_invariant.jl
index 03cb125..3bf900f 100644
--- a/src/kernel/dig_shift_invariant.jl
+++ b/src/kernel/dig_shift_invariant.jl
@@ -87,9 +87,11 @@ function compute_kernel_eigenvalues(kernel::KernelDigShiftInvar, x::AbstractMatr
     @assert ispow2(n) "Number of points must be a power of 2 for WHT, got $n"
     # Build first column: K(x_i ⊕ x_1)
     first_col = ones(n)
-    for i in 1:n
-        for j in 1:d
-            diff = _xor_float(x[i, j], x[1, j])
+    # Column-major reduction: update all rows for one coordinate at a time.
+    for j in 1:d
+        x1j = x[1, j]
+        for i in 1:n
+            diff = _xor_float(x[i, j], x1j)
             first_col[i] *= _dig_kernel_1d(kernel.order, diff)
         end
     end
diff --git a/src/kernel/shift_invariant.jl b/src/kernel/shift_invariant.jl
index d414174..530d4c1 100644
--- a/src/kernel/shift_invariant.jl
+++ b/src/kernel/shift_invariant.jl
@@ -65,9 +65,11 @@ function compute_kernel_eigenvalues(kernel::KernelShiftInvar, x::AbstractMatrix)
     # Build the first column of the circulant kernel matrix:
     #   K_col[i] = K(x_i - x_1) = ∏_j K_1d(x_{i,j} - x_{1,j})
     first_col = ones(n)
-    for i in 1:n
-        for j in 1:d
-            diff = x[i, j] - x[1, j]
+    # Column-major reduction: update all rows for one coordinate at a time.
+    for j in 1:d
+        x1j = x[1, j]
+        for i in 1:n
+            diff = x[i, j] - x1j
             first_col[i] *= _kernel_1d_shift(kernel.order, diff)
         end
     end
diff --git a/src/stopping_criterion/cub_qmc_lattice_g.jl b/src/stopping_criterion/cub_qmc_lattice_g.jl
index 0cc32a7..6a53415 100644
--- a/src/stopping_criterion/cub_qmc_lattice_g.jl
+++ b/src/stopping_criterion/cub_qmc_lattice_g.jl
@@ -72,9 +72,35 @@ function integrate(sc::CubQMCLatticeG; resume::Union{Nothing, Dict{Symbol, Any}}
         n_iter += 1
         estimates = Vector{Float64}(undef, R)
 
-        for r in 1:R
-            y = sample_and_evaluate(sc.integrand, n)
-            estimates[r] = mean(y)
+        dd = sc.integrand.true_measure.dd
+        # Batch the R replicates' transform + evaluate in GROUPS of `group_size`
+        # rather than one `sample_and_evaluate` call per replicate. The same R
+        # randomizations are drawn in the same order (gen_samples advances the
+        # same RNG), so per-replicate means, mu_hat, and the error bound are
+        # unchanged. Grouping does one large BLAS GEMM per group instead of R
+        # tiny ones while bounding the dense-transform temporaries to
+        # group_size*n rows.
+        group_size = 4
+        r0 = 1
+        while r0 <= R
+            g = min(group_size, R - r0 + 1)
+            first = gen_samples(dd, n)
+            first =
+                ndims(first) == 3 ?
+                reshape(first, size(first, 1) * size(first, 2), size(first, 3)) : first
+            m = size(first, 1)
+            x_group = Matrix{Float64}(undef, g * m, size(first, 2))
+            @inbounds x_group[1:m, :] .= first
+            @inbounds for k in 2:g
+                xu = gen_samples(dd, n)
+                xu = ndims(xu) == 3 ? reshape(xu, size(xu, 1) * size(xu, 2), size(xu, 3)) : xu
+                x_group[((k - 1) * m + 1):(k * m), :] .= xu
+            end
+            y_group = evaluate(sc.integrand, transform(sc.integrand.true_measure, x_group))
+            @inbounds for k in 1:g
+                estimates[r0 + k - 1] = mean(@view y_group[((k - 1) * m + 1):(k * m)])
+            end
+            r0 += g
         end
 
         mu_hat = mean(estimates)
diff --git a/src/stopping_criterion/cub_qmc_net_g.jl b/src/stopping_criterion/cub_qmc_net_g.jl
index 4d469e0..2798299 100644
--- a/src/stopping_criterion/cub_qmc_net_g.jl
+++ b/src/stopping_criterion/cub_qmc_net_g.jl
@@ -1,18 +1,34 @@
 """
     CubQMCNetG(integrand; abs_tol=0.01, rel_tol=0.0,
-               n_init=2^10, n_max=2^30, n_reps=16, alpha=0.01,
+               n_init=2^10, n_max=2^30, r_lag=4,
                trace_iterations=false)
 
-Guaranteed QMC cubature using replicated randomized digital nets.
+Guaranteed QMC cubature on a **single randomized digital net** (no replications),
+porting QMCPy's `CubQMCNetG`. The error bound comes from the decay of the
+orthonormal Walsh coefficients of the integrand rather than a replication-based
+confidence interval, so it samples far fewer points than the replicated
+[`CubQMCNetGRep`](@ref).
 
-Supports **resume** and optional **iteration logging** (`trace_iterations=true`).
+At each step the net of size `n = 2^m` is drawn, the integrand is evaluated, the
+scaled Walsh coefficients `ytilde` are formed, their indices are decay-ordered
+(`kappanumap`), and the half-width bound is
+
+    fudge(m) · Σ |ytilde[kappanumap[nllstart : 2·nllstart]]|,
+    fudge(m) = 5·2⁻ᵐ,   nllstart = 2^(m − r_lag − 1).
+
+Sampling doubles until the bound meets the tolerance or `n_max` is reached.
+
+The discrete distribution must be a non-replicated [`DigitalNetB2`](@ref) in
+**natural (radical-inverse) order**, i.e. constructed with `graycode=false` and
+`replications=nothing` (matching QMCPy's "RADICAL INVERSE" requirement). A
+randomized net (e.g. `randomize="LMS_DS"`) is expected.
 
 # Example
 ```julia
-dd = DigitalNetB2(3; randomize="LMS_DS")
-tm = Uniform(dd)
-f = Genz(tm; kind=:oscillatory)
-sc = CubQMCNetG(f; abs_tol=1e-4)
+dd = DigitalNetB2(3; randomize="LMS_DS", graycode=false, seed=7)
+tm = Gaussian(dd; covariance=0.5)
+f = Keister(tm)
+sc = CubQMCNetG(f; abs_tol=1e-3)
 result = integrate(sc)
 ```
 """
@@ -22,8 +38,7 @@ mutable struct CubQMCNetG{I <: AbstractIntegrand} <: AbstractStoppingCriterion
     rel_tol::Float64
     n_init::Int
     n_max::Int
-    n_reps::Int
-    alpha::Float64
+    r_lag::Int
     trace_iterations::Bool
 end
 
@@ -33,68 +48,76 @@ function CubQMCNetG(
     rel_tol::Float64=0.0,
     n_init::Int=2^10,
     n_max::Int=2^30,
-    n_reps::Int=16,
-    alpha::Float64=0.01,
+    r_lag::Int=4,
     trace_iterations::Bool=false,
 )
-    return CubQMCNetG(
-        integrand,
-        abs_tol,
-        rel_tol,
-        n_init,
-        n_max,
-        n_reps,
-        alpha,
-        trace_iterations,
+    dd = integrand.true_measure.dd
+    dd isa DigitalNetB2 || error("CubQMCNetG requires a DigitalNetB2 discrete distribution.")
+    dd.replications === nothing || error(
+        "CubQMCNetG requires a non-replicated net (construct DigitalNetB2 with replications=nothing).",
+    )
+    !dd.graycode || error(
+        "CubQMCNetG requires natural (radical-inverse) order: construct DigitalNetB2 with graycode=false.",
     )
+    ispow2(n_init) || error("CubQMCNetG: n_init must be a power of two.")
+    n_floor = 2^max(r_lag + 1, 8)
+    if n_init < n_floor
+        @warn "CubQMCNetG: n_init bumped up to $n_floor (the guaranteed bound needs n_init ≥ 2^max(r_lag+1, 8))."
+        n_init = n_floor
+    end
+    return CubQMCNetG(integrand, abs_tol, rel_tol, n_init, n_max, r_lag, trace_iterations)
 end
 
 function integrate(sc::CubQMCNetG; resume::Union{Nothing, Dict{Symbol, Any}}=nothing)
-    R = sc.n_reps
-    t_crit = quantile(TDist(R - 1), 1.0 - sc.alpha / 2.0)
-
-    if resume !== nothing
-        n_prev = haskey(resume, :n_per_rep) ? Int(resume[:n_per_rep]) : Int(resume[:n])
-        n = 2 * n_prev
-        prev_time = Float64(get(resume, :time_integrate, 0.0))
-    else
-        n = sc.n_init
-        prev_time = 0.0
-    end
+    r_lag = sc.r_lag
+    n = resume !== nothing ? 2 * Int(get(resume, :n, sc.n_init)) : sc.n_init
+    prev_time = resume !== nothing ? Float64(get(resume, :time_integrate, 0.0)) : 0.0
 
     t_start = time()
+    dd = sc.integrand.true_measure.dd
+
     mu_hat = 0.0
     err = Inf
     n_iter = 0
+    n_final = n
     log = IterationLog()
 
     while n <= sc.n_max
         n_iter += 1
-        estimates = Vector{Float64}(undef, R)
+        m = round(Int, log2(n))
 
-        for r in 1:R
-            y = sample_and_evaluate(sc.integrand, n)
-            estimates[r] = mean(y)
+        # One fresh randomized net per doubling (single net, no replications).
+        x_unit = gen_samples(dd, n)
+        if ndims(x_unit) == 3
+            size(x_unit, 1) == 1 || error("CubQMCNetG requires a non-replicated net.")
+            x_unit = reshape(x_unit, size(x_unit, 2), size(x_unit, 3))
         end
+        y = evaluate(sc.integrand, transform(sc.integrand.true_measure, x_unit))
 
-        mu_hat = mean(estimates)
-        sigma_reps = std(estimates; corrected=true)
-        err = t_crit * sigma_reps / sqrt(R)
-        tol = max(sc.abs_tol, sc.rel_tol * abs(mu_hat))
+        # Estimate and guaranteed half-width from the Walsh-coefficient decay.
+        ytilde = _ytilde_init(y)
+        mu_hat = mean(y)
+        kappanumap = _update_kappanumap!(collect(0:(n - 1)), ytilde, m - 1, 0, m)
+
+        mllstart = m - r_lag - 1
+        nllstart = 2^mllstart
+        fudge = 5.0 * 2.0^(-m)
+        s = 0.0
+        @inbounds for p in (nllstart + 1):(2 * nllstart)
+            s += abs(ytilde[kappanumap[p] + 1])
+        end
+        err = fudge * s
 
+        tol = max(sc.abs_tol, sc.rel_tol * abs(mu_hat))
+        n_final = n
         if sc.trace_iterations
-            push!(
-                log;
-                n=n*R,
-                solution=mu_hat,
-                error_bound=err,
-                tol=tol,
-                elapsed=time() - t_start,
-            )
+            push!(log; n=n, solution=mu_hat, error_bound=err, tol=tol, elapsed=time() - t_start)
         end
 
         err <= tol && break
-        n = min(2n, sc.n_max + 1)
+        # Only ever sample at powers of two; the `while` guard stops once the
+        # next doubling would exceed n_max (so n_max need not be a power of two).
+        n *= 2
     end
 
     converged = err <= max(sc.abs_tol, sc.rel_tol * abs(mu_hat))
@@ -105,10 +128,10 @@ function integrate(sc::CubQMCNetG; resume::Union{Nothing, Dict{Symbol, Any}}=not
     t_elapsed = time() - t_start
     n_per_rep = n > sc.n_max ? sc.n_max : n
     data = Dict{Symbol, Any}(
-        :n => n_per_rep,
-        :n_per_rep => n_per_rep,
-        :n_total => n_per_rep * R,
-        :n_reps => R,
+        :n => n_final,
+        :n_per_rep => n_final,
+        :n_total => n_final,
+        :n_reps => 1,
         :error_bound => err,
         :n_iterations => n_iter,
         :converged => converged,
@@ -120,6 +143,7 @@ function integrate(sc::CubQMCNetG; resume::Union{Nothing, Dict{Symbol, Any}}=not
     return QMCResult(mu_hat, data)
 end
 
-function Base.show(io::IO, sc::CubQMCNetG)
-    print(io, "CubQMCNetG(abs_tol=$(sc.abs_tol), n_reps=$(sc.n_reps))")
-end
+Base.show(io::IO, sc::CubQMCNetG) = print(io, "CubQMCNetG(abs_tol=$(sc.abs_tol))")
+
+# Backward-compatible alias for the earlier explicit single-net name.
+const CubQMCNetGSingle = CubQMCNetG
diff --git a/src/stopping_criterion/cub_qmc_net_g_rep.jl b/src/stopping_criterion/cub_qmc_net_g_rep.jl
new file mode 100644
index 0000000..b7e9335
--- /dev/null
+++ b/src/stopping_criterion/cub_qmc_net_g_rep.jl
@@ -0,0 +1,156 @@
+"""
+    CubQMCNetGRep(integrand; abs_tol=0.01, rel_tol=0.0,
+                  n_init=2^10, n_max=2^30, n_reps=16, alpha=0.01,
+                  trace_iterations=false)
+
+Guaranteed QMC cubature using replicated randomized digital nets.
+
+This is Julia's replicated companion to the single-net [`CubQMCNetG`](@ref).
+It estimates the mean from `n_reps` independently randomized digital nets and
+uses a Student's t half-width across the replication means as its stopping rule.
+
+Supports **resume** and optional **iteration logging** (`trace_iterations=true`).
+
+# Example
+```julia
+dd = DigitalNetB2(3; randomize="LMS_DS", seed=7)
+tm = Uniform(dd)
+f = Genz(tm; kind=:oscillatory)
+sc = CubQMCNetGRep(f; abs_tol=1e-4)
+result = integrate(sc)
+```
+"""
+mutable struct CubQMCNetGRep{I <: AbstractIntegrand} <: AbstractStoppingCriterion
+    integrand::I
+    abs_tol::Float64
+    rel_tol::Float64
+    n_init::Int
+    n_max::Int
+    n_reps::Int
+    alpha::Float64
+    trace_iterations::Bool
+end
+
+function CubQMCNetGRep(
+    integrand::AbstractIntegrand;
+    abs_tol::Float64=0.01,
+    rel_tol::Float64=0.0,
+    n_init::Int=2^10,
+    n_max::Int=2^30,
+    n_reps::Int=16,
+    alpha::Float64=0.01,
+    trace_iterations::Bool=false,
+)
+    return CubQMCNetGRep(
+        integrand,
+        abs_tol,
+        rel_tol,
+        n_init,
+        n_max,
+        n_reps,
+        alpha,
+        trace_iterations,
+    )
+end
+
+function integrate(sc::CubQMCNetGRep; resume::Union{Nothing, Dict{Symbol, Any}}=nothing)
+    R = sc.n_reps
+    t_crit = quantile(TDist(R - 1), 1.0 - sc.alpha / 2.0)
+
+    if resume !== nothing
+        n_prev = haskey(resume, :n_per_rep) ? Int(resume[:n_per_rep]) : Int(resume[:n])
+        n = 2 * n_prev
+        prev_time = Float64(get(resume, :time_integrate, 0.0))
+    else
+        n = sc.n_init
+        prev_time = 0.0
+    end
+
+    t_start = time()
+    mu_hat = 0.0
+    err = Inf
+    n_iter = 0
+    log = IterationLog()
+
+    while n <= sc.n_max
+        n_iter += 1
+        estimates = Vector{Float64}(undef, R)
+
+        dd = sc.integrand.true_measure.dd
+        # Batch the R replicates' transform + evaluate in GROUPS of `group_size`
+        # rather than all at once. Each group's points are still the same draws in
+        # the same order (gen_samples advances the same RNG), so the per-replicate
+        # means — and hence mu_hat and the error bound — are unchanged. Grouping
+        # keeps one large BLAS GEMM per group (far better than R tiny ones) while
+        # bounding the dense-transform temporaries (`z`, `y`) to group_size·n rows
+        # instead of R·n: the full-stack version allocated R×-larger transform
+        # buffers, which dominated memory for the GBM integrands.
+        group_size = 4
+        r0 = 1
+        while r0 <= R
+            g = min(group_size, R - r0 + 1)
+            first = gen_samples(dd, n)
+            first =
+                ndims(first) == 3 ?
+                reshape(first, size(first, 1) * size(first, 2), size(first, 3)) : first
+            m = size(first, 1)
+            x_group = Matrix{Float64}(undef, g * m, size(first, 2))
+            @inbounds x_group[1:m, :] .= first
+            @inbounds for k in 2:g
+                xu = gen_samples(dd, n)
+                xu = ndims(xu) == 3 ? reshape(xu, size(xu, 1) * size(xu, 2), size(xu, 3)) : xu
+                x_group[((k - 1) * m + 1):(k * m), :] .= xu
+            end
+            y_group = evaluate(sc.integrand, transform(sc.integrand.true_measure, x_group))
+            @inbounds for k in 1:g
+                estimates[r0 + k - 1] = mean(@view y_group[((k - 1) * m + 1):(k * m)])
+            end
+            r0 += g
+        end
+
+        mu_hat = mean(estimates)
+        sigma_reps = std(estimates; corrected=true)
+        err = t_crit * sigma_reps / sqrt(R)
+        tol = max(sc.abs_tol, sc.rel_tol * abs(mu_hat))
+
+        if sc.trace_iterations
+            push!(
+                log;
+                n=n*R,
+                solution=mu_hat,
+                error_bound=err,
+                tol=tol,
+                elapsed=time() - t_start,
+            )
+        end
+
+        err <= tol && break
+        n = min(2n, sc.n_max + 1)
+    end
+
+    converged = err <= max(sc.abs_tol, sc.rel_tol * abs(mu_hat))
+    if !converged
+        @warn "CubQMCNetGRep: did not converge within n_max=$(sc.n_max)."
+    end
+
+    t_elapsed = time() - t_start
+    n_per_rep = n > sc.n_max ? sc.n_max : n
+    data = Dict{Symbol, Any}(
+        :n => n_per_rep,
+        :n_per_rep => n_per_rep,
+        :n_total => n_per_rep * R,
+        :n_reps => R,
+        :error_bound => err,
+        :n_iterations => n_iter,
+        :converged => converged,
+        :time_integrate => prev_time + t_elapsed,
+    )
+    if sc.trace_iterations
+        data[:iteration_log] = log
+    end
+    return QMCResult(mu_hat, data)
+end
+
+function Base.show(io::IO, sc::CubQMCNetGRep)
+    print(io, "CubQMCNetGRep(abs_tol=$(sc.abs_tol), n_reps=$(sc.n_reps))")
+end
diff --git a/src/true_measure/acceptance_rejection.jl b/src/true_measure/acceptance_rejection.jl
index bccbd6b..ac8f57e 100644
--- a/src/true_measure/acceptance_rejection.jl
+++ b/src/true_measure/acceptance_rejection.jl
@@ -91,9 +91,13 @@ function transform(tm::AcceptanceRejection, x::AbstractMatrix)
     d = tm.dimension
     M = tm.envelope_multiplier
 
-    accepted = Vector{Vector{Float64}}()
-    sizehint!(accepted, div(n, 2))
-
+    # Accumulate accepted candidates row-major into one flat buffer instead of a
+    # Vector{Vector{Float64}} + reduce(vcat, [a' …]): that avoided a per-row
+    # `collect`, an array of adjoints, and a quadratic-ish vcat. The candidate
+    # computation and acceptance test are unchanged.
+    flat = Float64[]
+    sizehint!(flat, (d * n) ÷ 2)
+    m = 0
     @inbounds for i in 1:n
         # Candidate point (first d coordinates)
         candidate = tm.proposal_sample_func(@view(x[i, 1:d]))
@@ -103,14 +107,17 @@ function transform(tm::AcceptanceRejection, x::AbstractMatrix)
         # Accept if u ≤ pdf(candidate) / (M * proposal_pdf(candidate))
         ratio = tm.pdf_func(candidate) / (M * tm.proposal_pdf_func(candidate))
         if u <= ratio
-            push!(accepted, collect(candidate))
+            for k in 1:d
+                push!(flat, candidate[k])
+            end
+            m += 1
         end
     end
 
-    if isempty(accepted)
-        return Matrix{Float64}(undef, 0, d)
-    end
-    return reduce(vcat, [a' for a in accepted])
+    m == 0 && return Matrix{Float64}(undef, 0, d)
+    # `flat` holds m rows of length d in row-major order; reshape to d×m (each
+    # column is an accepted row) and transpose to the (m, d) result.
+    return permutedims(reshape(flat, d, m))
 end
 
 function Base.show(io::IO, tm::AcceptanceRejection)
@@ -226,25 +233,32 @@ function transform(tm::AcceptanceRejectionReal, x::AbstractMatrix)
     L = tm.upper_bound
     epsq = 1e-8
 
-    accepted = Vector{Vector{Float64}}()
-    sizehint!(accepted, ceil(Int, n * tm.acceptance_rate))
-    z = Vector{Float64}(undef, d)
-
-    @inbounds for i in 1:n
-        for j in 1:d
+    # Julia matrices are column-major, so map each proposal coordinate in
+    # column-major order first, then perform the row-wise acceptance test.
+    zmat = Matrix{Float64}(undef, n, d)
+    @inbounds for j in 1:d
+        invcdf = tm.inv_cdfs[j]
+        @simd for i in 1:n
             uj = clamp(x[i, j], epsq, 1.0 - epsq)
-            z[j] = tm.inv_cdfs[j](uj)
+            zmat[i, j] = invcdf(uj)
         end
+    end
+
+    # Acceptance test reuses the already-materialised `zmat`: record the indices
+    # of accepted rows, then gather them in one shot (`zmat[keep, :]`). This drops
+    # the per-row `copy(z)`, the Vector{Vector{Float64}}, and the reduce(vcat).
+    keep = Int[]
+    sizehint!(keep, ceil(Int, n * tm.acceptance_rate))
+    @inbounds for i in 1:n
+        z = @view zmat[i, :]
         u = x[i, d + 1]
         if tm.target_density(z) >= L * tm.H_func(z) * u
-            push!(accepted, copy(z))
+            push!(keep, i)
         end
     end
 
-    if isempty(accepted)
-        return Matrix{Float64}(undef, 0, d)
-    end
-    return reduce(vcat, [a' for a in accepted])
+    isempty(keep) && return Matrix{Float64}(undef, 0, d)
+    return zmat[keep, :]
 end
 
 function Base.show(io::IO, tm::AcceptanceRejectionReal)
diff --git a/src/true_measure/geometric_brownian_motion.jl b/src/true_measure/geometric_brownian_motion.jl
index 65b870a..280626c 100644
--- a/src/true_measure/geometric_brownian_motion.jl
+++ b/src/true_measure/geometric_brownian_motion.jl
@@ -88,7 +88,8 @@ function transform(tm::GeometricBrownianMotion, x::AbstractMatrix)
     bm_samples = transform(tm._bm, x)
     # S(t) = S₀ exp[(γ - σ²/2) t + BM(t)]
     # where BM(t) already has variance σ² t from the covariance
-    # offsets is 1×d; broadcast fuses + and exp into one SIMD pass over the n×d array.
+    # offsets is 1×d; broadcast fuses + and exp into one 
+    # Single Instruction, Multiple Data (SIMD) pass over the n×d array.
     offsets = transpose((tm.drift - 0.5 * tm.diffusion) .* tm.time_vector)
     return @. tm.initial_value * exp(offsets + bm_samples)
 end
diff --git a/src/true_measure/johnsons_su.jl b/src/true_measure/johnsons_su.jl
index 2732994..2fd62ca 100644
--- a/src/true_measure/johnsons_su.jl
+++ b/src/true_measure/johnsons_su.jl
@@ -35,16 +35,17 @@ function JohnsonsSU(dd::AbstractDiscreteDistribution; xi=0.0, lambda=1.0, gamma=
 end
 
 function transform(tm::JohnsonsSU, x::AbstractMatrix)
-    n, d = size(x)
-    ndist = Distributions.Normal()
-    y = Matrix{Float64}(undef, n, d)
-    @inbounds for j in 1:d
-        for i in 1:n
-            z = quantile(ndist, _open_unit_interval(x[i, j]))
-            y[i, j] = tm.xi[j] + tm.lambda[j] * sinh((z - tm.gamma[j]) / tm.delta[j])
-        end
-    end
-    return y
+    # Phi^-1(u) = sqrt(2) * erfinv(2u - 1): the same fast standard-normal quantile the Gaussian
+    # transform uses, applied over the whole n x d input as a fused broadcast rather
+    # than a per-element Distributions.quantile call. The Johnson S_u map
+    # y = xi + lambda * sinh((Phi^-1(u) - gamma)/delta) is then broadcast column-wise
+    # (parameters are 1 x d rows). Numerically equivalent to the scalar form.
+    xi = transpose(tm.xi)
+    lam = transpose(tm.lambda)
+    gam = transpose(tm.gamma)
+    del = transpose(tm.delta)
+    z = @. sqrt(2.0) * SpecialFunctions.erfinv(2.0 * _open_unit_interval(x) - 1.0)
+    return @. xi + lam * sinh((z - gam) / del)
 end
 
 Base.show(io::IO, tm::JohnsonsSU) = print(io, "JohnsonsSU(d=$(tm.dimension))")
diff --git a/src/true_measure/student_t.jl b/src/true_measure/student_t.jl
index fb9c544..a5fa973 100644
--- a/src/true_measure/student_t.jl
+++ b/src/true_measure/student_t.jl
@@ -33,11 +33,27 @@ end
 
 function transform(tm::StudentT, x::AbstractMatrix)
     n, d = size(x)
-    tdist = Distributions.TDist(tm.df)
     y = Matrix{Float64}(undef, n, d)
-    @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]))
+    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
+        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
     end
     return y
diff --git a/src/util/fwht.jl b/src/util/fwht.jl
index 97151b7..1b2d117 100644
--- a/src/util/fwht.jl
+++ b/src/util/fwht.jl
@@ -127,3 +127,85 @@ function fwht_2d!(X::AbstractMatrix{T}) where {T <: Real}
     end
     return X
 end
+
+"""
+    _fwht_ortho(x::AbstractVector) -> Vector{Float64}
+
+Orthonormal (norm-preserving) natural-order fast Walsh-Hadamard transform:
+`H_n*x / sqrt(n)`. This is the convention used by QMCPy's `fwht` (each
+butterfly step divides by `sqrt(2)`), as opposed to the unnormalised [`fwht`](@ref)
+(`H_n*x`). It is the fast transform underlying the guaranteed digital-net
+cubature error bound.
+
+The doubling identity holds with all-ones weights: for `x = [x_1; x_2]` with
+equal halves, `_fwht_ortho(x) == [y_1 .+ y_2; y_1 .- y_2] ./ sqrt(2)` where
+`y_i = _fwht_ortho(x_i)`.
+"""
+function _fwht_ortho(x::AbstractVector)
+    n = length(x)
+    return fwht(x) ./ sqrt(n)
+end
+
+"""
+    _ytilde_init(y) -> Vector{Float64}
+
+Scaled orthonormal Walsh coefficients of the function values `y` over the first
+`length(y)` net points: `_fwht_ortho(y) / √length(y)`. The first entry equals
+`mean(y)`. Used to initialise the guaranteed digital-net cubature.
+"""
+_ytilde_init(y::AbstractVector) = _fwht_ortho(y) ./ sqrt(length(y))
+
+"""
+    _ytilde_double(ytilde_prev, ynext) -> Vector{Float64}
+
+Incrementally extend the scaled Walsh coefficients when the sample size doubles,
+reusing `ytilde_prev` (coefficients of the existing points) and the function
+values `ynext` at the newly added points (equal in number). Uses the all-ones
+FWHT doubling merge, so the result equals `_ytilde_init([y_prev; ynext])` without
+recomputing the transform over the existing points.
+"""
+function _ytilde_double(ytilde_prev::AbstractVector, ynext::AbstractVector)
+    N = length(ynext)
+    yt_omega = _fwht_ortho(ynext) ./ sqrt(N)
+    return vcat((ytilde_prev .+ yt_omega) ./ 2, (ytilde_prev .- yt_omega) ./ 2)
+end
+
+"""
+    _update_kappanumap!(kappanumap, ytilde, mfrom, mto, m) -> kappanumap
+
+In-place port of QMCPy's coefficient decay-ordering used by the guaranteed
+digital-net cubature. `kappanumap` is a permutation of `0:2^m-1` (0-based values,
+1-based positions); `ytilde` are the orthonormal Walsh coefficients. For each
+level `l` from `mfrom` down to `mto+1`, it compares the magnitudes of paired
+coefficients and swaps their indices so larger-magnitude coefficients move to the
+earlier (low-frequency) block - building the ordering the error bound sums over.
+The first index (a power of two) is never moved. Scalar (single-output) version.
+"""
+function _update_kappanumap!(
+    kappanumap::Vector{Int},
+    ytilde::AbstractVector,
+    mfrom::Int,
+    mto::Int,
+    m::Int,
+)
+    for l in mfrom:-1:(mto + 1)
+        nl = 2^l
+        flip = Int[]
+        for j in 0:(nl - 2)
+            v_old = kappanumap[j + 2]
+            v_new = kappanumap[nl + j + 2]
+            if abs(ytilde[v_new + 1]) > abs(ytilde[v_old + 1])
+                push!(flip, j + 1)
+            end
+        end
+        isempty(flip) && continue
+        for off in 0:(2 ^ (l + 1)):(2 ^ m - 2)
+            for f in flip
+                p_old = f + off + 1
+                p_new = nl + f + off + 1
+                kappanumap[p_old], kappanumap[p_new] = kappanumap[p_new], kappanumap[p_old]
+            end
+        end
+    end
+    return kappanumap
+end
diff --git a/test/run_notebooks.jl b/test/run_notebooks.jl
index f589522..62f687b 100644
--- a/test/run_notebooks.jl
+++ b/test/run_notebooks.jl
@@ -68,24 +68,47 @@ end
 errors = String[]
 times = Dict{String, Float64}()   # notebook name → wall-clock seconds
 clogger = CountingLogger()
+verbose = get(ENV, "QMC_NOTEBOOK_VERBOSE", "0") == "1"
 
 for nb in notebooks
-    println("\n", "="^60)
-    println("Running: ", nb)
-    println("="^60)
+    print("Running ", nb, " ... ")
     clogger.current_nb[] = nb
+    captured = ""
     elapsed = @elapsed try
-        with_logger(clogger) do
+        runner = () -> with_logger(clogger) do
             @nbinclude(joinpath(demos_dir, nb))
         end
+        if verbose
+            runner()
+        else
+            mktemp() do path, io
+                redirect_stdout(runner, io)
+                flush(io)
+                seekstart(io)
+                captured = read(io, String)
+            end
+        end
     catch e
         push!(errors, nb)
+        println("FAILED")
+        if !isempty(captured)
+            println("---- captured notebook stdout ----")
+            print(captured)
+            endswith(captured, '\n') || println()
+            println("---- end captured stdout ----")
+        end
         println("  ✗ FAILED: ", sprint(showerror, e))
         showerror(stdout, e, catch_backtrace())
         println()
     end
     times[nb] = elapsed
-    println("  ⏱  elapsed: ", fmt_duration(elapsed))
+    if nb in errors
+        println("  ⏱  elapsed: ", fmt_duration(elapsed))
+    else
+        w = get(clogger.counts, nb, 0)
+        wtxt = w > 0 ? ", $w warning(s)" : ""
+        println("ok [$(fmt_duration(elapsed))]$wtxt")
+    end
 end
 
 # ── Summary ──────────────────────────────────────────────
@@ -96,12 +119,11 @@ n_pass = length(notebooks) - length(errors)
 println("\n", "="^60)
 println("Summary")
 println("="^60)
-for nb in notebooks
-    status = nb in errors ? "✗ FAILED" : "✓ ok"
+for nb in errors
     w = get(clogger.counts, nb, 0)
     t = get(times, nb, 0.0)
     wtxt = w > 0 ? ", $w warning(s)" : ""
-    println("  $status  $nb  [$(fmt_duration(t))]$wtxt")
+    println("  ✗ FAILED  $nb  [$(fmt_duration(t))]$wtxt")
 end
 println("-"^60)
 println(
diff --git a/test/runtests.jl b/test/runtests.jl
index 283364d..4b67798 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -14,6 +14,7 @@ const TEST_FILES = [
     "test_stopping_criteria.jl",
     "test_integration.jl",
     "test_multilevel.jl",
+    "test_utils.jl",
     "test_aqua.jl",
 ]
 
diff --git a/test/test_discrete_distributions.jl b/test/test_discrete_distributions.jl
index 6947bb5..d3faa15 100644
--- a/test/test_discrete_distributions.jl
+++ b/test/test_discrete_distributions.jl
@@ -96,6 +96,19 @@
             @test_throws ArgumentError gen_samples(dd_file, 17)
         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_rand_a = Lattice(4; randomize=false, order="linear", seed=11, generating_vector=6)
         dd_rand_b = Lattice(4; randomize=false, order="linear", seed=11, generating_vector=6)
         @test dd_rand_a.gen_vector == dd_rand_b.gen_vector
diff --git a/test/test_integrands.jl b/test/test_integrands.jl
index 72181a9..e8281bf 100644
--- a/test/test_integrands.jl
+++ b/test/test_integrands.jl
@@ -198,6 +198,17 @@
         y = sample_and_evaluate(f, 10000)
         # Exact mean = a/2 = 3.5
         @test abs(mean(y) - 3.5) < 0.3
+
+        # Closed-form analytical reference values (validated vs pick-freeze MC).
+        ref = ishigami_exact(a=7.0, b=0.1)
+        @test ref.mean == 3.5
+        @test ref.variance ≈ 13.844588 atol = 1e-3
+        @test ref.closed ≈ [0.31391, 0.44241, 0.0] atol = 1e-3
+        @test ref.total ≈ [0.55759, 0.44241, 0.24368] atol = 1e-3
+        # structural: 0 ≤ Sᵢ ≤ Tᵢ ≤ 1, x₃ has no main effect, x₂ has no interactions
+        @test all(0.0 .<= ref.closed .<= ref.total .<= 1.0)
+        @test ref.closed[3] == 0.0
+        @test ref.closed[2] ≈ ref.total[2]
     end
 
     @testset "Hartmann6D" begin
diff --git a/test/test_integration.jl b/test/test_integration.jl
index 5fca038..ad932b1 100644
--- a/test/test_integration.jl
+++ b/test/test_integration.jl
@@ -21,11 +21,11 @@
     end
 
     @testset "Digital Net QMC: Genz continuous" begin
-        dd = DigitalNetB2(2; randomize="DS", seed=3000)
+        dd = DigitalNetB2(2; randomize="LMS_DS", graycode=false, seed=3000)
         tm = Uniform(dd)
         f = Genz(tm; kind=:continuous, a=[1.0, 1.0], u=[0.5, 0.5])
         exact = genz_exact(f)
-        sc = CubQMCNetG(f; abs_tol=0.1, n_init=2^10, n_reps=16)
+        sc = CubQMCNetG(f; abs_tol=0.1, n_init=2^10)
         result = integrate(sc)
         @test abs(result.solution - exact) < 1.0
         @test result.data[:n] >= 2^10
diff --git a/test/test_stopping_criteria.jl b/test/test_stopping_criteria.jl
index 2e4af4e..97d49f5 100644
--- a/test/test_stopping_criteria.jl
+++ b/test/test_stopping_criteria.jl
@@ -25,11 +25,31 @@
     end
 
     @testset "CubQMCNetG" begin
+        # Single-net guaranteed cubature (QMCPy CubQMCNetG). Requires a
+        # non-replicated DigitalNetB2 in natural (radical-inverse) order.
+        dd = DigitalNetB2(3; randomize="LMS_DS", graycode=false, seed=2024)
+        tm = Gaussian(dd; covariance=0.5)
+        f = Keister(tm)
+        exact = keister_exact(3)
+        sc = CubQMCNetG(f; abs_tol=0.01, n_init=2^10)
+        result = integrate(sc)
+        @test abs(result.solution - exact) < 0.05
+        @test result.data[:error_bound] <= 0.01 + 1e-9
+        @test result.data[:converged]
+        @test result.data[:n_reps] == 1
+        @test result.data[:n_total] == result.data[:n]
+
+        # Guard: a default (graycode=true) net is rejected at construction.
+        dd_gc = DigitalNetB2(3; randomize="LMS_DS", seed=1)
+        @test_throws ErrorException CubQMCNetG(Keister(Gaussian(dd_gc; covariance=0.5)))
+    end
+
+    @testset "CubQMCNetGRep" begin
         dd = DigitalNetB2(2; randomize="DS", seed=700)
         tm = Uniform(dd)
         f = Genz(tm; kind=:gaussian_peak, a=[1.0, 1.0], u=[0.5, 0.5])
         exact = genz_exact(f)
-        sc = CubQMCNetG(f; abs_tol=0.1, n_init=2^10, n_reps=16)
+        sc = CubQMCNetGRep(f; abs_tol=0.1, n_init=2^10, n_reps=16)
         result = integrate(sc)
         @test abs(result.solution - exact) < 1.0
         @test result.data[:n] >= 2^10
@@ -37,6 +57,15 @@
         @test result.data[:n_total] == result.data[:n] * result.data[:n_reps]
     end
 
+    @testset "CubQMCNetGSingle alias" begin
+        dd = DigitalNetB2(3; randomize="LMS_DS", graycode=false, seed=11)
+        tm = Gaussian(dd; covariance=0.5)
+        f = Keister(tm)
+        sc = CubQMCNetGSingle(f; abs_tol=0.01, n_init=2^10)
+        @test sc isa CubQMCNetG
+        @test integrate(sc).data[:n_reps] == 1
+    end
+
     @testset "CubQMCBayesLatticeG (smoke)" begin
         dd = Lattice(2; randomize=true, seed=800)
         tm = Uniform(dd)
@@ -188,6 +217,64 @@
         @test scml.rmse_tol > 0
         set_tolerance!(scml; rmse_tol=0.05)
         @test scml.rmse_tol == 0.05
+
+        scmlqmc = CubMLQMC(FinancialOptionML(gbm; d_coarsest=4); abs_tol=0.1)
+        @test scmlqmc.rmse_tol > 0
+        set_tolerance!(scmlqmc; rmse_tol=0.05)
+        @test scmlqmc.rmse_tol == 0.05
+
+        scmlcont = CubMLMCCont(FinancialOptionML(gbm; d_coarsest=4); abs_tol=0.1)
+        @test scmlcont.target_tol > 0
+        set_tolerance!(scmlcont; rmse_tol=0.05)
+        @test scmlcont.target_tol == 0.05
+
         @test_throws ArgumentError set_tolerance!(scml; rel_tol=0.01)
     end
+
+    @testset "QMC criterion contract (characterization)" begin
+        # Value-free regression guards that any future rework of the replicated
+        # QMC criteria must preserve. They assert the convergence *contract* and
+        # cross-criterion consistency rather than hard-coded estimates, so they
+        # remain valid across legitimate refactors (e.g. an incremental-net or
+        # Walsh-coefficient rewrite) while catching a criterion that silently
+        # stops converging, drops a schema key, becomes non-reproducible, or
+        # disagrees with the other QMC rule.
+        exact = keister_exact(3)
+        net =
+            tol -> CubQMCNetGRep(
+                Keister(
+                    Gaussian(DigitalNetB2(3; randomize="LMS_DS", seed=2024); covariance=0.5),
+                );
+                abs_tol=tol,
+                n_init=2^10,
+                n_reps=16,
+            )
+        lat =
+            tol -> CubQMCLatticeG(
+                Keister(Gaussian(Lattice(3; randomize=true, seed=2024); covariance=0.5));
+                abs_tol=tol,
+                n_init=2^10,
+                n_reps=16,
+            )
+
+        for make in (net, lat)
+            res = integrate(make(0.05))
+            # Converges within a safe (5x) margin of the known exact value.
+            @test abs(res.solution - exact) < 0.25
+            # Reports having met the requested tolerance.
+            @test res.data[:error_bound] <= 0.05 + 1e-9
+            # Standardized result schema is complete and self-consistent.
+            for k in (:n, :n_per_rep, :n_total, :n_reps, :error_bound)
+                @test haskey(res.data, k)
+            end
+            @test res.data[:n_total] == res.data[:n_per_rep] * res.data[:n_reps]
+            # Reproducible under a fixed seed.
+            @test integrate(make(0.05)).solution == res.solution
+            # A tighter tolerance never uses fewer samples.
+            @test integrate(make(0.005)).data[:n_total] >= res.data[:n_total]
+        end
+
+        # The two QMC rules estimate the same integral.
+        @test abs(integrate(net(0.05)).solution - integrate(lat(0.05)).solution) < 0.2
+    end
 end
diff --git a/test/test_true_measures.jl b/test/test_true_measures.jl
index 462f782..a670bf2 100644
--- a/test/test_true_measures.jl
+++ b/test/test_true_measures.jl
@@ -226,4 +226,17 @@
         @test size(y) == (100, 2)
         @test abs(mean(y[:, 1])) < 0.5
     end
+
+    @testset "Open unit interval helper" begin
+        @test QMC._open_unit_interval(0.0f0) == eps(Float32)
+        @test QMC._open_unit_interval(1.0f0) == 1.0f0 - eps(Float32)
+        @test QMC._open_unit_interval(0) == eps(Float64)
+
+        dd = IIDStdUniform(2; seed=101)
+        tm = JohnsonsSU(dd)
+        x = Float32.(gen_samples(dd, 8))
+        y = transform(tm, x)
+        @test size(y) == size(x)
+        @test all(isfinite, y)
+    end
 end
diff --git a/test/test_utils.jl b/test/test_utils.jl
index 74303f2..ad39a39 100644
--- a/test/test_utils.jl
+++ b/test/test_utils.jl
@@ -40,4 +40,119 @@
         QMC.fwht_2d!(X)
         @test X == [4.0 0.0; 0.0 0.0]
     end
+
+    @testset "to_bin / to_float" begin
+        # Binary expansions of dyadic rationals.
+        @test to_bin(0.5, 4) == [1, 0, 0, 0]
+        @test to_bin(0.25, 4) == [0, 1, 0, 0]
+        @test to_bin(0.75, 4) == [1, 1, 0, 0]
+        @test to_bin(0.625, 4) == [1, 0, 1, 0]
+        @test to_bin(0.0, 4) == [0, 0, 0, 0]
+        # to_float inverts to_bin for values with finite expansions ≤ m bits.
+        for x in (0.0, 0.25, 0.5, 0.625, 0.75)
+            @test to_float(to_bin(x, 8)) ≈ x
+        end
+        # Out-of-range inputs are rejected.
+        @test_throws AssertionError to_bin(1.0, 4)
+        @test_throws AssertionError to_bin(-0.1, 4)
+    end
+
+    @testset "bro_fft / bro_ifft" begin
+        # bro_fft = FFT of the bit-reversed input; a constant gives a pure DC term.
+        @test bro_fft([1.0, 1.0, 1.0, 1.0]) ≈ ComplexF64[4, 0, 0, 0]
+        @test bro_fft([1.0, 2.0, 3.0, 4.0]) ≈ ComplexF64[10, -1 + 1im, -4, -1 - 1im]
+        # Bit reversal is an involution, so bro_ifft inverts bro_fft.
+        x = Float64[3, 1, 4, 1, 5, 9, 2, 6]
+        @test bro_ifft(bro_fft(x)) ≈ x
+    end
+
+    @testset "periodize" begin
+        x = [0.0 0.25; 0.5 0.75]
+        @test periodize(x, :NONE) == x
+        # All smooth transforms fix the endpoints and map the midpoint to itself.
+        for pt in (:C1, :C1SIN, :C2SIN, :C3)
+            @test periodize(x, pt)[2, 1] ≈ 0.5     # x = 0.5 ↦ 0.5
+        end
+        @test periodize([0.0 1.0], :C3) ≈ [0.0 1.0]
+        # Baker's (tent) transform: 1 - |2x - 1|.
+        @test periodize([0.0 0.25 0.5 0.75], :BAKER) ≈ [0.0 0.5 1.0 0.5]
+        # Unknown transforms error.
+        @test_throws ErrorException periodize(x, :BOGUS)
+    end
+
+    @testset "FWHT (orthonormal)" begin
+        # Matches QMCPy's `fwht` convention: H_n*x / sqrt(n) (orthonormal).
+        @test QMC._fwht_ortho([1.0, 1.0]) ≈ [sqrt(2.0), 0.0]
+        @test QMC._fwht_ortho([1.0, 0.0, 0.0, 0.0]) ≈ [0.5, 0.5, 0.5, 0.5]
+        @test QMC._fwht_ortho([1.0, 2.0, 3.0, 4.0]) ≈ [5.0, -1.0, -2.0, 0.0]
+        @test QMC._fwht_ortho([1.0, 1.0, 1.0, 1.0]) ≈ [2.0, 0.0, 0.0, 0.0]
+
+        x = Float64[3, 1, 4, 1, 5, 9, 2, 6]
+        # Orthonormal => norm-preserving, and the DC coefficient equals mean*sqrt(n).
+        @test norm(QMC._fwht_ortho(x)) ≈ norm(x)
+        @test QMC._fwht_ortho(x)[1] ≈ sum(x) / sqrt(length(x))
+
+        # Doubling identity (all-ones weights): _fwht_ortho([x_1; x_2]) for equal
+        # halves equals [y_1.+y_2; y_1.-y_2]/sqrt(2). This is the incremental-update
+        # relation the guaranteed cubature relies on.
+        x1 = Float64[1, 2, 3, 4]
+        x2 = Float64[5, 6, 7, 8]
+        y1 = QMC._fwht_ortho(x1)
+        y2 = QMC._fwht_ortho(x2)
+        @test QMC._fwht_ortho(vcat(x1, x2)) ≈ vcat(y1 .+ y2, y1 .- y2) ./ sqrt(2.0)
+    end
+
+    @testset "ytilde incremental doubling" begin
+        # The incremental update must equal the orthonormal transform of the full
+        # point set - the relation that lets the cubature reuse work across
+        # doublings - and the first coefficient must track the running mean.
+        y = Float64[3, 1, 4, 1, 5, 9, 2, 6, 5, 3, 5, 8, 9, 7, 9, 3]  # 16 = 2^4
+        N = 4
+        yt = QMC._ytilde_init(y[1:N])
+        off = N
+        cur = N
+        for _ in 1:2  # 4 -> 8 -> 16
+            yt = QMC._ytilde_double(yt, y[(off + 1):(off + cur)])
+            off += cur
+            cur *= 2
+        end
+        @test yt ≈ QMC._ytilde_init(y[1:off])
+        @test yt[1] ≈ sum(y[1:off]) / off
+    end
+
+    @testset "kappanumap (coefficient ordering)" begin
+        # Reference outputs computed from QMCPy's own _update_kappanumap (scalar).
+        # The comparison uses only the magnitude ordering of ytilde, so rounded
+        # literals reproduce the exact permutation.
+        yt3 = [-0.535669, 0.361595, 1.304, 0.947081, -0.703735, -1.265421, -0.623274, 0.041326]
+        k3 = QMC._update_kappanumap!(collect(0:7), yt3, 2, 0, 3)
+        @test k3 == [0, 5, 2, 3, 4, 1, 6, 7]
+
+        yt4 = [
+            -2.325,
+            -0.2188,
+            -1.2459,
+            -0.7323,
+            -0.5443,
+            -0.3163,
+            0.4116,
+            1.0425,
+            -0.1285,
+            1.3665,
+            -0.6652,
+            0.3515,
+            0.9035,
+            0.094,
+            -0.7435,
+            -0.9217,
+        ]
+        k4 = QMC._update_kappanumap!(collect(0:15), yt4, 3, 0, 4)
+        @test k4 == [0, 9, 2, 7, 12, 5, 14, 3, 8, 1, 10, 15, 4, 13, 6, 11]
+
+        # Structural invariants: output is a permutation of 0:n-1, and the first
+        # (power-of-two) index is never moved.
+        @test sort(k3) == collect(0:7)
+        @test sort(k4) == collect(0:15)
+        @test k3[1] == 0 && k4[1] == 0
+    end
 end