pf_random_exp_at Interface

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.


Module Procedures

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 IntentOptional Attributes 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 IntentOptional Attributes 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]