pf_random_stream Derived Type

type, public :: pf_random_stream

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

Tier 0 answers "what is the value at this coordinate?". Some programs cannot ask that, because how many values they need is data-dependent -- a rejection sampler, a random walk, a resample of unknown length. This type carries the position so the caller does not have to, and hands out consecutive values of one stream.

It changes no value. A freshly seeded stream's k-th %uniform is exactly pf_random_at(seed, stream, k), its k-th %uniform32 exactly pf_random32_at(seed, stream, k). The stream is a different way to reach the same grid, never a second generator.

Reproducibility is per iteration, and that is the discipline to follow: seed at the top of each loop iteration from a run-invariant label, then draw as many values as that iteration needs.

type(pf_random_stream) :: rng
!$omp parallel do schedule(dynamic) private(...)
do i = 1, n
    call rng%seed(seed, i)              ! O(1), no warm-up
    do while (...)                      ! however many draws this iteration turns out to need
        call rng%uniform(x)
    end do
end do

What can never be reproducible is one long-lived stream consumed across the iterations of a dynamically scheduled loop -- the value an iteration receives then depends on how many draws ran before it, which depends on the schedule. That is a property of every stateful generator, not a limitation of this one; the answer to it is that %seed costs nothing.

Position is measured in 32-bit words, 1-based, and the word cost of each producer is contract -- %position and %rewind are denominated in it: %uniform 2, %uniform32 1, %bits 2, %int_range 2, %exp 2, %exp_portable 2. %int_range and the two %exp forms additionally start on a word PAIR boundary, advancing one word first if a %uniform32 has left the cursor odd. This is what keeps them equal to the pf_random_int_at/pf_random_exp_at at the same coordinate rather than re-reading a word a previous draw already used. %int_range cost 4 words plus up to 3 of alignment until pf_random_algorithm reached /v2, when the integer generic's stride became 2.

A producer whose cost is VARIABLE cannot be predicted, only observed. Every producer listed above has a fixed cost, so a caller can compute where the stream will be after a known sequence of draws. %normal and %normal_portable cannot: each runs a rejection loop and consumes a number of words that depends on the values it drew -- 2 words in 98.5 % of %normal's draws and more otherwise, a multiple of 4 for %normal_portable. %position stays exact for those, because it reports where the stream is; only predicting it in advance becomes impossible. Checkpoint and restart (%position then %rewind) therefore keep working for every producer, while arithmetic on positions stops.

A variable-cost producer's value is NOT the coordinate-addressed one at the same position, and that is the price of a bulk fill being splittable at all -- see pf_random_normal_at, which explains why no implementation can have both.

The type is plain scalars: no allocatable components, no FINAL, deliberately and permanently. gfortran does not reliably default-initialise an OpenMP private() copy of a finalizable type, and ifx segfaults on a block-local instance of a type with allocatable components inside a parallel region (feature_risks.md Risk-45). A type that is neither is safe in both shapes, which is what makes a per-thread instance usable at all. Adding either to this type would break every parallel use of it, on one compiler or the other.

It holds the block it last enciphered, keyed by that block's index. One enciphering carries two real64 values or four real32s, so keeping it is worth 1.70x (gfortran) / 1.78x (ifx) on real64 and 2.88x / 3.35x on real32 -- measured, and enough to take the stream from slower than a tier-0 loop to faster than one. Because the cache is keyed rather than consumed, it reaches nothing in the contract: %position still means a word index, and a stream saved and restored through %position alone is exact.

pf_random_fill_draws is still cheaper again (1.87x / 1.67x against this type), so bulk work whose length is known in advance belongs there, not in a loop over a stream.


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.

  • private pure subroutine stream_seed_base(self, seed)

    %seed with no stream index: stream 0.

    Arguments

    Type IntentOptional Attributes Name
    class(pf_random_stream), intent(inout) :: self

    the stream to reseed
    integer(kind=int64), intent(in) :: seed

    the stream family's seed

  • private pure subroutine stream_seed_i32(self, seed, stream)

    %seed with an integer(int32) stream index; sign-extends, so any value is valid.

    Arguments

    Type IntentOptional Attributes Name
    class(pf_random_stream), intent(inout) :: self

    the stream to reseed
    integer(kind=int64), intent(in) :: seed

    the stream family's seed
    integer(kind=int32), intent(in) :: stream

    which stream of that family

  • private pure subroutine stream_seed_i64(self, seed, stream)

    %seed with an integer(int64) stream index; every value is valid.

    Arguments

    Type IntentOptional Attributes Name
    class(pf_random_stream), intent(inout) :: self

    the stream to reseed
    integer(kind=int64), intent(in) :: seed

    the stream family's seed
    integer(kind=int64), intent(in) :: stream

    which stream of that family

