pic_rng Module

Portable, reproducible pseudo-random number generators.

Two engines are provided:

  • splitmix64_t - SplitMix64 (Steele, Lea and Flood, 2014; reference implementation by Sebastiano Vigna, public domain). A 64-bit counter plus a strong finalising mix. Tiny, fast, and ideal for seeding other generators.
  • pcg32_t - PCG-XSH-RR 64/32 (O’Neill, 2014). A 64-bit LCG whose output is a permuted 32-bit word. Supports 2**63 distinct streams through its increment.

Both types expose seed, next and next_u64, and both work with the module-level generics next_real_dp, next_below and stream_for.

Portability notes (the whole reason this module exists)

  1. No BOZ literals anywhere. The 64-bit constants of SplitMix64 have their top bit set, so z'9E3779B97F4A7C15' is not representable as a positive integer(int64). Compilers disagree on BOZ literals appearing in initialisation expressions (gfortran, ifx and nvfortran all differ in what they accept and what they warn about). Every constant below is therefore written as an ordinary signed decimal literal equal to the two’s-complement interpretation of the intended bit pattern, with the hexadecimal value given in a comment.
  2. No reliance on signed overflow. Fortran has no unsigned integer type, and signed overflow is not defined by the standard: a build with -ftrapv or -fsanitize=signed-integer-overflow may abort on it. All modular arithmetic goes through u64_add and u64_mul in pic_uint64, which work on limbs so that no intermediate ever leaves the signed 64-bit range, and reassemble with ior/ishft, which are bit-model intrinsics and exactly defined for every bit pattern, sign bit included. See that module for the limb widths and the bounds; they are stated once, where they are implemented.
  3. No transfer between reals and integers. nvfortran and LFortran do not agree on the result of bit-casting a double to an integer. Reals are built arithmetically instead (see next_real_dp).
  4. No right shift by a negative count and no F2008 shift intrinsics. ishft(x, -n) is standard but its interaction with the sign bit is a recurring source of confusion, and shiftr/shifta are Fortran 2008 and not uniformly available on the older compilers PIC targets. Every logical right shift goes through u64_shr in pic_uint64, which is built from ibits: Fortran 90, zero-filling by definition, and correct for negative x. ishft is used only with non-negative (left-shift) counts.

Uses

  • module~~pic_rng~~UsesGraph module~pic_rng pic_rng module~pic_error pic_error module~pic_rng->module~pic_error module~pic_types pic_types module~pic_rng->module~pic_types module~pic_uint64 pic_uint64 module~pic_rng->module~pic_uint64 module~pic_error->module~pic_types module~pic_global_definitions pic_global_definitions module~pic_error->module~pic_global_definitions iso_c_binding iso_c_binding module~pic_types->iso_c_binding iso_fortran_env iso_fortran_env module~pic_types->iso_fortran_env module~pic_uint64->module~pic_types module~pic_global_definitions->module~pic_types

Used by

  • module~~pic_rng~~UsedByGraph module~pic_rng pic_rng module~pic_random_dist pic_random_dist module~pic_random_dist->module~pic_rng module~pic_random_dist_real pic_random_dist_real module~pic_random_dist_real->module~pic_rng

Variables

Type Visibility Attributes Name Initial
integer(kind=int64), private, parameter :: GOLDEN_GAMMA = -7046029254386353131_int64

SplitMix64 additive constant, bit pattern 0x9E3779B97F4A7C15.
integer(kind=int64), private, parameter :: MASK32 = 4294967295_int64

Low 32-bit mask (2**32 - 1).
integer(kind=int64), private, parameter :: MIX_C1 = -4658895280553007687_int64

SplitMix64 first multiplier, bit pattern 0xBF58476D1CE4E5B9.
integer(kind=int64), private, parameter :: MIX_C2 = -7723592293110705685_int64

SplitMix64 second multiplier, bit pattern 0x94D049BB133111EB.
integer(kind=int64), private, parameter :: PCG_DEFAULT_INCREMENT = 1442695040888963407_int64

