pf_random_resample Interface

public interface pf_random_resample

Fills idx with size(idx) values drawn from 1 .. m with replacement.

The third member of the resampling trio, and the one whose construction is not a construction at all: drawing with replacement means size(idx) independent uniform integers in [1, m], with no dedup structure, no permutation and no sort. So this is the draw-axis integer bulk fill under another name, and the identity is exact and is asserted by the suite:

call pf_random_resample(idx, m, seed, stream)
call pf_random_fill_draws(seed, stream, idx, 1_int64, m)   ! the SAME values

What the name buys is the 1.4x-1.6x a caller loses by writing the obvious loop. Without it the natural code is do k = 1, n; idx(k) = pf_random_int_at(seed, i, 1, m, k); end do, which re-enciphers a Philox block for every value where the bulk form serves two draws from each one -- measured at 25.7 against 15.7 ns per value on machine B (gfortran 15.2.1, --profile release, 10M values) and 16.1 against 10.1 on machine A. The values are identical either way; only the route to them differs.

idx is a rank-1 integer(int32) or integer(int64) array, intent(out); m is integer(int32) or integer(int64); seed is integer(int64). stream is optional, integer(int32) or integer(int64), and defaults to 1 -- it selects which replicate this is, so replicate b is reproducible from (seed, b) alone, independent of how many replicates were asked for or in what order they ran. A zero-sized idx is a defined no-op and is not validated: it asks for nothing, so no precondition applies to it.

Two preconditions, both aborting rather than truncating or wrapping, and only when idx is non-empty: m >= 1; and -- for an integer(int32) idx only -- m <= huge(int32), since an element may be any value in [1, m] and one above huge(int32) has nowhere to go. The second is why the integer(int32)-array/integer(int64)-m pairing is accepted at compile time and rejected at run time, exactly as on pf_random_subset.

There is deliberately NO size(idx) <= m precondition, and its absence is the clearest statement of how this differs from its sibling. That bound belongs to a subset drawn without replacement; here size(idx) == m is the most ordinary bootstrap there is, and size(idx) far beyond m is perfectly meaningful. A guard copied across from pf_random_subset would refuse the procedure's main use.

stream is this procedure's replicate axis, and the siblings do not have one. pf_random_permutation and pf_random_subset are keyed by (seed, m) alone, so independent replicates of those come from pf_random_key(seed, b) instead. A resample is built on the draw axis, which carries a stream coordinate already, so it costs a caller one integer rather than a key derivation. Both routes are available here: stream = b and seed = pf_random_key(seed, b) are equally independent.

threads behaves exactly as it does on pf_random_permutation and pf_random_subset -- same parquet_set_random_threads default, same parquet_set_random_parallel_min_elements work floor, and the floor applies to an explicit request too, so a small resample stays serial however many workers are asked for. The result is bit-identical at every thread count, and by construction rather than by care: element k is a pure function of (seed, stream, k), so a thread filling elements lo .. hi is the serial fill started at draw lo. That is asserted rather than argued.

threads requires an explicit stream, and that is a language constraint rather than a choice. threads and stream are both integers in argument position 4, so a generic offering (idx, m, seed [, threads]) beside (idx, m, seed, stream) is rejected by the compiler -- "Ambiguous interfaces in generic interface 'pf_random_resample'" -- and giving the two dummies different keyword names does not rescue it, because keyword names do not make specifics distinguishable. Both spellings were compiled to confirm it. So write call pf_random_resample(idx, m, seed, 1, threads=8) for the default replicate; omitting stream produces a "no specific subroutine matches" error that does not explain itself. This is the same constraint the parquet_open_reader split answers one level along, where an optional argument differing only by kind could not disambiguate either.

Module Procedures