procedure, public :: uniform => stream_uniform

Next real64 in [0, 1); costs 2 words.

  • private pure subroutine stream_uniform(self, x)

    The next real64 in [0, 1), advancing two words.

    Arguments

    Type IntentOptional Attributes Name
    class(pf_random_stream), intent(inout) :: self

    the stream to advance
    real(kind=real64), intent(out) :: x

    a uniform draw in [0, 1)

procedure, public :: uniform32 => stream_uniform32

Next real32 in [0, 1); costs 1 word.

  • private pure subroutine stream_uniform32(self, x)

    The next real32 in [0, 1), advancing one word.

    Arguments

    Type IntentOptional Attributes Name
    class(pf_random_stream), intent(inout) :: self

    the stream to advance
    real(kind=real32), intent(out) :: x

    a uniform draw in [0, 1)

procedure, public :: bits => stream_bits

Next 64 raw bits; costs 2 words.

  • private pure subroutine stream_bits(self, b)

    The next 64 raw bits, advancing two words -- the same two %uniform would have read.

    Arguments

    Type IntentOptional Attributes Name
    class(pf_random_stream), intent(inout) :: self

    the stream to advance
    integer(kind=int64), intent(out) :: b

    64 raw bits

generic, public :: int_range => int_range_i32, int_range_i64

Next integer in [lo, hi], exactly unbiased; costs one word pair, taken pair-aligned.

  • private pure subroutine stream_int_range_i32(self, lo, hi, r)

    %int_range for integer(int32) bounds and result.

    The result is inside the closed range by construction, so the narrowing is exact -- the same argument pf_random_int_at_i32 rests on.

    Arguments

    Type IntentOptional Attributes Name
    class(pf_random_stream), intent(inout) :: self

    the stream to advance
    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
    integer(kind=int32), intent(out) :: r

    a uniform integer in the closed range

  • private pure subroutine stream_int_range_i64(self, lo, hi, r)

    %int_range for integer(int64) bounds and result.

    Arguments

    Type IntentOptional Attributes Name
    class(pf_random_stream), intent(inout) :: self

    the stream to advance
    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
    integer(kind=int64), intent(out) :: r

    a uniform integer in the closed range

procedure, public :: exp => stream_exp

Next Exp(1); costs one word pair, pair-aligned.

  • private pure subroutine stream_exp(self, x)

    %exp: the next Exp(1) draw, taking one pair-aligned word pair.

    Pair-aligned like %int_range, not free-running like %uniform, and the promise that buys is unconditional: the value is exactly pf_random_exp_at(seed, stream, d) for the draw d it consumed, whatever the cursor was beforehand. %uniform agrees with its coordinate-addressed twin only while the cursor stays even, which it does unless a %uniform32 has been mixed in. The price is at most one wasted word after such a mix.

    Arguments

    Type IntentOptional Attributes Name
    class(pf_random_stream), intent(inout) :: self

    the stream to advance
    real(kind=real64), intent(out) :: x

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

procedure, public :: exp_portable => stream_exp_portable

%exp through the frozen log; same cost.

  • private subroutine stream_exp_portable(self, x)

    %exp_portable: the same draw through the frozen logarithm. Not pure, per exp_key.

    Arguments

    Type IntentOptional Attributes Name
    class(pf_random_stream), intent(inout) :: self

    the stream to advance
    real(kind=real64), intent(out) :: x

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

procedure, public :: normal => stream_normal

Next standard normal; VARIABLE cost, pair-aligned.

  • private pure subroutine stream_normal(self, x)

    %normal: the next Ziggurat normal, consuming as many pairs as its rejection loop needed.

    The one producer here that cannot advance before it reads, because how far to advance is what the read decides. The invariant every other producer keeps -- that an exhausted stream aborts having written and moved nothing -- is preserved anyway by computing into a local and assigning x only after advance_by has accepted the cost.

    Arguments

    Type IntentOptional Attributes Name
    class(pf_random_stream), intent(inout) :: self

    the stream to advance
    real(kind=real64), intent(out) :: x

    a standard normal draw

procedure, public :: normal_portable => stream_normal_portable

%normal frozen; variable cost.

  • private subroutine stream_normal_portable(self, x)

    %normal_portable: the next polar normal. Not pure, per exp_key. See stream_normal.

    Arguments

    Type IntentOptional Attributes Name
    class(pf_random_stream), intent(inout) :: self

    the stream to advance
    real(kind=real64), intent(out) :: x

    a standard normal draw

procedure, public :: gamma => stream_gamma