Default PCG stream increment, bit pattern 0x14057B7EF767814F.
integer(kind=int64), private, parameter :: PCG_MULTIPLIER = 6364136223846793005_int64

PCG 64-bit LCG multiplier, bit pattern 0x5851F42D4C957F2D.
integer(kind=int64), private, parameter :: TWO_TO_31 = 2147483648_int64

2**31, the first value that does not fit in a signed 32-bit integer.
integer(kind=int64), private, parameter :: TWO_TO_32 = 4294967296_int64

2**32, used to fold an unsigned 32-bit value into int32.
real(kind=dp), private, parameter :: TWO_TO_MINUS_53 = 1.0_dp/9007199254740992.0_dp

2**-53, exactly representable in IEEE binary64.

Interfaces

public interface next_below

Draw a uniform integer in [0, n) by rejection sampling.

The naive mod(next(), n) is biased whenever n does not divide the generator’s period: the first 2**64 mod n residues each get one extra preimage, so for large n some outcomes are measurably more likely. Instead a mask of the smallest 2**k - 1 at least n - 1 is applied to a fresh 64-bit word and values at or above n are rejected and redrawn. Every accepted value has exactly one preimage per draw, so the result is exactly uniform. The acceptance probability is above 1/2, so the expected number of draws is below 2 and the loop terminates with probability one.

n == 1 returns 0 without consuming a draw. n <= 0 is a usage error: it returns 0 and, if error is present, sets ERROR_VALIDATION.

  • private function splitmix64_next_below(gen, n, error) result(r)

    Uniform integer in [0, n) from a SplitMix64 generator.

    Arguments

    Type IntentOptional Attributes Name
    type(splitmix64_t), intent(inout) :: gen
    integer(kind=default_int), intent(in) :: n
    type(error_t), intent(inout), optional :: error

    Return Value integer(kind=default_int)

  • private function pcg32_next_below(gen, n, error) result(r)

    Uniform integer in [0, n) from a PCG32 generator. Each attempt consumes two 32-bit outputs, so the stream stays word aligned no matter how many rejections occur.

    Arguments

    Type IntentOptional Attributes Name
    type(pcg32_t), intent(inout) :: gen
    integer(kind=default_int), intent(in) :: n
    type(error_t), intent(inout), optional :: error

    Return Value integer(kind=default_int)

public interface next_real_dp

Draw a uniform real(dp) in [0, 1).

Construction, identical for both engines and on every compiler: take a fresh 64-bit word w from the generator, keep its top 53 bits u = ibits(w, 11, 53) (an integer in [0, 253)), and return real(u, dp) * 2**-53. Every integer below 253 is exact in IEEE binary64 and 2-53 is a power of two, so the product is exact: no rounding, no compiler-dependent fused multiply-add, and the result can never reach 1.0. Values are therefore multiples of 2-53 in [0, 1 - 2**-53].

  • private function splitmix64_next_real_dp(gen) result(r)

    Uniform real(dp) in [0, 1) from a SplitMix64 generator.

    Arguments

    Type IntentOptional Attributes Name
    type(splitmix64_t), intent(inout) :: gen

    Return Value real(kind=dp)

  • private function pcg32_next_real_dp(gen) result(r)

    Uniform real(dp) in [0, 1) from a PCG32 generator. Consumes two 32-bit outputs.

    Arguments

    Type IntentOptional Attributes Name
    type(pcg32_t), intent(inout) :: gen

    Return Value real(kind=dp)

public interface stream_for

Derive an independent, reproducible generator for a stream id.

