Fix amrex::Random() returning 1.0 in SP & add amrex::RandomPositive

[ax3l]

Fixes #5638.

Two commits: the first restores the documented contract of amrex::Random(), the second adds the generator that user code actually keeps reaching for.

1. amrex::Random() can return exactly 1.0

amrex::Random() is documented as [0,1), but on CUDA/HIP it is 1 - curand_uniform(). hiprand/curand draw from (0,1], so the flip is right in exact arithmetic and wrong in floating point: when the draw is below half an ULP of one, 1 - draw rounds up to exactly 1.0. In single precision that is a draw <= 2^-25, i.e. ~3e-8 per call (double precision: 2^-54, ~5e-17).

The loud consequence is the one from #5638 (1 - amrex::Random() is not a zero-guard). The quieter and arguably worse one is that callers relying on the strict upper bound break silently: int(N*Random()) becomes an out-of-bounds index, and problo + Random()*dx puts a particle one cell over. No inf, no assert.

Fixed by clamping to the largest representable Real below one. The constant is 1 - eps/2, which is exactly std::nextafter(Real(1), Real(0)) in both precisions — I verified this for float and double, with and without -ffast-math. Since it is constexpr, it works in device code, where std::nextafter is not available, so neither of the two workarounds discussed in #5638 is needed: no 1 - epsilon (which is 2 ULPs below one and needlessly discards a representable value), and no host-computed value stashed in RandomEngine.

The lower bound needs no clamp — 0.0 is in [0,1) by design. Host and SYCL already produce [0,1) natively and are untouched.

2. amrex::RandomPositive(), a uniform real in (0,1]

The deeper issue behind #5638 is that user code needs (0,1] at least as often as [0,1) — every -log(u), sqrt(-2*log(u)), pow(u,-a) — and AMReX offered no way to ask for it. So downstream codes invent 1 - amrex::Random(), which is exactly the transform that does not survive single precision. Documenting the hazard is necessary but not sufficient; there has to be something to point people to.

RandomPositive() costs nothing, because the flip is safe in one direction and unsafe in the other:

backend	native interval	Random() → [0,1)	RandomPositive() → (0,1]
CUDA / HIP	(0,1]	flip + clamp (commit 1)	raw generator call, no arithmetic
CPU (uniform_real_distribution)	[0,1)	free	1 - u, safe
SYCL (oneMKL uniform)	[0,1)	free	1 - u, safe

On GPU it is cheaper than Random(). On host/SYCL the flip is the safe direction: u <= 1 - eps/2 implies 1 - u >= eps/2 > 0.

Also switched RandomGamma over to it — that is AMReX's own exposure. Marsaglia–Tsang needs u in (0,1) for both the std::log(u) rejection step and the std::pow(u, 1/alpha) boost when alpha < 1.

Safety under all math modes

Both guarantees hold regardless of how AMReX is built — desktop or HPC, IEEE or relaxed:

• the clamp is a comparison, not arithmetic, so there is nothing to reassociate or contract;
• almost_one is exactly representable, so constant folding is stable under fast-math;
• the smallest value RandomPositive() can return is eps/2 (2^-24 in single precision), a normal number — flush-to-zero and denormals-are-zero cannot touch it;
• there is only ever a single subtraction from an opaque generator result, so no algebraic identity applies.

This is asserted in the code comments and in the manual, and verified empirically below under --use_fast_math -ftz=true.

Does the clamp perturb the uniformity of the draw?

No, and the reason is stronger than "the probability is small": every value inside the documented range passes through bit-for-bit unchanged. The clamp fires only on an output of exactly 1.0, which was never in [0,1) to begin with, so the distribution restricted to the contract is identical to before — the change only redirects mass that was out-of-contract.

Quantitatively, it moves 2^-25 (~3e-8) of probability by one ULP, from 1.0 onto the adjacent value 1-2^-24. That bounds the shift in any moment by 2^-25 x 2^-24 = 2^-49 (~2e-15), i.e. ~7 orders of magnitude below the 2^-24 resolution of a single-precision sample itself. In double precision the moved mass is 2^-54.