Next Gamma(shape, 1); VARIABLE cost.

  • private pure subroutine stream_gamma(self, shape, r)

    %gamma: the next Gamma(shape, 1) draw, by Marsaglia-Tsang rejection.

    Tier 1 only, deliberately -- there is no pf_random_gamma_at and no bulk fill. The sub-stream construction ยง3.4 of feature_random_phase3.md describes would make them possible, and the surface is already wide; if they are ever wanted, the construction and its independence obligations are the same as the normal's.

    shape is the Gamma shape parameter a, strictly positive; the scale is 1, so a caller wanting Gamma(a, theta) multiplies by theta. The mean is shape and the variance is shape.

    Two rejection loops, not one. The inner one redraws a normal until 1 + c*x > 0, which happens for about Phi(-3*sqrt(a - 1/3)) of draws and so essentially never above a = 1; the outer one is the squeeze-then-log acceptance, which accepts about 95 % of candidates on the squeeze alone. Consumption is therefore variable and unpredictable -- see pf_random_stream's note on what %position still guarantees.

    shape < 1 is supported through a boost, Gamma(a) = Gamma(a+1) * u**(1/a), at the cost of one extra uniform and a libm pow. Restricting the domain to shape >= 1 instead was considered and rejected (Q3-4): a partial procedure is a worse API than a per-libm promise, and Gamma is per-libm regardless. 1 - u rather than u in the boost, so the draw cannot come back exactly 0.

    Arguments

    Type IntentOptional Attributes Name
    class(pf_random_stream), intent(inout) :: self

    the stream to advance
    real(kind=real64), intent(in) :: shape

    the shape parameter; must be > 0
    real(kind=real64), intent(out) :: r

    a Gamma(shape, 1) draw, strictly positive

generic, public :: poisson => poisson_i32, poisson_i64

Next Poisson(lambda) count, into an int32 or int64; VARIABLE cost.

  • private pure subroutine stream_poisson_i32(self, lambda, k)

    %poisson into an integer(int32).

    Refuses a count that does not fit rather than narrowing it, which would silently return a plausible wrong number. Reaching that needs a lambda of order 2**31, at which point the int64 specific is what the caller wants.

    Arguments

    Type IntentOptional Attributes Name
    class(pf_random_stream), intent(inout) :: self

    the stream to advance
    real(kind=real64), intent(in) :: lambda

    the mean; must be >= 0 and finite
    integer(kind=int32), intent(out) :: k

    a Poisson(lambda) count

  • private pure subroutine stream_poisson_i64(self, lambda, k)

    %poisson into an integer(int64).

    Arguments

    Type IntentOptional Attributes Name
    class(pf_random_stream), intent(inout) :: self

    the stream to advance
    real(kind=real64), intent(in) :: lambda

    the mean; must be >= 0 and finite
    integer(kind=int64), intent(out) :: k

    a Poisson(lambda) count

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.

  • private pure subroutine stream_fill_r64(self, v)

    %fill for a real64 array.

    Routes to fill_r64 whenever the position is pair-aligned, because the bulk fills walk blocks rather than values and beat even the cached stream by 1.67-1.87x. The values are identical either way; only the number of encipherings differs.

    Arguments

    Type IntentOptional Attributes Name
    class(pf_random_stream), intent(inout) :: self

    the stream to advance
    real(kind=real64), intent(out) :: v(:)

    filled with the next size(v) values

  • private pure subroutine stream_fill_r32(self, v)

    %fill for a real32 array. Every position is aligned for it: one value is one word.

    Arguments

    Type IntentOptional Attributes Name
    class(pf_random_stream), intent(inout) :: self

    the stream to advance
    real(kind=real32), intent(out) :: v(:)

    filled with the next size(v) values

  • private pure subroutine stream_fill_i32(self, v, lo, hi)

    %fill for an integer(int32) array.

    Arguments

    Type IntentOptional Attributes Name
    class(pf_random_stream), intent(inout) :: self

    the stream to advance
    integer(kind=int32), intent(out) :: v(:)

    filled with the next size(v) values
    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

  • private pure subroutine stream_fill_i64(self, v, lo, hi)

    %fill for an integer(int64) array; aligns to a pair first, exactly as %int_range does.

    Arguments

    Type IntentOptional Attributes Name
    class(pf_random_stream), intent(inout) :: self

    the stream to advance
    integer(kind=int64), intent(out) :: v(:)

    filled with the next size(v) values
    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

generic, public :: rewind => rewind_base, rewind_i32, rewind_i64

