parquet_random – parquet-fortran

Counter-based random numbers: reproducible under any OpenMP schedule, at any thread count.

An ordinary generator carries state, so which value a loop iteration receives depends on how many draws came before it -- which, in a parallel loop, depends on the schedule, the thread count and the machine's timing. This module removes the state instead of guarding it. Every value is a pure function of three coordinates: a seed, a stream index i, and a 1-based draw within that stream. Iteration 5000 gets the same number whether it ran first, last, alone, or on a machine with 384 cores.

!$omp parallel do schedule(dynamic)
do i = 1, n
    x(i) = pf_random_at(seed, i)          ! same value for this i, always
end do

The generator is Philox4x32-10, a counter-based cipher of established quality. Every value this module can produce is frozen by the contract pf_random_algorithm names: a change to any of it is a major version, and the identifier is how a program can tell.

Not cryptographic. Do not use it for keys, tokens, or anything an adversary should not be able to predict -- the seed is recoverable from a handful of outputs.

Uses

Variables

Type	Visibility	Attributes		Name		Initial
character(len=*),	public,	parameter	::	pf_random_algorithm	=	"philox4x32-10/v3"	Identifies the algorithm together with every mapping this module freezes -- the cipher, the key and counter layout, the word order, the integer rule and its retry key. Its value changes if and only if one of those changes, so a program that records it can tell whether a stored result is still reproducible. Identical on both sides of the route (e) fork, which changes how a product is formed and never what it equals. `/v2` differs from `/v1` in exactly one mapping: `pf_random_int_at`'s draw axis, which had stride 4 (a whole block per value, half of it unused) and now has stride 2, agreeing with `pf_random_at` and `pf_random_bits_at`. Draw 1 is unchanged; draws from 2 up moved. The cipher, the key and counter layout, the word order, the rejection rule and the retry key are all identical between the two.
character(len=*),	public,	parameter	::	pf_exp_algorithm	=	"exp:-log(1-u)/libm+expkey/v1"	Identifies the exponential mapping, covering both realisations: the transform, the `1 - u` convention that keeps a draw finite, which logarithm each realisation uses, and the two-word cost they share. One identifier for the pair, per the rule above -- a program recording a single string wants to know whether either form moved. The two halves after the slash name the two logarithms: `libm` for `pf_random_exp_at` (fast, bit-identical for a given libm) and `expkey` for `pf_random_exp_portable_at` (`parquet_expkey`'s frozen transform, identical on every platform and compiler).
character(len=*),	public,	parameter	::	pf_normal_algorithm	=	"normal:ziggurat256+polar/libm+expkey/v1"	Identifies the normal mapping, covering both realisations and every part of each that decides a value: the Ziggurat's layer count and the construction that solved for its tables, the polar variant, both rejection loops' draw order and word cost, and the sub-stream labels the coordinate-addressed forms derive. `ziggurat256` is `%normal`/`pf_random_normal_at` -- 256 equal-area layers from `parquet_ziggurat`, with libm `log`/`exp` in the tail and the wedge test, so bit-identical for a given libm. `polar` is `%normal_portable`/`pf_random_normal_portable_at` -- Marsaglia's polar method over `parquet_expkey`'s frozen logarithm and IEEE `sqrt`, identical on every platform. It moves if `parquet_ziggurat` is regenerated with a different layer count, because every value moves with it. It does NOT move when `pf_random_algorithm` does, and vice versa: a program recording one must not be told the other changed.
character(len=*),	public,	parameter	::	pf_gamma_algorithm	=	"gamma:marsaglia-tsang+boost/normal=ziggurat256/libm/v1"	Identifies the Gamma mapping: Marsaglia-Tsang's squeeze-and-log rejection, the `Gamma(a) = Gamma(a+1) * u(1/a)` boost that extends it below `shape = 1`, the draw order within a candidate, and -- named explicitly, because it decides every value -- which normal the inner loop consumes. The dependency has to be in the string. Gamma draws a normal per candidate, so a change to the normal changes every Gamma value; a program recording only this identifier would otherwise miss it. `normal=ziggurat256` is that statement, and it makes `pf_normal_algorithm`'s first component load-bearing here too. Bit-identical for a given libm only, and there is no portable form (Q3-9). The boost needs `u(1/a)` -- a `pow` this project has not frozen and does not intend to -- so a `%gamma_portable` could not keep the promise its name would imply, whatever the rest of the algorithm did. Do not add one for symmetry with `%normal`.
character(len=*),	public,	parameter	::	pf_poisson_algorithm	=	"poisson:knuth\|10\|ptrs/libm/v1"	Identifies the Poisson mapping: both algorithms, the crossover between them, and each one's draw order. The crossover is contract, not a tuning knob, which is why it is in the string: it decides which algorithm runs and therefore which value comes back. Below `lambda = 10` the draw is Knuth's product of uniforms, which is exact and terminates unconditionally; at or above it, Hoermann's transformed rejection with squeeze (PTRS), which is `O(1)` in `lambda` where Knuth is `O(lambda)`. Bit-identical for a given libm: PTRS needs `log` and `log_gamma`, and Knuth needs `exp`.

Interfaces

public interface pf_random_at

One uniform real64 in [0, 1): value draw (default 1) of stream i under seed.

i is integer(int32) or integer(int64); the two give identical values, because an int32 stream sign-extends before it reaches the counter. seed and draw are always integer(int64), so an explicit draw literal is written 3_int64.

private pure elemental function pf_random_at_i32(seed, i, draw) result(r)

pf_random_at for an integer(int32) stream index.

Arguments

Type	Intent	Optional		Name
integer(kind=int64),	intent(in)		::	seed	the stream family's seed
integer(kind=int32),	intent(in)		::	i	stream index; sign-extends, so any value is valid
integer(kind=int64),	intent(in),	optional	::	draw	1-based value index; absent means 1

Return Value real(kind=real64)

a uniform draw in [0, 1)

private pure elemental function pf_random_at_i64(seed, i, draw) result(r)

pf_random_at for an integer(int64) stream index.

Arguments

Type	Intent	Optional		Name
integer(kind=int64),	intent(in)		::	seed	the stream family's seed
integer(kind=int64),	intent(in)		::	i	stream index; every value is valid, negatives included
integer(kind=int64),	intent(in),	optional	::	draw	1-based value index; absent means 1

Return Value real(kind=real64)

a uniform draw in [0, 1)

public interface pf_random32_at

One uniform real32 in [0, 1): value draw (default 1) of stream i under seed.

This enumerates its OWN sequence, one word per value, and is deliberately not a narrowing of pf_random_at -- the two share a stream but not a value. i is integer(int32) or integer(int64); seed and draw are integer(int64).

private pure elemental function pf_random32_at_i32(seed, i, draw) result(r)

pf_random32_at for an integer(int32) stream index.

Arguments

Type	Intent	Optional		Name
integer(kind=int64),	intent(in)		::	seed	the stream family's seed
integer(kind=int32),	intent(in)		::	i	stream index; sign-extends, so any value is valid
integer(kind=int64),	intent(in),	optional	::	draw	1-based value index; absent means 1

Return Value real(kind=real32)

a uniform draw in [0, 1)

private pure elemental function pf_random32_at_i64(seed, i, draw) result(r)

pf_random32_at for an integer(int64) stream index.

Arguments

Type	Intent	Optional		Name
integer(kind=int64),	intent(in)		::	seed	the stream family's seed
integer(kind=int64),	intent(in)		::	i	stream index; every value is valid
integer(kind=int64),	intent(in),	optional	::	draw	1-based value index; absent means 1

Return Value real(kind=real32)

a uniform draw in [0, 1)

public interface pf_random_bits_at

64 raw bits: value draw (default 1) of stream i under seed, as integer(int64).

Reads the same two words as pf_random_at, so pf_random_at is exactly this pattern's top 53 bits scaled into [0, 1). i is integer(int32) or integer(int64); seed and draw are integer(int64). Every one of the 2**64 patterns is possible.

private pure elemental function pf_random_bits_at_i32(seed, i, draw) result(r)

pf_random_bits_at for an integer(int32) stream index.

Arguments

Type	Intent	Optional		Name
integer(kind=int64),	intent(in)		::	seed	the stream family's seed
integer(kind=int32),	intent(in)		::	i	stream index; sign-extends, so any value is valid
integer(kind=int64),	intent(in),	optional	::	draw	1-based value index; absent means 1

Return Value integer(kind=int64)

64 raw bits

private pure elemental function pf_random_bits_at_i64(seed, i, draw) result(r)

pf_random_bits_at for an integer(int64) stream index.

Arguments

Type	Intent	Optional		Name
integer(kind=int64),	intent(in)		::	seed	the stream family's seed
integer(kind=int64),	intent(in)		::	i	stream index; every value is valid
integer(kind=int64),	intent(in),	optional	::	draw	1-based value index; absent means 1

Return Value integer(kind=int64)

64 raw bits

public interface pf_random_int_at

A uniform integer in [lo, hi], exactly unbiased: value draw (default 1) of stream i.

i, lo and hi share one kind -- integer(int32) or integer(int64) -- and the result follows it; seed and draw are always integer(int64). Swaps internally when lo > hi, so the function is total. Every width is exact, including the widest: no modulo bias, at any range, on either side of the fork.

private pure elemental function pf_random_int_at_i32(seed, i, lo, hi, draw) result(r)

pf_random_int_at for integer(int32) stream index and bounds.

The result is inside [min(lo,hi), max(lo,hi)] by construction, so narrowing the int64 worker's answer back to int32 is exact and cannot overflow.

Arguments

Type	Intent	Optional		Name
integer(kind=int64),	intent(in)		::	seed	the stream family's seed
integer(kind=int32),	intent(in)		::	i	stream index; sign-extends, so any value is valid
integer(kind=int32),	intent(in)		::	lo	one end of the closed range
integer(kind=int32),	intent(in)		::	hi	the other end; `lo > hi` is swapped, not an error
integer(kind=int64),	intent(in),	optional	::	draw	1-based value index; absent means 1

Return Value integer(kind=int32)

a uniform integer in [min(lo,hi), max(lo,hi)]

private pure elemental function pf_random_int_at_i64(seed, i, lo, hi, draw) result(r)

pf_random_int_at for integer(int64) stream index and bounds.

Arguments

Type	Intent	Optional		Name
integer(kind=int64),	intent(in)		::	seed	the stream family's seed
integer(kind=int64),	intent(in)		::	i	stream index; every value is valid
integer(kind=int64),	intent(in)		::	lo	one end of the closed range
integer(kind=int64),	intent(in)		::	hi	the other end; `lo > hi` is swapped, not an error
integer(kind=int64),	intent(in),	optional	::	draw	1-based value index; absent means 1

Return Value integer(kind=int64)

a uniform integer in [min(lo,hi), max(lo,hi)]

public interface pf_random_fill_draws

Fills v with consecutive values of one stream, starting at draw (default 1).

v is a rank-1 real(real64) or real(real32) array, intent(out); i is integer(int32) or integer(int64); seed and draw are integer(int64). The values are exactly what the matching scalar draws would give at those positions, so a prefix is a prefix: v(1:3) filled alone equals the first three of v(1:6). A zero-sized v is a defined no-op.

Precondition on the draw axis: draw + size(v) - 1 must not exceed huge(int64). The last element's position has to be representable, because there is no value at a position that cannot be named -- a fill that runs past the end is asking for draws that do not exist, and it silently receives wrapped ones. Everything up to and including the boundary is exact: a fill whose final position is huge(int64) itself is correct, and is tested. The scalar entry points have no such limit, since every representable draw is a valid one.

An INTEGER v takes two further required arguments, lo and hi, which share v's kind -- call pf_random_fill_draws(seed, i, v, lo, hi [, draw]). Element k is then exactly pf_random_int_at(seed, i, lo, hi, draw+k-1), the same identity the real forms have with pf_random_at. The specifics are distinguishable on v's type alone, so the generic resolves without ambiguity, and lo > hi is swapped rather than refused, exactly as in the scalar draw.

The integer form amortises much as real64 does: an integer draw has stride 2, so consecutive draws pair up two to a block and the fill enciphers once for each pair. Measured on machine B (gfortran 15.2.1, -O3 -funroll-loops, 4M values, best of three alternating rounds against the committed /v1 build, with the real64 fill flat at 8.66-8.69 ns as a cross-build control): 25.21-25.32 ns per value before, 15.74-15.76 after, 1.61x. It was not always so -- see pf_random_int_at's own note on the stride change that made it possible.

private pure subroutine pf_random_fill_draws_r64_i32(seed, i, v, draw)

pf_random_fill_draws filling real64 from an integer(int32) stream index.

Arguments

Type	Intent	Optional		Name
integer(kind=int64),	intent(in)		::	seed	the stream family's seed
integer(kind=int32),	intent(in)		::	i	stream index; sign-extends, so any value is valid
real(kind=real64),	intent(out)		::	v(:)	filled with values `draw .. draw+size(v)-1`
integer(kind=int64),	intent(in),	optional	::	draw	1-based starting value index; absent means 1

private pure subroutine pf_random_fill_draws_r64_i64(seed, i, v, draw)

pf_random_fill_draws filling real64 from an integer(int64) stream index.

Arguments

Type	Intent	Optional		Name
integer(kind=int64),	intent(in)		::	seed	the stream family's seed
integer(kind=int64),	intent(in)		::	i	stream index; every value is valid
real(kind=real64),	intent(out)		::	v(:)	filled with values `draw .. draw+size(v)-1`
integer(kind=int64),	intent(in),	optional	::	draw	1-based starting value index; absent means 1

private pure subroutine pf_random_fill_draws_r32_i32(seed, i, v, draw)

pf_random_fill_draws filling real32 from an integer(int32) stream index.

Arguments

Type	Intent	Optional		Name
integer(kind=int64),	intent(in)		::	seed	the stream family's seed
integer(kind=int32),	intent(in)		::	i	stream index; sign-extends, so any value is valid
real(kind=real32),	intent(out)		::	v(:)	filled with values `draw .. draw+size(v)-1`
integer(kind=int64),	intent(in),	optional	::	draw	1-based starting value index; absent means 1

private pure subroutine pf_random_fill_draws_r32_i64(seed, i, v, draw)

pf_random_fill_draws filling real32 from an integer(int64) stream index.

Arguments

Type	Intent	Optional		Name
integer(kind=int64),	intent(in)		::	seed	the stream family's seed
integer(kind=int64),	intent(in)		::	i	stream index; every value is valid
real(kind=real32),	intent(out)		::	v(:)	filled with values `draw .. draw+size(v)-1`
integer(kind=int64),	intent(in),	optional	::	draw	1-based starting value index; absent means 1

private pure subroutine pf_random_fill_draws_i32_i32(seed, i, v, lo, hi, draw)

pf_random_fill_draws filling integer(int32) from an integer(int32) stream index.

Arguments

Type	Intent	Optional		Name
integer(kind=int64),	intent(in)		::	seed	the stream family's seed
integer(kind=int32),	intent(in)		::	i	stream index; sign-extends, so any value is valid
integer(kind=int32),	intent(out)		::	v(:)	filled with values `draw .. draw+size(v)-1`
integer(kind=int32),	intent(in)		::	lo	one end of the closed range
integer(kind=int32),	intent(in)		::	hi	the other end; `lo > hi` is swapped, not an error
integer(kind=int64),	intent(in),	optional	::	draw	1-based starting value index; absent means 1

private pure subroutine pf_random_fill_draws_i32_i64(seed, i, v, lo, hi, draw)

pf_random_fill_draws filling integer(int32) from an integer(int64) stream index.

Arguments

Type	Intent	Optional		Name
integer(kind=int64),	intent(in)		::	seed	the stream family's seed
integer(kind=int64),	intent(in)		::	i	stream index; every value is valid
integer(kind=int32),	intent(out)		::	v(:)	filled with values `draw .. draw+size(v)-1`
integer(kind=int32),	intent(in)		::	lo	one end of the closed range
integer(kind=int32),	intent(in)		::	hi	the other end; `lo > hi` is swapped, not an error
integer(kind=int64),	intent(in),	optional	::	draw	1-based starting value index; absent means 1

private pure subroutine pf_random_fill_draws_i64_i32(seed, i, v, lo, hi, draw)

pf_random_fill_draws filling integer(int64) from an integer(int32) stream index.

Arguments

Type	Intent	Optional		Name
integer(kind=int64),	intent(in)		::	seed	the stream family's seed
integer(kind=int32),	intent(in)		::	i	stream index; sign-extends, so any value is valid
integer(kind=int64),	intent(out)		::	v(:)	filled with values `draw .. draw+size(v)-1`
integer(kind=int64),	intent(in)		::	lo	one end of the closed range
integer(kind=int64),	intent(in)		::	hi	the other end; `lo > hi` is swapped, not an error
integer(kind=int64),	intent(in),	optional	::	draw	1-based starting value index; absent means 1

private pure subroutine pf_random_fill_draws_i64_i64(seed, i, v, lo, hi, draw)

pf_random_fill_draws filling integer(int64) from an integer(int64) stream index.

Arguments

Type	Intent	Optional		Name
integer(kind=int64),	intent(in)		::	seed	the stream family's seed
integer(kind=int64),	intent(in)		::	i	stream index; every value is valid
integer(kind=int64),	intent(out)		::	v(:)	filled with values `draw .. draw+size(v)-1`
integer(kind=int64),	intent(in)		::	lo	one end of the closed range
integer(kind=int64),	intent(in)		::	hi	the other end; `lo > hi` is swapped, not an error
integer(kind=int64),	intent(in),	optional	::	draw	1-based starting value index; absent means 1

public interface pf_random_fill_streams

Fills v with ONE draw of each of size(v) consecutive streams, starting at stream i0.

The other axis. Where pf_random_fill_draws fixes the stream and walks the draws, this fixes the draw and walks the streams -- so it is the bulk form of the loop this module's own documentation opens with, x(i) = pf_random_at(seed, i). Element k is exactly pf_random_at(seed, i0 + k - 1 [, draw]), so the two forms are interchangeable and a prefix is a prefix.

v is a rank-1 real(real64) or real(real32) array, intent(out); i0 is integer(int32) or integer(int64); seed and draw are integer(int64), and draw (default 1) is which draw of every one of those streams to take. A zero-sized v is a defined no-op.

Precondition, on the STREAM axis this time: i0 + size(v) - 1 must not exceed huge(int64). Same reasoning as the draw-axis fill's own precondition -- the last element's stream has to be nameable. Note i0 may be negative, and usually the whole range is nowhere near the boundary.

It cannot be as cheap per value as pf_random_fill_draws, and that is the contract rather than the implementation. Each stream needs its own block, and one real64 value consumes two of that block's four words -- the other two belong to draw 2 of the same stream, which this call is not asking for. A real32 value consumes one of four. So where the draw-axis fill amortises one enciphering over two (or four) values, this one enciphers per value and wins only by removing the per-element call. Measured on machine B against the scalar loop it replaces: 1.28x on gfortran and 1.48x on ifx for real64, 1.47x and 1.60x for real32. When several values per stream are wanted, pf_random_fill_draws remains much the cheaper shape.

An INTEGER v takes two further required arguments, lo and hi, which share v's kind -- call pf_random_fill_streams(seed, i0, v, lo, hi [, draw]). Element k is exactly pf_random_int_at(seed, i0+k-1, lo, hi [, draw]). This axis was already one enciphering per value for the real kinds, so unlike the draw-axis integer form there is nothing given up here at all: it is the same work with the per-element call removed.

private pure subroutine pf_random_fill_streams_r64_i32(seed, i0, v, draw)

pf_random_fill_streams filling real64 from an integer(int32) first stream index.

Arguments

Type	Intent	Optional		Name
integer(kind=int64),	intent(in)		::	seed	the stream family's seed
integer(kind=int32),	intent(in)		::	i0	first stream index; sign-extends, so any value is valid
real(kind=real64),	intent(out)		::	v(:)	filled from streams `i0 .. i0+size(v)-1`
integer(kind=int64),	intent(in),	optional	::	draw	which draw of each stream; absent means 1

private pure subroutine pf_random_fill_streams_r64_i64(seed, i0, v, draw)

pf_random_fill_streams filling real64 from an integer(int64) first stream index.

Arguments

Type	Intent	Optional		Name
integer(kind=int64),	intent(in)		::	seed	the stream family's seed
integer(kind=int64),	intent(in)		::	i0	first stream index; every value is valid
real(kind=real64),	intent(out)		::	v(:)	filled from streams `i0 .. i0+size(v)-1`
integer(kind=int64),	intent(in),	optional	::	draw	which draw of each stream; absent means 1

private pure subroutine pf_random_fill_streams_r32_i32(seed, i0, v, draw)

pf_random_fill_streams filling real32 from an integer(int32) first stream index.

Arguments

Type	Intent	Optional		Name
integer(kind=int64),	intent(in)		::	seed	the stream family's seed
integer(kind=int32),	intent(in)		::	i0	first stream index; sign-extends, so any value is valid
real(kind=real32),	intent(out)		::	v(:)	filled from streams `i0 .. i0+size(v)-1`
integer(kind=int64),	intent(in),	optional	::	draw	which draw of each stream; absent means 1

private pure subroutine pf_random_fill_streams_r32_i64(seed, i0, v, draw)

pf_random_fill_streams filling real32 from an integer(int64) first stream index.

Arguments

Type	Intent	Optional		Name
integer(kind=int64),	intent(in)		::	seed	the stream family's seed
integer(kind=int64),	intent(in)		::	i0	first stream index; every value is valid
real(kind=real32),	intent(out)		::	v(:)	filled from streams `i0 .. i0+size(v)-1`
integer(kind=int64),	intent(in),	optional	::	draw	which draw of each stream; absent means 1

private pure subroutine pf_random_fill_streams_i32_i32(seed, i0, v, lo, hi, draw)

pf_random_fill_streams filling integer(int32) from an integer(int32) first stream index.

Arguments

Type	Intent	Optional		Name
integer(kind=int64),	intent(in)		::	seed	the stream family's seed
integer(kind=int32),	intent(in)		::	i0	first stream index; sign-extends, so any value is valid
integer(kind=int32),	intent(out)		::	v(:)	filled from streams `i0 .. i0+size(v)-1`
integer(kind=int32),	intent(in)		::	lo	one end of the closed range
integer(kind=int32),	intent(in)		::	hi	the other end; `lo > hi` is swapped, not an error
integer(kind=int64),	intent(in),	optional	::	draw	which draw of each stream; absent means 1

private pure subroutine pf_random_fill_streams_i32_i64(seed, i0, v, lo, hi, draw)

pf_random_fill_streams filling integer(int32) from an integer(int64) first stream index.

Arguments

Type	Intent	Optional		Name
integer(kind=int64),	intent(in)		::	seed	the stream family's seed
integer(kind=int64),	intent(in)		::	i0	first stream index; every value is valid
integer(kind=int32),	intent(out)		::	v(:)	filled from streams `i0 .. i0+size(v)-1`
integer(kind=int32),	intent(in)		::	lo	one end of the closed range
integer(kind=int32),	intent(in)		::	hi	the other end; `lo > hi` is swapped, not an error
integer(kind=int64),	intent(in),	optional	::	draw	which draw of each stream; absent means 1

private pure subroutine pf_random_fill_streams_i64_i32(seed, i0, v, lo, hi, draw)

pf_random_fill_streams filling integer(int64) from an integer(int32) first stream index.

Arguments

Type	Intent	Optional		Name
integer(kind=int64),	intent(in)		::	seed	the stream family's seed
integer(kind=int32),	intent(in)		::	i0	first stream index; sign-extends, so any value is valid
integer(kind=int64),	intent(out)		::	v(:)	filled from streams `i0 .. i0+size(v)-1`
integer(kind=int64),	intent(in)		::	lo	one end of the closed range
integer(kind=int64),	intent(in)		::	hi	the other end; `lo > hi` is swapped, not an error
integer(kind=int64),	intent(in),	optional	::	draw	which draw of each stream; absent means 1

private pure subroutine pf_random_fill_streams_i64_i64(seed, i0, v, lo, hi, draw)

pf_random_fill_streams filling integer(int64) from an integer(int64) first stream index.

Arguments

Type	Intent	Optional		Name
integer(kind=int64),	intent(in)		::	seed	the stream family's seed
integer(kind=int64),	intent(in)		::	i0	first stream index; every value is valid
integer(kind=int64),	intent(out)		::	v(:)	filled from streams `i0 .. i0+size(v)-1`
integer(kind=int64),	intent(in)		::	lo	one end of the closed range
integer(kind=int64),	intent(in)		::	hi	the other end; `lo > hi` is swapped, not an error
integer(kind=int64),	intent(in),	optional	::	draw	which draw of each stream; absent means 1

public interface pf_random_exp_at

One Exp(1) draw: value draw (default 1) of stream i under seed.

i is integer(int32) or integer(int64); seed and draw are integer(int64). The result is in [0, 36.7368] -- -log of the smallest uniform this generator can produce.

The mapping is -log(1 - u) where u is exactly pf_random_at(seed, i, draw), so the draw is bit-identical for a given libm: the logarithm is the intrinsic one, and libm is not part of any contract this project controls. pf_random_exp_portable_at is the same value computed through a frozen logarithm, identical on every platform -- and about 3x the cost, which is why this is the default rather than that one.

1 - u, not u, and that is contract. pf_random_at can return exactly 0 and can never return 1, so -log(u) would be infinite once in 2**53 draws while -log(1 - u) is finite always. The cost is that the draw can be exactly 0, which is correct and harmless.

Two words, the same two pf_random_at reads at that coordinate, so tier 0, tier 1 and the bulk fill are the same value computed the same way -- which the normal's forms deliberately are not, its consumption being variable (pf_random_normal_at says why).

"The same way" is not quite "bit for bit", and the reason is the libm promise showing its teeth inside one program. A compiler may serve log from a VECTOR libm inside the bulk fill's loop and a scalar one in this elemental call, and the two variants need not agree in the last bit. Measured on machine B: gfortran gives identical values on all three tiers, while ifx at its default -fp-model=fast differs on 21 of 64 by at most 2 ulp. pf_random_exp_portable_at has no such exposure -- its logarithm is this library's own and there is no vector variant to substitute -- and its three tiers agree exactly everywhere.

What is exact on every tier and every compiler is CHUNK INVARIANCE, which is the property that matters: a fill split at any boundary, in any order, on any number of threads gives the same values as one whole fill. That is asserted directly, and it holds because a given tier uses one code path for every element regardless of how many there are.

private pure elemental function pf_random_exp_at_i32(seed, i, draw) result(r)

pf_random_exp_at for an integer(int32) stream index.

Arguments

Type	Intent	Optional		Name
integer(kind=int64),	intent(in)		::	seed	the stream family's seed
integer(kind=int32),	intent(in)		::	i	stream index; sign-extends, so any value is valid
integer(kind=int64),	intent(in),	optional	::	draw	1-based value index; absent means 1

Return Value real(kind=real64)

an Exp(1) draw in [0, 36.7368]

private pure elemental function pf_random_exp_at_i64(seed, i, draw) result(r)

pf_random_exp_at for an integer(int64) stream index.

Arguments

Type	Intent	Optional		Name
integer(kind=int64),	intent(in)		::	seed	the stream family's seed
integer(kind=int64),	intent(in)		::	i	stream index; every value is valid
integer(kind=int64),	intent(in),	optional	::	draw	1-based value index; absent means 1

Return Value real(kind=real64)

an Exp(1) draw in [0, 36.7368]

public interface pf_random_exp_portable_at

pf_random_exp_at's value computed through a frozen logarithm: identical on every platform, compiler and flag set, not merely for a given libm.

Same arguments, same range, same two words, same 1 - u convention. The only difference is which logarithm: parquet_expkey's transform, built from IEEE + - * / with rounding barriers no compiler may reorder, rather than libm's. The two agree to about 2 ulp, which is eleven orders of magnitude below the Monte Carlo error of anything that could measure the difference -- so this is a choice about reproducibility, never about accuracy.

It costs about 3x, measured on machine B (gfortran 15.2.1, -O3 -funroll-loops): the bulk fill 39.4 ns per value against 11.2, the scalar draw 81.3 against 27.0. The frozen transform is twelve barriered Horner steps, each a store and a reload, so it neither vectorises nor pipelines while a libm log does both. Reach for it when a stored result must reproduce across machines; take the default otherwise.

Not pure, because the frozen transform's rounding barrier is a volatile local and a pure procedure may not have one. Elemental use still works, but a caller's own pure procedure and a do concurrent body cannot reach it -- pf_random_exp_at can, and pf_random_fill_exp_portable works outside the construct.

private impure elemental function pf_random_exp_portable_at_i32(seed, i, draw) result(r)

pf_random_exp_portable_at for an integer(int32) stream index.

impure because exp_key is: the frozen transform's rounding barrier is a volatile local, and a pure procedure may not have one. Elemental use is unaffected.

Arguments

Type	Intent	Optional		Name
integer(kind=int64),	intent(in)		::	seed	the stream family's seed
integer(kind=int32),	intent(in)		::	i	stream index; sign-extends, so any value is valid
integer(kind=int64),	intent(in),	optional	::	draw	1-based value index; absent means 1

Return Value real(kind=real64)

an Exp(1) draw in [0, 36.7368]

private impure elemental function pf_random_exp_portable_at_i64(seed, i, draw) result(r)

pf_random_exp_portable_at for an integer(int64) stream index.

Arguments

Type	Intent	Optional		Name
integer(kind=int64),	intent(in)		::	seed	the stream family's seed
integer(kind=int64),	intent(in)		::	i	stream index; every value is valid
integer(kind=int64),	intent(in),	optional	::	draw	1-based value index; absent means 1

Return Value real(kind=real64)

an Exp(1) draw in [0, 36.7368]

public interface pf_random_fill_exp

Fills v with consecutive Exp(1) draws of one stream, starting at draw (default 1).

v is a rank-1 real(real64) array, intent(out); i is integer(int32) or integer(int64); seed and draw are integer(int64). Element k is exactly pf_random_exp_at(seed, i, draw+k-1) -- to within the libm caveat that entry describes, and exactly as far as chunking is concerned: a prefix is a prefix and a fill split at any boundary agrees with a whole one, on every compiler. A zero-sized v is a defined no-op.

Precondition on the draw axis: draw + size(v) - 1 must not exceed huge(int64) -- the same bound, for the same reason, as pf_random_fill_draws.

No threads=, deliberately. Every element is a pure function of its coordinates, so the caller wraps their own !$omp parallel do around any chunking they like and gets the same answer at any thread count. An internal thread count would be less flexible and no faster.

private pure subroutine pf_random_fill_exp_i32(seed, i, v, draw)

pf_random_fill_exp from an integer(int32) stream index.

Arguments

Type	Intent	Optional		Name
integer(kind=int64),	intent(in)		::	seed	the stream family's seed
integer(kind=int32),	intent(in)		::	i	stream index; sign-extends, so any value is valid
real(kind=real64),	intent(out)		::	v(:)	filled with values `draw .. draw+size(v)-1`
integer(kind=int64),	intent(in),	optional	::	draw	1-based starting value index; absent means 1

private pure subroutine pf_random_fill_exp_i64(seed, i, v, draw)

pf_random_fill_exp from an integer(int64) stream index.

Arguments

Type	Intent	Optional		Name
integer(kind=int64),	intent(in)		::	seed	the stream family's seed
integer(kind=int64),	intent(in)		::	i	stream index; every value is valid
real(kind=real64),	intent(out)		::	v(:)	filled with values `draw .. draw+size(v)-1`
integer(kind=int64),	intent(in),	optional	::	draw	1-based starting value index; absent means 1

public interface pf_random_fill_exp_portable

pf_random_fill_exp through the frozen logarithm; element k is exactly pf_random_exp_portable_at(seed, i, draw+k-1).

Same arguments and same preconditions as pf_random_fill_exp, and the same reasons for taking no threads=. Not pure, for the reason pf_random_exp_portable_at gives.

private subroutine pf_random_fill_exp_portable_i32(seed, i, v, draw)

pf_random_fill_exp_portable from an integer(int32) stream index.

Arguments

Type	Intent	Optional		Name
integer(kind=int64),	intent(in)		::	seed	the stream family's seed
integer(kind=int32),	intent(in)		::	i	stream index; sign-extends, so any value is valid
real(kind=real64),	intent(out)		::	v(:)	filled with values `draw .. draw+size(v)-1`
integer(kind=int64),	intent(in),	optional	::	draw	1-based starting value index; absent means 1

private subroutine pf_random_fill_exp_portable_i64(seed, i, v, draw)

pf_random_fill_exp_portable from an integer(int64) stream index.

Arguments

Type	Intent	Optional		Name
integer(kind=int64),	intent(in)		::	seed	the stream family's seed
integer(kind=int64),	intent(in)		::	i	stream index; every value is valid
real(kind=real64),	intent(out)		::	v(:)	filled with values `draw .. draw+size(v)-1`
integer(kind=int64),	intent(in),	optional	::	draw	1-based starting value index; absent means 1

public interface pf_random_normal_at

One standard normal draw: value draw (default 1) of stream i under seed.

i is integer(int32) or integer(int64); seed and draw are integer(int64). Mean 0, variance 1, and the whole real line is reachable up to what the tail algorithm can produce.

Ziggurat over 256 equal-area layers (parquet_ziggurat), which accepts 98.5 % of draws after one 64-bit read and falls back to a wedge test or a tail walk otherwise. Bit-identical for a given libm: the wedge test needs exp and the tail needs log. pf_random_normal_portable_at is the same distribution through a frozen logarithm, identical on every platform.

This does NOT equal %normal at the same coordinate, and that is deliberate. A rejection algorithm consumes a number of words that depends on the values it drew, so no caller can say where value k of a stream walk begins without having drawn the preceding k-1 -- which would make a chunked or threaded fill impossible. The coordinate-addressed forms sidestep that by giving each (i, draw) its own derived sub-stream, so this value is a pure function of its coordinates and a bulk fill splits anywhere. The stream walk cannot do that, because being a walk is its purpose. See pf_normal_algorithm, and pf_random_exp_at for the one distribution where the two DO agree.

Independent of pf_random_normal_portable_at at the same coordinate, by construction: the two derive their sub-streams through different labels, so they are two normals rather than two functions of one uniform. Same for either one against its own stream walk.

private pure elemental function pf_random_normal_at_i32(seed, i, draw) result(r)

pf_random_normal_at for an integer(int32) stream index.

Arguments

Type	Intent	Optional		Name
integer(kind=int64),	intent(in)		::	seed	the stream family's seed
integer(kind=int32),	intent(in)		::	i	stream index; sign-extends, so any value is valid
integer(kind=int64),	intent(in),	optional	::	draw	1-based value index; absent means 1

Return Value real(kind=real64)

a standard normal draw

private pure elemental function pf_random_normal_at_i64(seed, i, draw) result(r)

pf_random_normal_at for an integer(int64) stream index.

Arguments

Type	Intent	Optional		Name
integer(kind=int64),	intent(in)		::	seed	the stream family's seed
integer(kind=int64),	intent(in)		::	i	stream index; every value is valid
integer(kind=int64),	intent(in),	optional	::	draw	1-based value index; absent means 1

Return Value real(kind=real64)

a standard normal draw

public interface pf_random_normal_portable_at

A standard normal that is identical on every platform, compiler and flag set, not merely for a given libm.

Same arguments and same distribution as pf_random_normal_at; a different algorithm, and so a different value at the same coordinate. Marsaglia's polar method: draw a point in the square, reject it unless it lands in the unit disc, and map the survivor through sqrt(-2*log(s)/s). Every step is IEEE + - * / and sqrt -- both correctly rounded by the standard -- over parquet_expkey's frozen logarithm, with rounding barriers on the two squares and on the argument of the sqrt so that no compiler may fuse or regroup them.

It costs more than the Ziggurat, in two ways. It rejects 21.5 % of candidate pairs rather than 1.5 % of single draws, it discards the second variate the method produces (a stream that kept it would have hidden state, so %rewind would stop being exact), and its logarithm is about 3x libm's. Reach for it when a stored result must reproduce across machines; take pf_random_normal_at otherwise.

private impure elemental function pf_random_normal_portable_at_i32(seed, i, draw) result(r)

pf_random_normal_portable_at for an integer(int32) stream index.

impure because exp_key is: the frozen transform's rounding barrier is a volatile local, and a pure procedure may not have one. Elemental use is unaffected.

Arguments

Type	Intent	Optional		Name
integer(kind=int64),	intent(in)		::	seed	the stream family's seed
integer(kind=int32),	intent(in)		::	i	stream index; sign-extends, so any value is valid
integer(kind=int64),	intent(in),	optional	::	draw	1-based value index; absent means 1

Return Value real(kind=real64)

a standard normal draw

private impure elemental function pf_random_normal_portable_at_i64(seed, i, draw) result(r)

pf_random_normal_portable_at for an integer(int64) stream index.

Arguments

Type	Intent	Optional		Name
integer(kind=int64),	intent(in)		::	seed	the stream family's seed
integer(kind=int64),	intent(in)		::	i	stream index; every value is valid
integer(kind=int64),	intent(in),	optional	::	draw	1-based value index; absent means 1

Return Value real(kind=real64)

a standard normal draw

public interface pf_random_fill_normal

Fills v with consecutive standard normal draws of one stream, starting at draw.

v is a rank-1 real(real64) array, intent(out); i is integer(int32) or integer(int64); seed and draw are integer(int64). Element k is exactly pf_random_normal_at(seed, i, draw+k-1), so a prefix is a prefix and a chunked fill agrees with a whole one. A zero-sized v is a defined no-op.

This is why the coordinate-addressed realisation exists. Each element runs its own rejection loop in its own sub-stream, so the fill splits at any boundary and gives the same answer at any thread count -- which a stream walk cannot, at any speed. It takes no threads= for exactly that reason: the caller wraps !$omp parallel do around whatever chunking they like.

It does not agree with a loop of %normal. See pf_random_normal_at.

private pure subroutine pf_random_fill_normal_i32(seed, i, v, draw)

pf_random_fill_normal from an integer(int32) stream index.

Arguments

Type	Intent	Optional		Name
integer(kind=int64),	intent(in)		::	seed	the stream family's seed
integer(kind=int32),	intent(in)		::	i	stream index; sign-extends, so any value is valid
real(kind=real64),	intent(out)		::	v(:)	filled with values `draw .. draw+size(v)-1`
integer(kind=int64),	intent(in),	optional	::	draw	1-based starting value index; absent means 1

private pure subroutine pf_random_fill_normal_i64(seed, i, v, draw)

pf_random_fill_normal from an integer(int64) stream index.

Arguments

Type	Intent	Optional		Name
integer(kind=int64),	intent(in)		::	seed	the stream family's seed
integer(kind=int64),	intent(in)		::	i	stream index; every value is valid
real(kind=real64),	intent(out)		::	v(:)	filled with values `draw .. draw+size(v)-1`
integer(kind=int64),	intent(in),	optional	::	draw	1-based starting value index; absent means 1

public interface pf_random_fill_normal_portable

pf_random_fill_normal through the polar method and the frozen logarithm; element k is exactly pf_random_normal_portable_at(seed, i, draw+k-1).

Same arguments and same preconditions, and the same reasons for taking no threads=. Not pure, for the reason pf_random_normal_portable_at gives.

private subroutine pf_random_fill_normal_portable_i32(seed, i, v, draw)

pf_random_fill_normal_portable from an integer(int32) stream index.

Arguments

Type	Intent	Optional		Name
integer(kind=int64),	intent(in)		::	seed	the stream family's seed
integer(kind=int32),	intent(in)		::	i	stream index; sign-extends, so any value is valid
real(kind=real64),	intent(out)		::	v(:)	filled with values `draw .. draw+size(v)-1`
integer(kind=int64),	intent(in),	optional	::	draw	1-based starting value index; absent means 1

private subroutine pf_random_fill_normal_portable_i64(seed, i, v, draw)

pf_random_fill_normal_portable from an integer(int64) stream index.

Arguments

Type	Intent	Optional		Name
integer(kind=int64),	intent(in)		::	seed	the stream family's seed
integer(kind=int64),	intent(in)		::	i	stream index; every value is valid
real(kind=real64),	intent(out)		::	v(:)	filled with values `draw .. draw+size(v)-1`
integer(kind=int64),	intent(in),	optional	::	draw	1-based starting value index; absent means 1

public interface pf_random_key

Derives an independent seed from a seed and a label, so one seed can fan out into families.

label is integer(int32) or integer(int64); seed and the result are integer(int64). A derived key IS a seed, so derivations compose by nesting.

private pure elemental function pf_random_key_i32(seed, label) result(r)

pf_random_key for an integer(int32) label.

Arguments

Type	Intent	Optional	Attributes		Name
integer(kind=int64),	intent(in)			::	seed	the seed to derive from
integer(kind=int32),	intent(in)			::	label	which derived family; sign-extends

Return Value integer(kind=int64)

an independent seed

private pure elemental function pf_random_key_i64(seed, label) result(r)

pf_random_key for an integer(int64) label.

Arguments

Type	Intent	Optional	Attributes		Name
integer(kind=int64),	intent(in)			::	seed	the seed to derive from
integer(kind=int64),	intent(in)			::	label	which derived family

Return Value integer(kind=int64)

an independent seed

Derived Types

type, public :: pf_random_stream

A walk along one stream: the same values tier 0 addresses, reached in sequence.

Type-Bound Procedures

generic, public :: seed => seed_base, seed_i32, seed_i64	(Re)seeds the stream to position 1. O(1), with no warm-up; `stream` defaults to 0.
procedure, public :: uniform => stream_uniform	Next `real64` in `[0, 1)`; costs 2 words.
procedure, public :: uniform32 => stream_uniform32	Next `real32` in `[0, 1)`; costs 1 word.
procedure, public :: bits => stream_bits	Next 64 raw bits; costs 2 words.
generic, public :: int_range => int_range_i32, int_range_i64	Next integer in `[lo, hi]`, exactly unbiased; costs one word pair, taken pair-aligned.
procedure, public :: exp => stream_exp	Next `Exp(1)`; costs one word pair, pair-aligned.
procedure, public :: exp_portable => stream_exp_portable	`%exp` through the frozen log; same cost.
procedure, public :: normal => stream_normal	Next standard normal; VARIABLE cost, pair-aligned.
procedure, public :: normal_portable => stream_normal_portable	`%normal` frozen; variable cost.
procedure, public :: gamma => stream_gamma	Next `Gamma(shape, 1)`; VARIABLE cost.
generic, public :: poisson => poisson_i32, poisson_i64	Next `Poisson(lambda)` count, into an `int32` or `int64`; VARIABLE cost.
generic, public :: fill => fill_arr_r64, fill_arr_r32, fill_arr_i32, fill_arr_i64	Fills `v` with the next `size(v)` values; an integer `v` also takes `lo` and `hi`.
generic, public :: rewind => rewind_base, rewind_i32, rewind_i64	Sets the position; with no argument, back to 1. Accepts any value `%position` gave.
procedure, public :: position => stream_position	Current 1-based word position.

Functions

public function pf_random_seed() result(s)

A fresh, nondeterministic seed, in [1, huge(int64)].

Arguments

None

Return Value integer(kind=int64)

a nondeterministic seed in [1, huge(int64)]

public pure function parquet_debug_random_uses_int128() result(r)

Reports which side of the route (e) fork was compiled. Test-only.

Arguments

None

Return Value logical

.true. if the 128-bit multiply was compiled

public pure function parquet_debug_random_uses_safe64() result(r)

Reports whether the OVERFLOW-FREE 64-bit arm was compiled. Test-only.

Arguments

None

Return Value logical

.true. if the overflow-free arm was compiled

Subroutines

public pure subroutine parquet_debug_gamma_path(rng, shape, r, path, squeeze_ok)

Runs the Philox block function directly, on raw counter and key words. Test-only.

Arguments

Type	Intent	Optional		Name
type(pf_random_stream),	intent(inout)		::	rng	the stream to advance, as `%gamma` would
real(kind=real64),	intent(in)		::	shape	the shape parameter; must be > 0
real(kind=real64),	intent(out)		::	r	the draw, equal to `%gamma`'s
integer(kind=int32),	intent(out)		::	path	1/2 squeeze/log; +2 when boosted
logical,	intent(out),	optional	::	squeeze_ok	`.false.` iff the squeeze accepted a candidate the full test would reject

public pure subroutine parquet_debug_poisson_path(rng, lambda, k, path, squeeze_ok)

Reports which algorithm and branch a %poisson draw took. Test-only; see parquet_debug_normal_path.

Arguments

Type	Intent	Optional		Name
type(pf_random_stream),	intent(inout)		::	rng	the stream to advance, as `%poisson` would
real(kind=real64),	intent(in)		::	lambda	the mean; must be >= 0 and finite
integer(kind=int64),	intent(out)		::	k	the count, equal to `%poisson`'s
integer(kind=int32),	intent(out)		::	path	1 Knuth, 2 PTRS fast, 3 PTRS log test
logical,	intent(out),	optional	::	squeeze_ok	`.false.` iff the fast region accepted a candidate the full test would reject

public subroutine parquet_debug_normal_path(seed, i, draw, portable, x, path, pairs)

Arguments

Type	Intent		Name
integer(kind=int64),	intent(in)	::	seed	the stream family's seed
integer(kind=int64),	intent(in)	::	i	stream index
integer(kind=int64),	intent(in)	::	draw	1-based value index
logical,	intent(in)	::	portable	`.true.` for the polar form, `.false.` Ziggurat
real(kind=real64),	intent(out)	::	x	the draw, equal to the ordinary call's
integer(kind=int32),	intent(out)	::	path	1 rectangle, 2 wedge, 3 tail
integer(kind=int64),	intent(out)	::	pairs	word pairs the rejection loop consumed

public pure subroutine parquet_debug_random_block(key, stream, index, w0, w1, w2, w3)

Arguments

Type	Intent		Name
integer(kind=int64),	intent(in)	::	key	64-bit key: key words 0 and 1, low half first
integer(kind=int64),	intent(in)	::	stream	splits into counter words 2 and 3
integer(kind=int64),	intent(in)	::	index	0-based block index: counter words 0 and 1
integer(kind=int64),	intent(out)	::	w0	output word `c0`
integer(kind=int64),	intent(out)	::	w1	output word `c1`
integer(kind=int64),	intent(out)	::	w2	output word `c2`
integer(kind=int64),	intent(out)	::	w3	output word `c3`

parquet_random Module

Contents

Variables

Interfaces

Derived Types

Functions

Subroutines

Uses

Variables

Interfaces

public interface pf_random_at

private pure elemental function pf_random_at_i32(seed, i, draw) result(r)

Arguments

Return Value real(kind=real64)

private pure elemental function pf_random_at_i64(seed, i, draw) result(r)

Arguments

Return Value real(kind=real64)

public interface pf_random32_at

private pure elemental function pf_random32_at_i32(seed, i, draw) result(r)

Arguments

Return Value real(kind=real32)

private pure elemental function pf_random32_at_i64(seed, i, draw) result(r)

Arguments

Return Value real(kind=real32)

public interface pf_random_bits_at

private pure elemental function pf_random_bits_at_i32(seed, i, draw) result(r)

Arguments

Return Value integer(kind=int64)

private pure elemental function pf_random_bits_at_i64(seed, i, draw) result(r)

Arguments

Return Value integer(kind=int64)

public interface pf_random_int_at

private pure elemental function pf_random_int_at_i32(seed, i, lo, hi, draw) result(r)

Arguments

Return Value integer(kind=int32)

private pure elemental function pf_random_int_at_i64(seed, i, lo, hi, draw) result(r)

Arguments

Return Value integer(kind=int64)

public interface pf_random_fill_draws

private pure subroutine pf_random_fill_draws_r64_i32(seed, i, v, draw)

Arguments

private pure subroutine pf_random_fill_draws_r64_i64(seed, i, v, draw)

Arguments

private pure subroutine pf_random_fill_draws_r32_i32(seed, i, v, draw)

Arguments

private pure subroutine pf_random_fill_draws_r32_i64(seed, i, v, draw)

Arguments

private pure subroutine pf_random_fill_draws_i32_i32(seed, i, v, lo, hi, draw)

Arguments

private pure subroutine pf_random_fill_draws_i32_i64(seed, i, v, lo, hi, draw)

Arguments

private pure subroutine pf_random_fill_draws_i64_i32(seed, i, v, lo, hi, draw)

Arguments

private pure subroutine pf_random_fill_draws_i64_i64(seed, i, v, lo, hi, draw)

Arguments

public interface pf_random_fill_streams

private pure subroutine pf_random_fill_streams_r64_i32(seed, i0, v, draw)

Arguments

private pure subroutine pf_random_fill_streams_r64_i64(seed, i0, v, draw)

Arguments

private pure subroutine pf_random_fill_streams_r32_i32(seed, i0, v, draw)

Arguments

private pure subroutine pf_random_fill_streams_r32_i64(seed, i0, v, draw)

Arguments

private pure subroutine pf_random_fill_streams_i32_i32(seed, i0, v, lo, hi, draw)

Arguments

private pure subroutine pf_random_fill_streams_i32_i64(seed, i0, v, lo, hi, draw)

Arguments

private pure subroutine pf_random_fill_streams_i64_i32(seed, i0, v, lo, hi, draw)

Arguments

private pure subroutine pf_random_fill_streams_i64_i64(seed, i0, v, lo, hi, draw)

Arguments

public interface pf_random_exp_at

private pure elemental function pf_random_exp_at_i32(seed, i, draw) result(r)

Arguments

Return Value real(kind=real64)

private pure elemental function pf_random_exp_at_i64(seed, i, draw) result(r)

Arguments

Return Value real(kind=real64)

public interface pf_random_exp_portable_at