It is worth putting that next to the discretization the transform already has. Because 1 - curand_uniform() rounds onto the coarse float grid near one, amrex::Random() returns exactly 0.0 with probability ~2^-25 — the run above measured 103 zeros in 3.9e9 draws — which is vastly more than an ideal uniform-over-floats would give. The endpoints of this distribution are already lumpy at exactly the 2^-25 scale; the clamp relocates one such atom by a single ULP. It is not a new artifact, and it is smaller than what is already there.

Two things it does not touch:

• The engine stream. The clamp is a pure post-map on the output value. It consumes no extra draws and does not alter engine state, so period, equidistribution and sequence-level test-suite behaviour are unchanged.
• RandomPositive. It contains no clamp at all. On CUDA/HIP it is the raw generator output, so it applies strictly less transformation than Random() does; on host/SYCL the u >= 0.5 half of the flip is exact by Sterbenz.

Empirically, at 4e8 draws in single precision on GPU, binned into 1024 bins:

  Random()         mean=0.49999260  var=0.08333193  chi2=1035.5/1023 dof  z=+0.28
  RandomPositive() mean=0.50000655  var=0.08333565  chi2=1022.9/1023 dof  z=-0.00
                                    (exact: mean 0.5, var 0.0833334)

Both are consistent with uniform; neither shows any bias attributable to the change.

Docs

New Random Numbers section in the Basics chapter of the user's guide (sec:basics:random): how to draw inside a GPU kernel with ParallelForRNG, seeding, thread safety, a table of both uniform generators against what each is for, worked exponential and power-law examples, a table of the other distributions (RandomNormal, RandomPoisson, RandomGamma, Random_int, Random_long), and a warning:: block on the 1 - amrex::Random() anti-pattern. Cross-referenced from the FAQ entry on random numbers, and Doxygen updated on both functions.

While documenting intervals I also fixed the FillRandom doxygen, which reported its interval with typos ("SYCl", "CUADA") — see the note below.

Why not just tell people to use amrex::RandomNormal?

For a plain normal deviate, we should — the manual's distribution table lists it directly below these examples, and the doxygen carries a \see cross-reference. RandomNormal is faster (on GPU it is the vendor's own normal generator, and it does not burn two uniforms per deviate), and since it never evaluates std::log it is untouched by this whole class of bug. Worth noting that this is not a coincidence: curand_uniform excludes zero precisely because curand's own normal generator needs log(u). The 1 - c flip discarded exactly the property the vendor put there on purpose.

Hand-written Box-Muller still earns its place when the radial variable itself is needed, which RandomNormal cannot give you:

• Truncation. ImpactX's Gaussian samples a radially truncated beam by inverting the radial CDF directly, r^2 = -2 log(1 - u(1-e^-f)) with f = cut^2/2, plus a variance renormalization. Rejection sampling on RandomNormal would be an unbounded loop with divergent warps on GPU, would consume a variable number of draws (breaking reproducibility), and would give a square cut in (x,px) rather than the intended circular one — a different distribution, not just a slower route to the same one.
• Isotropic directions. Waterbag normalizes a 6-vector of deviates onto the unit sphere; there the Gaussians are a means to an end.

So Box-Muller is a legitimate RandomPositive consumer, but a poor headline example, since it invites hand-rolling a normal that AMReX already provides. The doxygen and manual examples now lead with the exponential and power-law cases, which have no AMReX alternative, and Box-Muller appears only in the note explaining when it is actually the right tool.

Verification

All on an RTX A2000 (sm_86), CUDA 13.2.

The bug, reproduced on hardware. Direct comparison of the old and new transform over real curand_uniform draws:

