[wenzeslaus]

Tools which draw random numbers in parallel, run ensembles, or keep several model components reproducible on their own cannot do that with the one generator shared by the whole program. This adds generators a program owns and layouts which place one generator state per unit of work, such as a raster row or a particle, on the generator's cycle, so that units never share values and a result does not depend on the number of threads. The shared generator and its functions do not change. r.sim.water was the primary motivation, with r.mapcalc as a secondary consideration together with the other existing uses of random numbers in the code base. The design comes from my deep review of the reentrant drand48 streams (#7884) and is meant as their replacement. The documentation page and r.sim.water implementation will come in separate PRs.

Details (AI generated)

This description and the code were written with Claude Code and reviewed by the author of the PR.

What is added

• struct G_random_state, an 8-byte generator state a program owns, and struct G_random_layout, a seed plus the way the generator's span is cut into streams of stride draws, one per unit of work and batch. A generator state is a uint64_t, the counts of a layout are int64_t and its flag is a bool, for which include/grass/gis.h includes <stdint.h> and <stdbool.h>; the public functions take and return int64_t as well, and the messages print it with PRId64, as the vector library does.
• Three initializers: G_random_init_layout(), the basic one, gives the units the whole span, so every unit gets the largest stream the span allows, in a single batch, and it holds at most 2^20 units; G_random_init_layout_exact() takes the number of draws per unit as given (a tool which drew a known number of values per row from the shared generator keeps its earlier values, and gets them for any number of threads); G_random_init_layout_bounded() takes an upper bound and rounds it up to odd.
• Two queries, G_random_layout_batches() (how many batches fit, 0 when none) and G_random_layout_length() (the stride).
• Two placements, G_random_state_for_unit() and G_random_state_for_batch(); a batch is one stream for every unit, and batches are appended, never reserved, so an ensemble member needs only its batch number. Batches start units * stride draws apart, or one more when that number is even: as with the odd strides, the same draw of different batches is then as unrelated as the generator allows, while the relation remains between draws whose numbers differ.
• G_random_advance(), which moves a state by any number of draws in one jump, G_random_double(), and G_random_state_from_seed(), which gives a state the sequence G_srand48() gives the shared generator; such a state is the single unit of a layout of the whole span.
• G_random_generate_seed(), which reads GRASS_RANDOM_SEED, or else SOURCE_DATE_EPOCH, as an integer: an empty value counts as not set, a value which is not an integer or does not fit is a fatal error, a value from -2^31 to 2^32 - 1 is used as it is, and a value outside that range is reduced to its low 32 bits with a warning; without either variable the seed comes from the time and the process ID. G_srand48_auto() now seeds with it and returns that value; before, it returned an environment value unreduced although it had used only the low 32 bits. GRASS_RANDOM_SEED is added to the list of environment variables. In the diff, the time and process ID part of the old G_srand48_auto() body (struct timeval tv;, the gettimeofday() check and the sums) shows as re-indented: it moved out of the else branch of the inline environment check, which is replaced by an early return through the new seed_from_environment() helper with its checks, and its arithmetic widened to uint64_t; the logic of that block is unchanged. The definitions of the multiplier, the increment and the mask also show as moved: they were inside the branch for builds with atomics and are now above it, since the generators a program owns need them in every build, and the longer names of the new constants realign them.

Why the span is a quarter of the cycle

The generator is drand48's 48-bit linear congruential generator. Its multiplier has order 2^46, so two positions 2^46, 2^47 or 3 * 2^46 draws apart give values which differ by a constant (0.25, 0.5, 0.75) at every draw, and the four quarters of the cycle hold the same values shifted. Every layout therefore lives within the first 2^46 draws after the seed, where no two positions are that far apart. Finer relations exist at distances divisible by smaller powers of two (two constants in turn at 2^45, four at 2^44), which is why bounded strides are rounded to odd and a layout of the whole span divides it into an odd number of parts. The stride of such a layout is rounded down to odd and so loses up to two draws per unit, which accumulate along the span; the limit of 2^20 units keeps the units nearest to 2^45 and 2^44 at least 0.48 and 0.24 of a stride away from them. The documentation page gives the table of these relations and how likely a layout meets them.