Sets the position; with no argument, back to 1. Accepts any value %position gave.

  • private pure subroutine stream_rewind_base(self)

    %rewind with no argument: back to position 1.

    Arguments

    Type IntentOptional Attributes Name
    class(pf_random_stream), intent(inout) :: self

    the stream to reposition

  • private pure subroutine stream_rewind_i32(self, pos)

    %rewind to an integer(int32) position.

    Arguments

    Type IntentOptional Attributes Name
    class(pf_random_stream), intent(inout) :: self

    the stream to reposition
    integer(kind=int32), intent(in) :: pos

    1-based word position

  • private pure subroutine stream_rewind_i64(self, pos)

    %rewind to an integer(int64) position.

    Arguments

    Type IntentOptional Attributes Name
    class(pf_random_stream), intent(inout) :: self

    the stream to reposition
    integer(kind=int64), intent(in) :: pos

    1-based word position

procedure, public :: position => stream_position

Current 1-based word position.

  • private pure function stream_position(self) result(p)

    The current 1-based word position.

    Arguments

    Type IntentOptional Attributes Name
    class(pf_random_stream), intent(in) :: self

    the stream to query

    Return Value integer(kind=int64)

    1-based word position

Source Code

    type, public :: pf_random_stream
        private
        integer(int64) :: key = 0_int64             !! the stream family's seed
        integer(int64) :: stream = 0_int64          !! which stream of that family
        integer(int64) :: pos = 0_int64             !! 0-based word position; `%position` reports `pos+1`
        integer(int64) :: blk = -1_int64            !! block index held below, or -1 when none is
        integer(int64) :: c0 = 0_int64              !! held word 0
        integer(int64) :: c1 = 0_int64              !! held word 1
        integer(int64) :: c2 = 0_int64              !! held word 2
        integer(int64) :: c3 = 0_int64              !! held word 3
    contains
        procedure, private :: seed_base => stream_seed_base   !! `%seed` with no stream index
        procedure, private :: seed_i32 => stream_seed_i32     !! `%seed` with an `int32` stream index
        procedure, private :: seed_i64 => stream_seed_i64     !! `%seed` with an `int64` stream index
        !> (Re)seeds the stream to position 1. O(1), with no warm-up; `stream` defaults to 0.
        generic :: seed => seed_base, seed_i32, seed_i64
        procedure :: uniform => stream_uniform      !! Next `real64` in `[0, 1)`; costs 2 words.
        procedure :: uniform32 => stream_uniform32  !! Next `real32` in `[0, 1)`; costs 1 word.
        procedure :: bits => stream_bits            !! Next 64 raw bits; costs 2 words.
        procedure, private :: int_range_i32 => stream_int_range_i32  !! `%int_range`, `int32`
        procedure, private :: int_range_i64 => stream_int_range_i64  !! `%int_range`, `int64`
        !> Next integer in `[lo, hi]`, exactly unbiased; costs one word pair, taken pair-aligned.
        generic :: int_range => int_range_i32, int_range_i64
        procedure :: exp => stream_exp              !! Next `Exp(1)`; costs one word pair, pair-aligned.
        procedure :: exp_portable => stream_exp_portable  !! `%exp` through the frozen log; same cost.
        procedure :: normal => stream_normal        !! Next standard normal; VARIABLE cost, pair-aligned.
        procedure :: normal_portable => stream_normal_portable  !! `%normal` frozen; variable cost.
        procedure :: gamma => stream_gamma          !! Next `Gamma(shape, 1)`; VARIABLE cost.
        procedure, private :: poisson_i32 => stream_poisson_i32 !! `%poisson`, `int32` result
        procedure, private :: poisson_i64 => stream_poisson_i64 !! `%poisson`, `int64` result
        !> Next `Poisson(lambda)` count, into an `int32` or `int64`; VARIABLE cost.
        generic :: poisson => poisson_i32, poisson_i64
        procedure, private :: fill_arr_r64 => stream_fill_r64        !! `%fill`, `real64`
        procedure, private :: fill_arr_r32 => stream_fill_r32        !! `%fill`, `real32`
        procedure, private :: fill_arr_i32 => stream_fill_i32        !! `%fill`, `int32`
        procedure, private :: fill_arr_i64 => stream_fill_i64        !! `%fill`, `int64`
        !> Fills `v` with the next `size(v)` values; an integer `v` also takes `lo` and `hi`.
        generic :: fill => fill_arr_r64, fill_arr_r32, fill_arr_i32, fill_arr_i64
        procedure, private :: rewind_base => stream_rewind_base      !! `%rewind` to position 1
        procedure, private :: rewind_i32 => stream_rewind_i32        !! `%rewind`, `int32`
        procedure, private :: rewind_i64 => stream_rewind_i64        !! `%rewind`, `int64`
        !> Sets the position; with no argument, back to 1. Accepts any value `%position` gave.
        generic :: rewind => rewind_base, rewind_i32, rewind_i64
        procedure :: position => stream_position    !! Current 1-based word position.
    end type pf_random_stream