3932160000 draws
  Random(old)    == 1.0f : 109   (contract violation, predicted ~118)
  Random(new)    == 1.0f : 0
  Random(new)    >= 1.0f : 0
  Random(new)    == 0.0f : 103   (legal: 0 is in [0,1))
  RandomPositive == 0.0f : 0
  RandomPositive >  1.0f : 0

109 violations in 3.9e9 draws against ~118 predicted by 2^-25. Byte-identical results under -O3 and under --use_fast_math -ftz=true. The pathological sweep also shows the downstream effect directly — for curand = 2^-25, 1 - Random() goes from 0 (old) to 5.96e-08 (new), so the WarpX-style guard starts working again:

curand         Random(old)    Random(new)    1-Random(new)  RandomPositive
5.960464e-08   0.99999994     0.99999994     5.96046448e-08 5.96046448e-08
2.980232e-08   1              0.99999994     5.96046448e-08 2.98023224e-08   <== old broke here
1.490116e-08   1              0.99999994     5.96046448e-08 1.49011612e-08   <== old broke here
2.328306e-10   1              0.99999994     5.96046448e-08 2.32830644e-10   <== old broke here

Against the real API, single precision + CUDA, via ParallelForRNG:

sizeof(Real) = 4, draws = 400000000
  Random()         >= 1  : 0    (must be 0)
  Random()         == 0  : 15   (legal, and exactly why RandomPositive exists)
  RandomPositive() <= 0  : 0    (must be 0)
  RandomPositive() >  1  : 0    (must be 0)
  host draws = 20000000 -- same, all 0
  RandomGamma non-finite/non-positive : 0
RESULT: PASS

Builds and tests

# single precision + CUDA, full library build: 0 warnings, 0 errors
cmake -S . -B build-sp-cuda -DAMReX_GPU_BACKEND=CUDA -DAMReX_PRECISION=SINGLE \
      -DAMReX_PARTICLES_PRECISION=SINGLE -DAMReX_CUDA_ARCH=86 -DCMAKE_BUILD_TYPE=Release
cmake --build build-sp-cuda -j6

# CPU + tests: 0 warnings, 0 errors; 8/8 ctest passed
cmake -S . -B build-cpu-check -DAMReX_ENABLE_TESTS=ON -DAMReX_TEST_TYPE=Small -DAMReX_MPI=OFF
cmake --build build-cpu-check -j6
ctest --test-dir build-cpu-check --output-on-failure -j4     # 100% tests passed, 0 failed out of 8

# both GPU precision branches, both math modes -- compile clean
nvcc -std=c++20 -arch=sm_86 [--use_fast_math -ftz=true] ...  # DOUBLE and SINGLE: exit 0

Sphinx builds clean; the only warning is the pre-existing missing amrex.pdf download.

Notes / open questions

• Naming is bikesheddable — RandomPositive fits the Random* prefix, RandomNonZero reads more literally. Happy to change.
• Should Random() keep returning 0.0? It is legal for [0,1) and this PR keeps it. A (0,1) variant excluding both endpoints would cover everything with one function (it would also serve logistic/Cauchy inversion, which (0,1] does not), at the cost of a clamp and an epsilon on every backend. I left it out; easy to add if wanted.
• FillRandom is still backend-dependent ([0,1) on CPU/SYCL, (0,1] on CUDA/HIP), since it forwards to the vendor host-side generators. Making it consistent would need a full extra pass over the array, which is a real cost on a bulk-fill path and a call for maintainers — so this PR only documents it accurately in both the doxygen and the manual. Worth a follow-up decision.
• Downstream: ImpactX [WIP] Guard single-precision amrex::Random() edges in beam & spin distributions BLAST-ImpactX/impactx#1627 can drop its amrex::max(u1, numeric_limits::min()) clamps in favor of RandomPositive, and WarpX has four sites that want it (SampleGaussianFluxDistribution.H:49, InjectorMomentum.H:320,322, DefaultInitialization.H:71).

To Do

• draft 🤖 Generated with Claude Code
• self-review