Seeds

Seeds are int64_t and must be in -2^31 to 2^32 - 1, which is what the shared generator accepts; other values are a fatal error instead of being reduced silently. Seeds are related on this generator (values 2^31 or 2^30 apart give shifted sequences), so the layouts derive everything from one seed rather than from many.

Tests

lib/gis/tests/lib_gis_lrand48_test.py gains 136 tests (143 in the file, about 4 s): the identity of G_random_state_from_seed() with G_srand48(), the three layouts against a reference implementation of the recurrence, pinned values, the batches, every refusal (in subprocesses, since a fatal error ends the process), the absence of constant-shift twins within the span for any unit count, the distance of the units of a layout of the whole span from the two- and four-constant levels for every unit count it accepts, the low-bits rule for odd strides, the advance, the seed from each environment source including the invalid and out-of-range cases, and independence from the shared generator and between threads.

Suggested squash commit message

libgis: Add reproducible random numbers for parallel tools

Tools which draw random numbers in parallel, or run ensembles, or keep
several model components reproducible on their own, cannot do so with
the one generator shared by the whole program: whichever thread asks
first gets the next value, and different seeds are related on this
generator (seeds 2^31 or 2^30 apart are shifted by a constant, and
consecutive seeds form a lattice). Add generators the program owns,
struct G_random_state, and layouts, struct G_random_layout, which place
one state per unit of work on the generator's cycle.

A layout is a seed plus the way the span is cut into streams of stride
draws, one stream per unit and batch. The span is the first quarter of
the cycle after the seed, 2^46 draws: the multiplier has order 2^46, so
positions 2^46, 2^47 or 3 * 2^46 apart give values differing by a
constant, and the other three quarters hold the same values plus 0.25,
0.5 and 0.75. Within the span no two positions are that far apart. The
basic layout, G_random_init_layout(), gives the units the whole span,
every unit the largest stream the span allows, in one batch, dividing
the span into an odd number of parts: with an even number, the unit
halfway along would start 2^45 draws after unit 0 and the unit a quarter
of the way along near 2^44, distances at which the values differ by two
or four constants in turn. It holds at most 2^20 units, because the
stride rounded down to odd loses up to two draws per unit, which
accumulate along the span; up to that count the units nearest to 2^45
and 2^44 start at least 0.48 and 0.24 of a stride away from them. An
exact layout takes the draws of a unit as given, so the streams are the
positions of a serial run and a tool keeps its earlier values; a bounded
layout takes a bound and rounds it up to odd, so that units a power of
two apart relate only a few low bits. A batch is one stream for every
unit, and batches are appended, never reserved, so an ensemble member
needs only its batch number. Batches start units * stride draws apart,
or one more when that number is even, so that the same draw of different
batches is as unrelated as the generator allows, as the odd strides do
for units; the relation remains between draws whose numbers differ, and
a layout gives streams which do not overlap, not independent ones. The
queries return the batches that fit and the stride; the two state calls
place a state on a unit's stream, or a unit's stream of a batch;
G_random_advance() moves a state by any number of draws in one jump;
G_random_state_from_seed() gives the shared generator's sequence, as the
single unit of a layout of the whole span.

Add G_random_generate_seed(), which reads GRASS_RANDOM_SEED, or else
SOURCE_DATE_EPOCH, as an integer: an empty value counts as not set, a
value which is not an integer or does not fit is a fatal error, a value
from -2^31 to 2^32 - 1 is used as it is, and a value outside that range
is reduced to its low 32 bits with a warning. Without either variable
the seed is a weak hash of the time and the process ID. G_srand48_auto()
seeds with that value and returns it; before, it returned an
environment value unreduced although it had used only its low 32 bits.
Add GRASS_RANDOM_SEED to the list of environment variables.