Parallel codes must not let the seed of a task depend on scheduling order. stream_for maps (master_seed, stream_id) to a generator through two SplitMix64 finalisers, so stream 7 is the same sequence whether it runs first, last or on another rank, and different ids give uncorrelated sequences.

  • private subroutine stream_for_splitmix64(master_seed, stream_id, gen)

    Build the SplitMix64 generator belonging to stream_id.

    Arguments

    Type IntentOptional Attributes Name
    integer(kind=int64), intent(in) :: master_seed
    integer(kind=default_int), intent(in) :: stream_id
    type(splitmix64_t), intent(out) :: gen

  • private subroutine stream_for_pcg32(master_seed, stream_id, gen)

    Build the PCG32 generator belonging to stream_id. The derived seed drives a local SplitMix64 which supplies both the LCG state and the stream selector.

    Arguments

    Type IntentOptional Attributes Name
    integer(kind=int64), intent(in) :: master_seed
    integer(kind=default_int), intent(in) :: stream_id
    type(pcg32_t), intent(out) :: gen

Derived Types

type, public ::  pcg32_t

PCG-XSH-RR 64/32 generator: 64-bit LCG state with a 32-bit output permutation, plus a per-stream odd increment.

Components

Type Visibility Attributes Name Initial
integer(kind=int64), public :: increment = PCG_DEFAULT_INCREMENT

Stream increment; always odd, which is what makes streams distinct.
integer(kind=int64), public :: state = 0_int64

Explicit 64-bit LCG state (fixed width, not default_int).

Type-Bound Procedures

procedure, public :: next => pcg32_next
procedure, public :: next_u64 => pcg32_next_u64
procedure, public :: seed => pcg32_seed

type, public ::  splitmix64_t

SplitMix64 generator: a 64-bit counter run through a strong finaliser.

Components

Type Visibility Attributes Name Initial
integer(kind=int64), public :: state = 0_int64

Explicit 64-bit state. Fixed width on purpose: the algorithm is defined on 64-bit words and must not follow default_int.

Type-Bound Procedures

procedure, public :: next => splitmix64_next
procedure, public :: next_u64 => splitmix64_next_u64
procedure, public :: seed => splitmix64_seed

Functions

private pure function mix64(x) result(r)

The SplitMix64 finalising mix (the “MurmurHash3-style” avalanche used by splitmix64.c). Used both as the generator output stage and as the seed derivation function for stream_for.

Arguments

Type IntentOptional Attributes Name
integer(kind=int64), intent(in) :: x

Return Value integer(kind=int64)

private function pcg32_next(self) result(r)

Next 32-bit output, returned as the int32 with the same bit pattern as the reference implementation’s uint32_t result. Roughly half of the returned values are therefore negative; use next_real_dp or next_below if a magnitude is wanted.

Arguments

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

Return Value integer(kind=int32)

private function pcg32_next_below(gen, n, error) result(r)

Uniform integer in [0, n) from a PCG32 generator. Each attempt consumes two 32-bit outputs, so the stream stays word aligned no matter how many rejections occur.

Arguments

Type IntentOptional Attributes Name
type(pcg32_t), intent(inout) :: gen
integer(kind=default_int), intent(in) :: n
type(error_t), intent(inout), optional :: error

Return Value integer(kind=default_int)

private function pcg32_next_real_dp(gen) result(r)

Uniform real(dp) in [0, 1) from a PCG32 generator. Consumes two 32-bit outputs.

Arguments

Type IntentOptional Attributes Name
type(pcg32_t), intent(inout) :: gen

Return Value real(kind=dp)

private function pcg32_next_u64(self) result(r)

Full 64-bit random word built from two consecutive 32-bit outputs; the first output supplies the high half. Fixing the order here is what makes next_real_dp and next_below reproducible for this engine.

Arguments

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

Return Value integer(kind=int64)

private function pcg32_step(self) result(r)

Advance the LCG one step and return the XSH-RR permutation of the previous state as an unsigned value in [0, 2**32).

Arguments

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

Return Value integer(kind=int64)

private pure function real_from_word(word) result(r)

Turn a 64-bit random word into a real(dp) in [0, 1) using its top 53 bits. See the next_real_dp interface documentation for the exactness argument.

Arguments

Type IntentOptional Attributes Name
integer(kind=int64), intent(in) :: word