The shared generator is unchanged. Tests cover the identity of
G_random_state_from_seed() with G_srand48(), the three layouts, the
batches, the refusals, the absence of constant-shift twins within the
span, the distance of the units of a layout of the whole span from the
two- and four-constant levels for every unit count it accepts, the
low-bits rule, the advance, the seed from each environment source, and
independence from the shared generator and between threads.

Written with the help of Claude Code.

[nilason]

Any strong reasons not to use int64_t instead of the long long in the new API?

(And uint64_t for the unsigned long long).

[wenzeslaus]

@nilason I don't have a strong opinion either way, but now changed the internal structures to the fixed-width ints. They don't seem to be particularly native in GRASS code base yet and I had to add stdint.h to gis.h. This fits with the RNG on main which uses atomic version of the int, although that one is the "least" version (not the fixed-width). Internally, it would not matter that much which one. The long long int has the same guarantees as the "least" version, and the computation relies only on the minimal length. However, I find fixed-width more readable in this context (you don't have to wonder about the rest of the implementation).

For the API, I kept more general long long. It is what we use everywhere and it makes the messages simpler.

To mirror the stdint.h addition to git.h, I also added stdbool.h to use bool for what was an int used as a boolean.

The use case here is exactly where using them at least internally makes sense, but it requires the standard headers to be included which we avoided so far. The reason I did not consider that before is simply that I'm not used to them and my AI reviewer didn't bring that up either (likely due to the general non-use of them and the need for the header files on top level).

[nilason]

@nilason I don't have a strong opinion either way, but now changed the internal structures to the fixed-width ints. They don't seem to be particularly native in GRASS code base yet and I had to add stdint.h to gis.h. This fits with the RNG on main which uses atomic version of the int, although that one is the "least" version (not the fixed-width). Internally, it would not matter that much which one. The long long int has the same guarantees as the "least" version, and the computation relies only on the minimal length. However, I find fixed-width more readable in this context (you don't have to wonder about the rest of the implementation).

I have currently been made aware of issues of ”unpredictable” size of types, from older standards, in my current work on addressing ”Implicit conversion loses precision” warnings. In the early 90’s when an int could be 16 bit, long was used for assure 32 bit. Well, that is not the case anymore, and what’s worse is that the size varies between data models. The situation with [unsigned] long long is currently more predictable, at least 64 bit practically everywhere. But who knows what’s to come? Internally, if size comes even close to matter, then we should use fixed widths types. But I also think we should not hesitate to use fixed width types in the official API if it makes things clearer (without casting back and forth). We’re C11 now, so including this or that std-header is what’s needed, is the new normal. (Just think of how C++ code looks like regarding fixed-width types).

For the API, I kept more general long long. It is what we use everywhere and it makes the messages simpler.

There are always the PRI• macros (https://cppreference.com/c/types/integer#Format_macro_constants)

To mirror the stdint.h addition to git.h, I also added stdbool.h to use bool for what was an int used as a boolean.

👍

The use case here is exactly where using them at least internally makes sense, but it requires the standard headers to be included which we avoided so far. The reason I did not consider that before is simply that I'm not used to them and my AI reviewer didn't bring that up either (likely due to the general non-use of them and the need for the header files on top level).

No reason to avoid including standard headers (in public GRASS API header files that is, what is really bad is that we include external ones such as GDAL and PROJ), that is unfortunately necessary at this stage. And tell the AI to encourage use and take advantage of C11 (and modern coding practice).

[nilason]

Thanks, looks nice!

(Important for caller code to include stdint.h and not rely on transitive includes via gis.h, as with bool)

[wenzeslaus]

Thank you for the review. Merged. We can still adjust with the documentation or with changes for r.sim.water, r.mapcalc or any other followup tools, but at least r.sim.water should be basically finalized (#7980) and I'm quite happy with the documentation (#8004).