Return Value real(kind=dp)

private pure function rejection_mask(n) result(mask)

Smallest 2**k - 1 that is at least n - 1, for n >= 1.

Arguments

Type IntentOptional Attributes Name
integer(kind=int64), intent(in) :: n

Return Value integer(kind=int64)

private pure function rotate_right32(x, rot) result(r)

Rotate the low 32 bits of x right by rot bits, 0 <= rot <= 31. Branch free: when rot is zero the left-shifted term is shifted clean out of the low 32 bits and the mask removes it.

Arguments

Type IntentOptional Attributes Name
integer(kind=int64), intent(in) :: x
integer(kind=int64), intent(in) :: rot

Return Value integer(kind=int64)

private function splitmix64_next(self) result(r)

Advance the state by the golden-ratio increment and return the mixed 64-bit output.

Arguments

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

Return Value integer(kind=int64)

private function splitmix64_next_below(gen, n, error) result(r)

Uniform integer in [0, n) from a SplitMix64 generator.

Arguments

Type IntentOptional Attributes Name
type(splitmix64_t), intent(inout) :: gen
integer(kind=default_int), intent(in) :: n
type(error_t), intent(inout), optional :: error

Return Value integer(kind=default_int)

private function splitmix64_next_real_dp(gen) result(r)

Uniform real(dp) in [0, 1) from a SplitMix64 generator.

Arguments

Type IntentOptional Attributes Name
type(splitmix64_t), intent(inout) :: gen

Return Value real(kind=dp)

private function splitmix64_next_u64(self) result(r)

Full 64-bit random word. For SplitMix64 this is simply next; it exists so that both engines present the same interface to the shared real and bounded-integer conversions.

Arguments

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

Return Value integer(kind=int64)

private pure function stream_seed(master_seed, stream_id) result(s)

Deterministic seed for stream_id under master_seed.

Read more…

Arguments

Type IntentOptional Attributes Name
integer(kind=int64), intent(in) :: master_seed
integer(kind=default_int), intent(in) :: stream_id

Return Value integer(kind=int64)

private pure function u32_to_int32(u) result(r)

Reinterpret an unsigned 32-bit value held in [0, 2**32) as the signed int32 with the same bit pattern. Done arithmetically rather than with transfer, and without ever passing an out-of-range value to int.

Arguments

Type IntentOptional Attributes Name
integer(kind=int64), intent(in) :: u

Return Value integer(kind=int32)


Subroutines

private subroutine pcg32_seed(self, state, seq)

Seed the generator, following the reference pcg32_srandom_r: the stream selector is shifted left and forced odd to become the LCG increment, then the state is warmed up, offset by the requested state and warmed up again.

Arguments

Type IntentOptional Attributes Name
class(pcg32_t), intent(inout) :: self
integer(kind=int64), intent(in) :: state
integer(kind=int64), intent(in) :: seq

private subroutine splitmix64_seed(self, s)

Set the generator state. Any 64-bit value is a valid seed, including zero: SplitMix64 has no degenerate state.

Arguments

Type IntentOptional Attributes Name
class(splitmix64_t), intent(inout) :: self
integer(kind=int64), intent(in) :: s

private subroutine stream_for_pcg32(master_seed, stream_id, gen)

Build the PCG32 generator belonging to stream_id. The derived seed drives a local SplitMix64 which supplies both the LCG state and the stream selector.

Arguments

Type IntentOptional Attributes Name
integer(kind=int64), intent(in) :: master_seed
integer(kind=default_int), intent(in) :: stream_id
type(pcg32_t), intent(out) :: gen

private subroutine stream_for_splitmix64(master_seed, stream_id, gen)

Build the SplitMix64 generator belonging to stream_id.

Arguments

Type IntentOptional Attributes Name
integer(kind=int64), intent(in) :: master_seed
integer(kind=default_int), intent(in) :: stream_id
type(splitmix64_t), intent(out) :: gen