From 9393faaa937bdc8553d702b7481f6cf549c23b3f Mon Sep 17 00:00:00 2001 From: Vaclav Petras Date: Wed, 30 Sep 2026 12:36:10 -0400 Subject: [PATCH 01/12] libgis: Add random number layouts and generator states a program owns 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 run. 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. 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 spread layout gives every unit the largest stream the span allows, for one run, 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. A spread layout 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. Runs are appended blocks of streams, never reserved, so an ensemble needs only the run number. The queries return the runs that fit and the stride; the two state calls place a state on a unit's stream, or a unit's stream of a run; G_random_advance() moves a state by any number of draws in one jump; G_random_seed() gives the shared generator's sequence. 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_seed() with G_srand48(), the three layouts, the runs, the refusals, the absence of constant-shift twins within the span, the distance of spread units from the two- and four-constant levels for every unit count a spread layout 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. --- include/grass/defs/gis.h | 17 + include/grass/gis.h | 20 + lib/gis/lrand48.c | 545 ++++++++++++++- lib/gis/tests/lib_gis_lrand48_test.py | 913 +++++++++++++++++++++++++- lib/init/variables.html | 12 + lib/init/variables.md | 12 + 6 files changed, 1483 insertions(+), 36 deletions(-) diff --git a/include/grass/defs/gis.h b/include/grass/defs/gis.h index 498418fb01e..bd9952d3202 100644 --- a/include/grass/defs/gis.h +++ b/include/grass/defs/gis.h @@ -518,6 +518,23 @@ long G_srand48_auto(void); long G_lrand48(void); long G_mrand48(void); double G_drand48(void); +void G_random_seed(struct G_random_state *, long long); +void G_random_init_layout_exact(struct G_random_layout *, long long, long long, + long long); +void G_random_init_layout_bounded(struct G_random_layout *, long long, + long long, long long); +void G_random_init_layout_spread(struct G_random_layout *, long long, + long long); +long long G_random_layout_runs(const struct G_random_layout *); +long long G_random_layout_length(const struct G_random_layout *); +void G_random_state_for_unit(struct G_random_state *, + const struct G_random_layout *, long long); +void G_random_state_for_run(struct G_random_state *, + const struct G_random_layout *, long long, + long long); +void G_random_advance(struct G_random_state *, long long); +double G_random_double(struct G_random_state *); +long long G_random_generate_seed(void); /* ls.c */ void G_set_ls_filter(int (*)(const char *, void *), void *); diff --git a/include/grass/gis.h b/include/grass/gis.h index 2319e9a8869..9a1e8838454 100644 --- a/include/grass/gis.h +++ b/include/grass/gis.h @@ -628,6 +628,26 @@ struct Popen { int pid; }; +/* A position of a pseudo-random number generator owned by the program; + * see G_random_seed() and G_random_state_for_unit(). The contents are + * private to the library and may change with the generator; read and + * advance the state only through the G_random_*() functions. It is a + * plain value: copying it copies the generator's position. */ +struct G_random_state { + unsigned long long state; +}; + +/* A layout: a seed and the way the generator's span is cut into streams + * of stride draws, one per unit of work and run; see + * G_random_init_layout_exact(). Private to the library. */ +struct G_random_layout { + unsigned long long start; /* generator state at the seed */ + long long units; + long long stride; + long long runs; /* runs that fit into the span, 0 when none */ + int spread; /* a spread layout allows a single run */ +}; + typedef int CELL; typedef double DCELL; typedef float FCELL; diff --git a/lib/gis/lrand48.c b/lib/gis/lrand48.c index 1491dfe63a7..0275f6441d4 100644 --- a/lib/gis/lrand48.c +++ b/lib/gis/lrand48.c @@ -18,21 +18,32 @@ * * The seeding functions are not thread-safe; see G_srand48(). * + * The G_random_*() functions instead advance generators owned by the + * program, one struct G_random_state per unit of work, never shared. They + * need neither atomics nor locks, behave identically on every build, and + * give each unit a sequence fixed by the seed and the unit's number, not + * by the thread that draws it. The library places the units within the + * span, the first 2^46 draws after the seed, a quarter of the generator's + * cycle: a layout gives every unit of every run a stream of stride draws + * in it, and a state is put at the start of its unit's stream. See \ref + * gislib_random_streams for the model. + * * SPDX-FileCopyrightText: 2014-2026 GRASS Development Team * SPDX-License-Identifier: GPL-2.0-or-later * * \authors Glynn Clements, Maris Nartiss (thread safety) */ +#include +#include +#include #include #include #include -#include #if defined(__STDC_VERSION__) && __STDC_VERSION__ >= 201112L && \ !defined(__STDC_NO_ATOMICS__) #include -#include #define LRAND48_ATOMIC 1 #else #define LRAND48_ATOMIC 0 @@ -54,6 +65,20 @@ typedef unsigned short uint16; typedef unsigned int uint32; typedef signed int int32; +#define LCG_A UINT64_C(0x5DEECE66D) +#define LCG_B UINT64_C(0xB) +#define MASK48 UINT64_C(0xFFFFFFFFFFFF) + +/* The multiplier has order 2^46 modulo 2^48, so two states 2^46, 2^47 or + * 3 * 2^46 draws apart give values which differ by a constant for as long + * as they run. The layouts therefore place all their units within the + * first 2^46 draws after the seed, the span. */ +#define LCG_SPAN (UINT64_C(1) << 46) + +/* The most units a spread layout takes; see G_random_init_layout_spread() + * for why. */ +#define SPREAD_MAX_UNITS (1LL << 20) + #if LRAND48_ATOMIC /* The whole 48-bit state is kept in one atomic integer so that it can be @@ -63,10 +88,6 @@ typedef signed int int32; * the result modulo 2^48. */ static atomic_uint_least64_t state; -#define LCG_A UINT64_C(0x5DEECE66D) -#define LCG_B UINT64_C(0xB) -#define MASK48 UINT64_C(0xFFFFFFFFFFFF) - #else static uint16 x0, x1, x2; @@ -106,48 +127,113 @@ void G_srand48(long seedval) seeded = 1; } +/* Read a seed from an environment variable. An unset or empty variable + * gives no seed, so that the caller can try the next source. */ +static int seed_from_environment(const char *name, long long *seed) +{ + const char *text = getenv(name); + char *end; + long long value; + + if (!text || !*text) + return 0; + errno = 0; + value = strtoll(text, &end, 10); + if (errno == ERANGE) + G_fatal_error( + _("Random number seed %s from %s is too large in magnitude " + "to be read"), + text, name); + if (end == text || *end != '\0') + G_fatal_error(_("Random number seed %s from %s is not an integer"), + text, name); + if (value < -(1LL << 31) || value > (1LL << 32) - 1) { + long long reduced = (long long)((unsigned long long)value & 0xFFFFFFFF); + + G_warning(_("Random number seed %s from %s is used as %lld, " + "its low 32 bits"), + text, name, reduced); + value = reduced; + } + *seed = value; + return 1; +} + /*! - * \brief Seed the pseudo-random number generator from the time and PID + * \brief Generate a seed for a random number generator * - * A weak hash of the current time and PID is generated and used to - * seed the PRNG + * The seed is the value of the environment variable GRASS_RANDOM_SEED, or + * of SOURCE_DATE_EPOCH when GRASS_RANDOM_SEED is not set or empty, and + * otherwise a weak hash of the current time and process ID. A value from + * the environment must be a decimal integer within the range of long long, + * with nothing after it (leading white space and a sign are allowed); + * anything else is a fatal error naming the variable. A value from -2^31 + * to 2^32 - 1 is returned as it is, and a value outside that range is + * reduced to its low 32 bits, between 0 and 2^32 - 1, with a warning. The + * result is therefore a seed G_random_seed(), the layout functions and + * G_srand48() accept. Record it, for example in the history of the output + * map, so that the run can be repeated with it as the seed. * - * This function is not thread-safe. In a multi-threaded program, call - * `G_srand48_auto()` once *before* starting the worker threads; it must - * not run concurrently with another thread seeding or generating values. + * The hash of the time and process ID lies between 0 and 2^32 - 1. Two + * calls in one process within the same microsecond, or within the same + * second on systems without gettimeofday(), return the same value, and + * two processes can get the same value too. * - * \return generated seed value passed to G_srand48() + * The function reads the environment and the clock and changes no + * generator; it may issue a warning or end with a fatal error. + * + * \return the seed */ -long G_srand48_auto(void) +long long G_random_generate_seed(void) { - unsigned long seed; - char *grass_random_seed = getenv("GRASS_RANDOM_SEED"); + long long given; + unsigned long long seed; - if (!grass_random_seed) - grass_random_seed = getenv("SOURCE_DATE_EPOCH"); - if (grass_random_seed) { - seed = strtoull(grass_random_seed, NULL, 10); - } - else { - seed = (unsigned long)getpid(); + if (seed_from_environment("GRASS_RANDOM_SEED", &given) || + seed_from_environment("SOURCE_DATE_EPOCH", &given)) + return given; + + seed = (unsigned long long)getpid(); #ifdef HAVE_GETTIMEOFDAY - { - struct timeval tv; + { + struct timeval tv; - if (gettimeofday(&tv, NULL) < 0) - G_fatal_error(_("gettimeofday failed: %s"), strerror(errno)); - seed += (unsigned long)tv.tv_sec; - seed += (unsigned long)tv.tv_usec; - } + if (gettimeofday(&tv, NULL) < 0) + G_fatal_error(_("gettimeofday failed: %s"), strerror(errno)); + seed += (unsigned long long)tv.tv_sec; + seed += (unsigned long long)tv.tv_usec; + } #else - { - time_t t = time(NULL); + { + time_t t = time(NULL); - seed += (unsigned long)t; - } -#endif + seed += (unsigned long long)t; } +#endif + + return (long long)(seed & 0xFFFFFFFF); +} + +/*! + * \brief Seed the pseudo-random number generator from the time and PID + * + * The seed is what G_random_generate_seed() returns: the value of + * GRASS_RANDOM_SEED or SOURCE_DATE_EPOCH, or a weak hash of the current + * time and PID. + * + * This function is not thread-safe. In a multi-threaded program, call + * `G_srand48_auto()` once *before* starting the worker threads; it must + * not run concurrently with another thread seeding or generating values. + * + * \return the seed passed to G_srand48(), between -2^31 and 2^32 - 1; + * where long has 32 bits, a value of 2^31 or more comes back + * negative, which G_srand48() and G_random_seed() read as the same + * seed + */ +long G_srand48_auto(void) +{ + long long seed = G_random_generate_seed(); G_srand48((long)seed); return (long)seed; @@ -272,6 +358,395 @@ double G_drand48(void) #endif } +/* Advance a generator of the program's own by one step. The + * multiplication may wrap around at 2^64; that does not change the result + * modulo 2^48. */ +static unsigned long long lcg_step(unsigned long long x) +{ + return (LCG_A * x + LCG_B) & MASK48; +} + +/* Turn a seed into a generator state the way G_srand48() does: only the + * low 32 bits are used, so a negative seed gives the state of its two's + * complement 32-bit value. */ +static unsigned long long lcg_seed(long long seed) +{ + return (((unsigned long long)seed & 0xFFFFFFFF) << 16) | 0x330E; +} + +/* Advance the generator by an arbitrary number of steps without taking + * them one at a time. One step is the affine map x -> a * x + c, and + * composing two such maps gives another, so the map for `steps` steps is + * built by repeated squaring, as an integer power would be. */ +static unsigned long long lcg_jump(unsigned long long x, + unsigned long long steps) +{ + unsigned long long a_total = 1, c_total = 0; /* the identity map */ + unsigned long long a = LCG_A, c = LCG_B; /* one step */ + + while (steps) { + if (steps & 1) { + c_total = (a * c_total + c) & MASK48; + a_total = (a * a_total) & MASK48; + } + /* Square the map, so a and c then describe twice as many steps. */ + c = (a * c + c) & MASK48; + a = (a * a) & MASK48; + steps >>= 1; + } + + return (a_total * x + c_total) & MASK48; +} + +/* The generator tells seeds apart by their low 32 bits only. Seeds read + * as signed or as unsigned 32-bit values are accepted, so -1 and + * 4294967295 are both accepted and are the same seed; anything else is + * rejected rather than silently losing its high bits. */ +static void check_seed(long long seed) +{ + if (seed < -(1LL << 31) || seed > (1LL << 32) - 1) + G_fatal_error(_("Random number seed %lld is outside the range from " + "-2147483648 to 4294967295 the generator can use"), + seed); +} + +static void check_units(long long units) +{ + if (units <= 0) + G_fatal_error(_("The number of units of a random number layout must " + "be positive, not %lld"), + units); +} + +/* Fill in a layout of units of the given stride, all of one run first, + * then those of the next run. The number of runs that fit is how many + * times units * stride draws fit into the span. The caller has checked + * that units * stride does not overflow. */ +static void init_layout(struct G_random_layout *layout, long long seed, + long long units, long long stride) +{ + layout->start = lcg_seed(seed); + layout->units = units; + layout->stride = stride; + layout->runs = (long long)(LCG_SPAN / ((unsigned long long)units * + (unsigned long long)stride)); + layout->spread = 0; +} + +/*! + * \brief Seed a pseudo-random number generator of the program's own + * + * Puts the state at the seed: it then produces the sequence G_srand48() + * followed by G_drand48() produces for the same seed, so code moving from + * the shared generator to one of its own reproduces its existing results. + * Use it for a single sequence; for one sequence per unit of work, use a + * layout and G_random_state_for_unit(). + * + * A seed outside -2^31 to 2^32 - 1 is a fatal error. A negative seed + * means its two's complement 32-bit value, as for G_srand48(). + * + * Thread-safe as long as no two threads seed the same state. + * + * \param[out] state generator state to seed + * \param[in] seed seed, from -2^31 to 2^32 - 1 + */ +void G_random_seed(struct G_random_state *state, long long seed) +{ + check_seed(seed); + state->state = lcg_seed(seed); +} + +/*! + * \brief Initialize a layout whose units draw an exact number of values + * + * The span is cut into streams of exactly \p draws_per_unit draws, one per + * unit, placed one after another from the seed, and then one after + * another for the following runs. Unit u of run 0 therefore draws what a + * single sequence from G_random_seed() draws at positions u * + * \p draws_per_unit to (u + 1) * \p draws_per_unit - 1, so the units + * together reproduce a serial run whatever order they are processed in. + * A unit which draws more than \p draws_per_unit values runs into the + * next unit's stream; when the number of draws is only bounded, use + * G_random_init_layout_bounded(). + * + * The stride is \p draws_per_unit. G_random_layout_length() returns it, + * and G_random_layout_runs() returns the runs that fit, 2^46 divided by + * \p units * \p draws_per_unit, which is 0 when one run does not fit into + * the span. See \ref gislib_random_streams. + * + * A seed outside -2^31 to 2^32 - 1, a number of units or of draws which + * is not positive, or a product of the two beyond the range of long long + * is a fatal error. + * + * \param[out] layout layout to initialize + * \param[in] seed seed, see G_random_seed() + * \param[in] units number of units of work, positive + * \param[in] draws_per_unit number of values each unit draws, positive + */ +void G_random_init_layout_exact(struct G_random_layout *layout, long long seed, + long long units, long long draws_per_unit) +{ + check_seed(seed); + check_units(units); + if (draws_per_unit <= 0) + G_fatal_error(_("The number of random numbers a unit draws must be " + "positive, not %lld"), + draws_per_unit); + if (units > LLONG_MAX / draws_per_unit) + G_fatal_error(_("A random number layout of %lld units of %lld draws " + "each is too large (the number of draws must not " + "exceed %lld)"), + units, draws_per_unit, LLONG_MAX); + init_layout(layout, seed, units, draws_per_unit); +} + +/*! + * \brief Initialize a layout whose units draw at most a number of values + * + * As G_random_init_layout_exact(), but each unit may draw any number of + * values up to \p max_draws: the streams are \p max_draws draws long, or + * \p max_draws + 1 when \p max_draws is even. With an odd stride, units + * 2^i apart in number start an odd multiple of 2^i draws apart, so the + * differences of their states are fixed in the low i + 2 bits only; an + * even stride would add its own power of two to that count. + * + * The stride is \p max_draws rounded up to odd. G_random_layout_length() + * returns it, and G_random_layout_runs() returns the runs that fit, 2^46 + * divided by \p units times the stride, which is 0 when one run does not + * fit into the span. See \ref gislib_random_streams. + * + * A seed outside -2^31 to 2^32 - 1, a number of units or a bound which is + * not positive, or a product of the units and the stride beyond the range + * of long long is a fatal error. + * + * \param[out] layout layout to initialize + * \param[in] seed seed, see G_random_seed() + * \param[in] units number of units of work, positive + * \param[in] max_draws most values any unit draws, positive + */ +void G_random_init_layout_bounded(struct G_random_layout *layout, + long long seed, long long units, + long long max_draws) +{ + long long stride; + + check_seed(seed); + check_units(units); + if (max_draws <= 0) + G_fatal_error(_("The most random numbers a unit draws must be " + "positive, not %lld"), + max_draws); + /* LLONG_MAX is odd, so an even bound can be rounded up. */ + stride = max_draws % 2 == 0 ? max_draws + 1 : max_draws; + if (units > LLONG_MAX / stride) + G_fatal_error(_("A random number layout of %lld units of at most " + "%lld draws each is too large (the number of units " + "times the bound rounded up to odd must not exceed " + "%lld)"), + units, max_draws, LLONG_MAX); + init_layout(layout, seed, units, stride); +} + +/*! + * \brief Initialize a layout which spreads the units over the whole span + * + * For units whose number of draws is not known in advance. The span is + * divided into parts, as many as there are units, or one more when that + * number is even. The stride is 2^46 / parts rounded down, and then down + * to odd when there is more than one part, and unit u starts u * stride + * draws after the seed. The odd number of parts keeps the unit halfway or + * a quarter of the way along from starting 2^45 or 2^44 draws after unit + * 0, distances at which this generator's values relate (see + * \ref gislib_random_streams): the unit nearest to 2^45 starts at least + * 0.48 of a stride away from it and the unit nearest to 2^44 at least + * 0.24, so a unit meets such a relation only after drawing about a + * quarter of its stride. The odd stride puts units 2^i apart in number an + * odd multiple of 2^i draws apart, as in a bounded layout. A single unit + * keeps the whole span of 2^46 draws. The layout holds a single run. + * + * The number of units is limited to 2^20: the rounded stride loses up to + * two draws per unit, which accumulate along the span, so beyond that + * count the units nearest to 2^45 and 2^44 would drift toward them. A + * bounded layout serves any number of units. + * + * The stride is the length of a part. G_random_layout_length() returns + * it, and G_random_layout_runs() returns 1. + * + * A seed outside -2^31 to 2^32 - 1, or a number of units which is not + * positive or above 2^20, is a fatal error. + * + * \param[out] layout layout to initialize + * \param[in] seed seed, see G_random_seed() + * \param[in] units number of units of work, from 1 to 2^20 + */ +void G_random_init_layout_spread(struct G_random_layout *layout, long long seed, + long long units) +{ + unsigned long long parts, stride; + + check_seed(seed); + check_units(units); + if (units > SPREAD_MAX_UNITS) + G_fatal_error(_("A spread random number layout takes at most %lld " + "units, not %lld; a bounded layout serves any number " + "of units"), + SPREAD_MAX_UNITS, units); + parts = (unsigned long long)units | 1; + stride = LCG_SPAN / parts; + if (parts > 1 && stride % 2 == 0) + stride--; + layout->start = lcg_seed(seed); + layout->units = units; + layout->stride = (long long)stride; + layout->runs = 1; + layout->spread = 1; +} + +/*! + * \brief Return the number of runs that fit into the span + * + * \param[in] layout an initialized layout + * + * \return the runs that fit, 1 for a spread layout, 0 when one run of the + * layout does not fit into the span + */ +long long G_random_layout_runs(const struct G_random_layout *layout) +{ + return layout->runs; +} + +/*! + * \brief Return the number of values a unit may draw + * + * \param[in] layout an initialized layout + * + * \return the stride of the layout, the number of values every unit may + * draw without running into the next unit's stream + */ +long long G_random_layout_length(const struct G_random_layout *layout) +{ + return layout->stride; +} + +/*! + * \brief Put a generator state at the start of a unit's stream in a run + * + * The stream of \p unit in \p run starts (\p run * units + \p unit) * stride + * draws after the seed. The state is set whatever it held before, so the + * same call always restarts the same sequence. See + * \ref gislib_random_streams. + * + * A \p unit outside 0 to units - 1, a negative \p run, a \p run other + * than 0 of a spread layout, or a \p run beyond the runs that fit is a + * fatal error. When no run fits, run 0 is still allowed, and only run 0: + * a tool which warned that its layout does not fit into the span may use + * it as earlier versions did. + * + * Thread-safe as long as no two threads use the same state; the layout + * is only read. + * + * \param[out] state generator state to set + * \param[in] layout an initialized layout + * \param[in] run number of the run, from 0 + * \param[in] unit number of the unit, from 0 to units - 1 + */ +void G_random_state_for_run(struct G_random_state *state, + const struct G_random_layout *layout, long long run, + long long unit) +{ + unsigned long long offset; + + if (unit < 0 || unit >= layout->units) + G_fatal_error(_("Random number unit %lld is out of range " + "(must be between 0 and %lld)"), + unit, layout->units - 1); + if (run < 0) + G_fatal_error(_("Random number run %lld is out of range " + "(must not be negative)"), + run); + if (run > 0 && layout->spread) + G_fatal_error(_("Random number run %lld is out of range " + "(a spread layout holds a single run)"), + run); + if (layout->runs >= 1 && run >= layout->runs) + G_fatal_error(n_("Random number run %lld is out of range " + "(%lld run fits into the generator's span)", + "Random number run %lld is out of range " + "(%lld runs fit into the generator's span)", + (unsigned long)layout->runs), + run, layout->runs); + if (layout->runs == 0 && run > 0) + G_fatal_error(_("Random number run %lld is out of range " + "(no run fits into the generator's span, only run 0 " + "can be used)"), + run); + + /* With run below the runs that fit, the offset is below 2^46; with + * run 0 of a layout which does not fit, it may reach past the span. + * It is below 2^63 either way, since units times stride is. */ + offset = ((unsigned long long)run * (unsigned long long)layout->units + + (unsigned long long)unit) * + (unsigned long long)layout->stride; + state->state = lcg_jump(layout->start, offset); +} + +/*! + * \brief Put a generator state at the start of a unit's stream + * + * The same as G_random_state_for_run() with run 0. + * + * \param[out] state generator state to set + * \param[in] layout an initialized layout + * \param[in] unit number of the unit, from 0 to units - 1 + */ +void G_random_state_for_unit(struct G_random_state *state, + const struct G_random_layout *layout, + long long unit) +{ + G_random_state_for_run(state, layout, 0, unit); +} + +/*! + * \brief Advance a generator of the program's own as if values had been drawn + * + * Moves the state by \p draws draws without taking them one at a time, + * so that the next G_random_double() returns what the draw after those + * would have returned. + * + * A negative \p draws is a fatal error. + * + * \param[in,out] state a seeded generator state + * \param[in] draws number of values to skip, not negative + */ +void G_random_advance(struct G_random_state *state, long long draws) +{ + if (draws < 0) + G_fatal_error(_("Cannot advance a random number generator by %lld " + "draws (the number must not be negative)"), + draws); + state->state = lcg_jump(state->state, (unsigned long long)draws); +} + +/*! + * \brief Generate a floating-point value in the range [0,1) from a + * generator of the program's own + * + * Thread-safe as long as no two threads share a state. Unlike + * G_drand48(), this needs no atomics and so behaves identically on every + * build. + * + * \param[in,out] state generator state, set with G_random_seed(), + * G_random_state_for_unit() or G_random_state_for_run() + * + * \return the generated value + */ +double G_random_double(struct G_random_state *state) +{ + state->state = lcg_step(state->state); + /* The state is below 2^53, so the conversion to double is exact. */ + return (double)state->state / 281474976710656.0; /* 2^48 */ +} + /* Test program diff --git a/lib/gis/tests/lib_gis_lrand48_test.py b/lib/gis/tests/lib_gis_lrand48_test.py index ace570ffe86..8fb3a02ffc0 100644 --- a/lib/gis/tests/lib_gis_lrand48_test.py +++ b/lib/gis/tests/lib_gis_lrand48_test.py @@ -6,11 +6,41 @@ is caught here, which makes this a guard for output compatibility. The reference values below were captured from the current implementation. + +The tests after those cover the generators a program owns, the +G_random_*() functions and their layouts. Their expectations are computed +with a closed-form jump written here, independent of the library's jump, +and the single sequence is cross-checked against the shared generator and +the reference values. """ +import subprocess +import sys +import threading +from ctypes import byref + import pytest -from grass.lib.gis import G_drand48, G_lrand48, G_mrand48, G_srand48 +from grass.lib.gis import ( + G_drand48, + G_lrand48, + G_mrand48, + G_random_advance, + G_random_double, + G_random_generate_seed, + G_random_init_layout_bounded, + G_random_init_layout_exact, + G_random_init_layout_spread, + G_random_layout_length, + G_random_layout_runs, + G_random_seed, + G_random_state_for_run, + G_random_state_for_unit, + G_srand48, + G_srand48_auto, + struct_G_random_layout, + struct_G_random_state, +) # First ten outputs of each generator after G_srand48(seed). The seeds cover # 0, 1, two ordinary values, and the largest 32-bit seed value. Matching the @@ -247,3 +277,884 @@ def test_srand48_is_reproducible(): G_srand48(1337) second = [G_lrand48() for _ in range(20)] assert first == second + + +# The generator constants, as in lrand48.c. Changing the generator changes +# every expectation below, so these tests are meant to fail when that +# happens. +LCG_A = 0x5DEECE66D +LCG_B = 0xB +LCG_MODULUS = 2**48 + +# The part of the cycle the layouts place their streams in: the multiplier +# has order 2^46, so states 2^46 apart give values differing by a constant. +SPAN = 2**46 + +# The most units a spread layout takes. +SPREAD_MAX_UNITS = 2**20 + +# Unit counts for the spread layout: small counts, powers of two, which +# give an even stride before rounding, a hundred thousand, and the two +# largest counts allowed. +SPREAD_UNITS = [2, 3, 4, 6, 64, 4096, 100000, SPREAD_MAX_UNITS - 1, SPREAD_MAX_UNITS] + + +def lcg_jump_reference(state, steps): + """Generator state after the given number of steps, in closed form + + Steps compose to a^n * x + b * (a^n - 1) / (a - 1) modulo 2^48. The + division is exact over the integers but a - 1 has no inverse modulo + 2^48, so a^n is computed modulo 2^48 * (a - 1), which keeps the + quotient correct modulo 2^48. Unlike the library, which composes the + affine map by repeated squaring, this takes the additive term in + closed form and reuses none of the library's code. + """ + power = pow(LCG_A, steps, LCG_MODULUS * (LCG_A - 1)) + return (power * state + LCG_B * ((power - 1) // (LCG_A - 1))) % LCG_MODULUS + + +def seed_state(seed): + """Generator state which seeding with the given value produces""" + return ((seed & 0xFFFFFFFF) << 16) | 0x330E + + +def reference_states(seed, offset, n): + """States of draws offset + 1 to offset + n after the seed, in closed form""" + start = lcg_jump_reference(seed_state(seed), offset) + return [lcg_jump_reference(start, i) for i in range(1, n + 1)] + + +def draw_states(state, n): + """Draw n values from a state and return the generator states behind them + + A value is its state divided by 2^48, which a double holds exactly, so + the multiplication recovers the state without rounding. + """ + return [int(G_random_double(byref(state)) * LCG_MODULUS) for _ in range(n)] + + +def seeded_states(seed, n): + """States of the first n draws of the single sequence of a seed""" + state = struct_G_random_state() + G_random_seed(byref(state), seed) + return draw_states(state, n) + + +def exact_layout(seed, units, draws): + layout = struct_G_random_layout() + G_random_init_layout_exact(byref(layout), seed, units, draws) + return layout + + +def bounded_layout(seed, units, max_draws): + layout = struct_G_random_layout() + G_random_init_layout_bounded(byref(layout), seed, units, max_draws) + return layout + + +def spread_layout(seed, units): + layout = struct_G_random_layout() + G_random_init_layout_spread(byref(layout), seed, units) + return layout + + +def unit_states(layout, unit, n, run=0): + """States of the first n draws of a unit in a run of a layout""" + state = struct_G_random_state() + G_random_state_for_run(byref(state), byref(layout), run, unit) + return draw_states(state, n) + + +def spread_stride(units): + """The spread stride: the span divided by the parts, rounded down to odd. + + The parts are the units, or one more when the units are even, so that + the unit halfway along does not start 2^45 draws after unit 0. + """ + parts = units | 1 + stride = SPAN // parts + return stride - 1 if parts > 1 and stride % 2 == 0 else stride + + +def constant_shift(first, second, max_lag=4): + """Lag at which second is first plus a constant, or None if there is none. + + A constant difference at some lag means the two sequences of states are + one sequence shifted in time and in value, which is what states a + multiple of 2^46 apart produce on this generator. + """ + n = len(first) + for lag in range(-max_lag, max_lag + 1): + pairs = [(first[i], second[i + lag]) for i in range(n) if 0 <= i + lag < n] + if len({(b - a) % LCG_MODULUS for a, b in pairs}) == 1: + return lag + return None + + +# Runs calls given as "function number..." in order, on one state and one +# layout, so that a test can check which call ends in a fatal error. +CALLS_SCRIPT = """ +import sys +from ctypes import byref + +from grass.lib import gis + +state = gis.struct_G_random_state() +layout = gis.struct_G_random_layout() +gis.G_random_seed(byref(state), 1337) +for call in sys.argv[1:]: + name, *numbers = call.split() + function = getattr(gis, name) + numbers = [int(number) for number in numbers] + if name.startswith("G_random_init_layout"): + function(byref(layout), *numbers) + elif name.startswith("G_random_state_for"): + function(byref(state), byref(layout), *numbers) + else: + function(byref(state), *numbers) +print("all calls returned") +""" + + +def run_calls(session, tmp_path, *calls): + """Run the calls in a subprocess and return the completed process. + + A fatal error exits the calling process, so the calls run in a + subprocess, and with a session environment because without GISBASE the + error message is not printed. + """ + script = tmp_path / "calls.py" + script.write_text(CALLS_SCRIPT, encoding="utf-8") + return subprocess.run( + [sys.executable, str(script), *calls], + env=session.env, + capture_output=True, + text=True, + timeout=60, + check=False, + ) + + +def assert_fatal(result, message): + """Check that the calls ended in a fatal error containing the message. + + The library wraps long messages across lines, so whitespace is + normalized before the comparison. + """ + assert result.returncode != 0, result.stdout + assert "all calls returned" not in result.stdout + assert message in " ".join(result.stderr.split()), result.stderr + + +def test_lcg_jump_reference_matches_stepping(): + """The closed form agrees with the generator drawn step by step. + + This ties the reference used below to the shared generator, whose + sequences are pinned in REFERENCE. + """ + G_srand48(1337) + stepped = [int(G_drand48() * LCG_MODULUS) for _ in range(1000)] + assert reference_states(1337, 0, 1000) == stepped + assert [ + state / LCG_MODULUS for state in reference_states(1337, 0, 10) + ] == REFERENCE[1337]["drand48"] + + +@pytest.mark.parametrize("seed", [0, 1, 42, 1337, 2147483647, -1, 4294967295]) +def test_random_seed_matches_shared_generator(seed): + """G_random_seed() gives the sequence G_srand48() and G_drand48() give. + + This is what lets code switch from the shared generator to one of its + own without changing its single-threaded results. + """ + G_srand48(seed) + shared = [G_drand48() for _ in range(100)] + state = struct_G_random_state() + G_random_seed(byref(state), seed) + assert [G_random_double(byref(state)) for _ in range(100)] == shared + if seed in REFERENCE: + assert shared[:10] == REFERENCE[seed]["drand48"] + + +@pytest.mark.parametrize( + ("negative", "equivalent"), [(-1, 4294967295), (-2147483648, 2147483648)] +) +def test_random_negative_seed_is_its_32_bit_value(negative, equivalent): + """A negative seed gives the sequence of its two's complement 32-bit value.""" + assert seeded_states(negative, 20) == seeded_states(equivalent, 20) + + +@pytest.mark.parametrize("seed", [-(2**31), 2**32 - 1]) +def test_random_seed_accepts_the_range_boundaries(seed): + """The seed range is inclusive at both ends.""" + assert seeded_states(seed, 5) == reference_states(seed, 0, 5) + + +@pytest.mark.parametrize("seed", [2**32, -(2**31) - 1]) +@pytest.mark.parametrize( + "call", + [ + "G_random_seed {}", + "G_random_init_layout_exact {} 5 7", + "G_random_init_layout_bounded {} 5 7", + "G_random_init_layout_spread {} 5", + ], + ids=["seed", "exact", "bounded", "spread"], +) +def test_random_seed_outside_range_is_fatal( + xy_session_for_module, tmp_path, seed, call +): + """A seed outside -2^31 to 2^32 - 1 is a fatal error, not used modulo 2^32.""" + result = run_calls(xy_session_for_module, tmp_path, call.format(seed)) + assert_fatal(result, f"seed {seed} is outside the range") + + +@pytest.mark.parametrize("seed", [0, 1337, -1]) +def test_exact_layout_units_follow_the_single_sequence(seed): + """Units of an exact layout draw one sequence in consecutive streams. + + Unit u's first value is draw u * 7 + 1 of the single sequence, and the + last draw of unit u is followed by the first draw of unit u + 1, so + the units together reproduce a serial run. + """ + units, draws = 5, 7 + layout = exact_layout(seed, units, draws) + single = seeded_states(seed, units * draws + 1) + for unit in range(units): + states = unit_states(layout, unit, draws + 1) + assert states == single[unit * draws : (unit + 1) * draws + 1], unit + + +@pytest.mark.parametrize( + ("units", "draws"), + [(5, 7), (1, 1), (1000, 100), (4096, 2**34), (2**23, 2**23), (2**23, 2**23 + 1)], +) +def test_exact_layout_runs_and_length(units, draws): + """The runs that fit are the span divided by units times draws. + + The last pair's product exceeds 2^46, so no run fits and runs() is 0. + The length is the number of draws exactly. + """ + layout = exact_layout(1337, units, draws) + assert G_random_layout_runs(byref(layout)) == SPAN // (units * draws) + assert G_random_layout_length(byref(layout)) == draws + + +def test_exact_layout_runs_follow_each_other(): + """Run r's unit u starts where the single sequence has drawn (r * units + u) * draws.""" + layout = exact_layout(1337, 5, 7) + assert unit_states(layout, 2, 3, run=4) == reference_states( + 1337, (4 * 5 + 2) * 7, 3 + ) + + +def test_exact_layout_which_does_not_fit_allows_run_zero( + xy_session_for_module, tmp_path +): + """With no run fitting, run 0 is still allowed but run 1 is fatal. + + A tool which warned that its layout does not fit into the span may + still use it, as earlier versions did. + """ + units, draws = 2**23, 2**23 + 1 + layout = exact_layout(1337, units, draws) + assert G_random_layout_runs(byref(layout)) == 0 + unit = units - 1 + assert unit_states(layout, unit, 3) == reference_states(1337, unit * draws, 3) + result = run_calls( + xy_session_for_module, + tmp_path, + f"G_random_init_layout_exact 1337 {units} {draws}", + "G_random_state_for_run 1 0", + ) + assert_fatal(result, "run 1 is out of range") + + +def test_exact_layout_run_past_the_runs_that_fit_is_fatal( + xy_session_for_module, tmp_path +): + """A run beyond those that fit is a fatal error naming both numbers.""" + units, draws = 4096, 2**34 + runs = SPAN // (units * draws) + assert runs == 1 + ok = run_calls( + xy_session_for_module, + tmp_path, + f"G_random_init_layout_exact 1337 {units} {draws}", + f"G_random_state_for_run {runs - 1} {units - 1}", + ) + assert ok.returncode == 0, ok.stderr + result = run_calls( + xy_session_for_module, + tmp_path, + f"G_random_init_layout_exact 1337 {units} {draws}", + f"G_random_state_for_run {runs} 0", + ) + assert_fatal(result, f"run {runs} is out of range ({runs} run fits") + + +@pytest.mark.parametrize( + ("call", "message"), + [ + ("G_random_init_layout_exact 1337 0 7", "must be positive, not 0"), + ("G_random_init_layout_exact 1337 -5 7", "must be positive, not -5"), + ("G_random_init_layout_exact 1337 5 0", "must be positive, not 0"), + ("G_random_init_layout_exact 1337 5 -7", "must be positive, not -7"), + ("G_random_init_layout_exact 1337 4294967296 4294967296", "too large"), + ("G_random_init_layout_bounded 1337 5 0", "must be positive, not 0"), + ("G_random_init_layout_bounded 1337 0 7", "must be positive, not 0"), + ("G_random_init_layout_bounded 1337 4294967296 4294967296", "too large"), + ("G_random_init_layout_spread 1337 0", "must be positive, not 0"), + ], +) +def test_layout_with_impossible_counts_is_fatal( + xy_session_for_module, tmp_path, call, message +): + """Counts which are not positive or whose product overflows are fatal.""" + assert_fatal(run_calls(xy_session_for_module, tmp_path, call), message) + + +@pytest.mark.parametrize(("bound", "stride"), [(100, 101), (101, 101), (1, 1), (2, 3)]) +def test_bounded_layout_rounds_an_even_bound_up_to_odd(bound, stride): + """An even bound becomes the next odd number, an odd bound stays.""" + layout = bounded_layout(1337, 1000, bound) + assert G_random_layout_length(byref(layout)) == stride + assert G_random_layout_runs(byref(layout)) == SPAN // (1000 * stride) + + +@pytest.mark.parametrize("seed", [0, 1337, -1]) +def test_bounded_layout_unit_starts(seed): + """Unit u starts u * stride draws after the seed, run 1's (units + u) * stride.""" + units, stride = 1000, 101 + layout = bounded_layout(seed, units, 100) + for unit in (0, 1, 2, 999): + assert unit_states(layout, unit, 3) == reference_states( + seed, unit * stride, 3 + ), unit + assert unit_states(layout, unit, 3, run=1) == reference_states( + seed, (units + unit) * stride, 3 + ), unit + + +@pytest.mark.parametrize("units", SPREAD_UNITS) +def test_spread_layout_stride(units): + """The stride is the span divided by the parts, rounded down to odd.""" + layout = spread_layout(1337, units) + stride = G_random_layout_length(byref(layout)) + assert stride == spread_stride(units) + assert stride % 2 == 1 + assert G_random_layout_runs(byref(layout)) == 1 + + +def test_spread_layout_of_one_unit_keeps_the_span(): + """A single unit keeps the whole span of 2^46 draws. + + Its stream is the single sequence of the seed. + """ + layout = spread_layout(1337, 1) + assert G_random_layout_length(byref(layout)) == SPAN + assert G_random_layout_runs(byref(layout)) == 1 + assert unit_states(layout, 0, 20) == seeded_states(1337, 20) + + +@pytest.mark.parametrize("seed", [0, 1337, -1]) +@pytest.mark.parametrize("units", [4, 4096, 100000, SPREAD_MAX_UNITS]) +def test_spread_layout_unit_starts(seed, units): + """Unit u starts u * stride draws after the seed. + + The first, a middle and the last unit are checked; the last one also + shows that the highest unit is accepted. + """ + layout = spread_layout(seed, units) + stride = spread_stride(units) + for unit in sorted({0, 1, units // 2, units - 1}): + assert unit_states(layout, unit, 2) == reference_states( + seed, unit * stride, 2 + ), unit + + +def test_spread_layout_run_one_is_fatal(xy_session_for_module, tmp_path): + """A spread layout holds a single run.""" + result = run_calls( + xy_session_for_module, + tmp_path, + "G_random_init_layout_spread 1337 4", + "G_random_state_for_run 1 0", + ) + assert_fatal(result, "a spread layout holds a single run") + + +@pytest.mark.parametrize("units", [SPREAD_MAX_UNITS + 1, 2**46 + 1, 2**62]) +def test_spread_layout_with_too_many_units_is_fatal( + xy_session_for_module, tmp_path, units +): + """A spread layout takes at most 2^20 units.""" + result = run_calls( + xy_session_for_module, tmp_path, f"G_random_init_layout_spread 1337 {units}" + ) + assert_fatal(result, f"at most {SPREAD_MAX_UNITS} units, not {units}") + + +def test_spread_layout_of_the_most_units(): + """2^20 units, the most allowed, get a stride of 67,108,799 draws. + + The span divided by 2^20 + 1 parts is just over 67,108,800, which is + even, so the stride is rounded down to the odd 67,108,799. The last + unit starts where the formula says. + """ + layout = spread_layout(1337, SPREAD_MAX_UNITS) + assert G_random_layout_length(byref(layout)) == 67108799 + assert spread_stride(SPREAD_MAX_UNITS) == 67108799 + last = SPREAD_MAX_UNITS - 1 + assert unit_states(layout, last, 3) == reference_states(1337, last * 67108799, 3) + + +def nearest_unit_distance(units, distance): + """Distance of the unit start nearest to a distance, in strides. + + Units start a multiple of the stride apart, so this is also the + nearest any two units come to being that distance apart. + """ + stride = spread_stride(units) + nearest = distance // stride + return min( + abs(k * stride - distance) / stride + for k in (nearest, nearest + 1) + if 1 <= k < units + ) + + +def test_spread_layout_keeps_units_off_the_shifted_distances(): + """No unit starts near 2^45 or 2^44 draws after another unit. + + At those distances the values of two positions differ by two or four + constants in turn. Dividing the span by an even number of parts would + put the unit halfway along 2^45 draws after unit 0 and the unit a + quarter of the way along 2^44 draws after it. With an odd number of + parts, the nearest units are half and a quarter of a stride away from + those distances, less the draws lost by rounding the stride down, + which accumulate along the span. For every allowed number of units, + the nearest unit is at least 0.48 of a stride from 2^45 and 0.24 of a + stride from 2^44. This checks the rule; test_spread_layout_stride ties + the rule to the library. + """ + closest_to_half = min( + nearest_unit_distance(units, 2**45) for units in range(2, SPREAD_MAX_UNITS + 1) + ) + closest_to_quarter = min( + nearest_unit_distance(units, 2**44) for units in range(4, SPREAD_MAX_UNITS + 1) + ) + assert closest_to_half >= 0.48 + assert closest_to_quarter >= 0.24 + + +def layouts_for_span_check(): + """Layouts with at least one run fitting, as (layout, units, stride) + + The stride is computed here from the documented rules, not read from the + layout. + """ + cases = [ + (spread_layout(1337, units), units, spread_stride(units)) + for units in [1, *SPREAD_UNITS] + ] + cases += [ + (exact_layout(1337, units, draws), units, draws) + for units, draws in [(5, 7), (4096, 2**34), (2**23, 2**23)] + ] + cases += [ + (bounded_layout(1337, units, bound), units, bound | 1) + for units, bound in [(1000, 100), (4096, 2**34 - 2), (2**23, 2**23 - 2)] + ] + return cases + + +@pytest.mark.parametrize(("layout", "units", "stride"), layouts_for_span_check()) +def test_layout_lies_within_the_span(layout, units, stride): + """The last stream of the last run that fits ends within the span. + + Every position in a stream is then less than 2^46 draws from every + other, so no two streams can be 2^46, 2^47 or 3 * 2^46 draws apart, + the distances at which this generator repeats its values plus a + constant. + """ + runs = G_random_layout_runs(byref(layout)) + length = G_random_layout_length(byref(layout)) + assert length == stride + assert runs >= 1 + last_start = (runs * units - 1) * stride + assert last_start + length <= SPAN + + +def test_spread_layout_units_have_no_constant_shift(): + """Two units of a spread layout are not one sequence shifted in value. + + The first and the last unit of the largest spread layout start less + than 2^46 draws apart, so no lag relates their values by a constant. + """ + layout = spread_layout(1337, SPREAD_MAX_UNITS) + first = unit_states(layout, 0, 100) + second = unit_states(layout, SPREAD_MAX_UNITS - 1, 100) + assert constant_shift(first, second) is None + + +def test_states_half_a_span_apart_alternate_between_two_shifts(): + """States 2^45 apart differ by one of two constants, alternately. + + This level of relation lies within the span, and the library makes no + promise about it. At draw n, the difference of two states d draws + apart is the multiplier to the power n times the difference at the + start, which for d = 2^45 is 2^45 times an odd number. The multiplier + is 5 modulo 8, so its powers alternate between 1 and 5 modulo 8, and + the difference alternates between two values 2^47 apart. For seed 42 + the values differ by 0.625 and 0.125 in turn. This is documented here, + not claimed away. + """ + plain = reference_states(42, 0, 50) + advanced = reference_states(42, 2**45, 50) + shifts = [ + ((a - p) % LCG_MODULUS) / LCG_MODULUS + for p, a in zip(plain, advanced, strict=True) + ] + assert set(shifts[0::2]) == {0.625} + assert set(shifts[1::2]) == {0.125} + state = struct_G_random_state() + G_random_seed(byref(state), 42) + G_random_advance(byref(state), 2**45) + assert draw_states(state, 50) == advanced + + +def low_bits_relation(first, second): + """Low bits in which the differences of two sequences of states are constant""" + differences = [(b - a) % LCG_MODULUS for a, b in zip(first, second, strict=True)] + bits = 0 + while bits < 48 and len({d % 2 ** (bits + 1) for d in differences}) == 1: + bits += 1 + return bits + + +def test_spread_layout_units_two_strides_apart_relate_three_low_bits(): + """Units 2 * stride apart, stride odd, share only the low 3 bits of a shift. + + A distance of 2^k times an odd number makes the differences of the + two units' states constant in their low k + 2 bits; above those bits + the difference varies from draw to draw. + """ + layout = spread_layout(1337, 4) + assert spread_stride(4) % 2 == 1 + first = unit_states(layout, 0, 200) + second = unit_states(layout, 2, 200) + assert low_bits_relation(first, second) == 3 + + +def test_bounded_layout_even_bound_relates_three_low_bits_too(): + """The stride of an even bound is rounded to odd, so units two apart relate 3 bits. + + Without the rounding, a stride of 100 = 4 * 25 would make units two + apart 8 * 25 draws apart and relate 5 low bits. + """ + layout = bounded_layout(1337, 1000, 100) + first = unit_states(layout, 0, 200) + second = unit_states(layout, 2, 200) + assert low_bits_relation(first, second) == 3 + unrounded = reference_states(1337, 200, 200) + assert low_bits_relation(first, unrounded) == 5 + + +def advanced_states(seed, draws, n): + state = struct_G_random_state() + G_random_seed(byref(state), seed) + G_random_advance(byref(state), draws) + return draw_states(state, n) + + +def test_random_advance_equals_drawing(): + """Advancing gives what drawing and discarding gives; 0 is a no-op.""" + assert advanced_states(1337, 60, 40) == seeded_states(1337, 100)[60:] + assert advanced_states(1337, 0, 10) == seeded_states(1337, 10) + + +def test_random_advance_composes(): + """Two advances are one advance by their sum.""" + state = struct_G_random_state() + G_random_seed(byref(state), 1337) + G_random_advance(byref(state), 25) + G_random_advance(byref(state), 35) + assert draw_states(state, 10) == advanced_states(1337, 60, 10) + + +def test_random_advance_beyond_the_span(): + """Advancing past the span gives the states the closed form gives.""" + draws = SPAN + 12345 + assert advanced_states(1337, draws, 10) == reference_states(1337, draws, 10) + + +def test_random_advance_by_the_cycle_returns_to_the_state(): + """After 2^48 draws, the generator's cycle, the state is the same again.""" + assert advanced_states(1337, LCG_MODULUS, 10) == seeded_states(1337, 10) + + +def test_random_advance_by_the_span_shifts_by_a_quarter(): + """After 2^46 draws the generator repeats its values plus one quarter.""" + plain = seeded_states(1337, 20) + shifted = advanced_states(1337, SPAN, 20) + assert shifted == reference_states(1337, SPAN, 20) + assert {(s - p) % LCG_MODULUS for p, s in zip(plain, shifted, strict=True)} == { + LCG_MODULUS // 4 + } + + +def test_random_advance_negative_is_fatal(xy_session_for_module, tmp_path): + """Advancing by a negative number of draws is a fatal error.""" + result = run_calls(xy_session_for_module, tmp_path, "G_random_advance -1") + assert_fatal(result, "by -1 draws") + + +@pytest.mark.parametrize( + ("layout_call", "unit", "message"), + [ + ("G_random_init_layout_spread 1337 4", 4, "unit 4 is out of range"), + ("G_random_init_layout_spread 1337 4", -1, "unit -1 is out of range"), + ("G_random_init_layout_exact 1337 5 7", 5, "unit 5 is out of range"), + ("G_random_init_layout_bounded 1337 5 7", -1, "unit -1 is out of range"), + ], +) +def test_random_state_for_unit_out_of_range_is_fatal( + xy_session_for_module, tmp_path, layout_call, unit, message +): + """A unit outside 0 to units - 1 is a fatal error.""" + result = run_calls( + xy_session_for_module, + tmp_path, + layout_call, + f"G_random_state_for_unit {unit}", + ) + assert_fatal(result, message) + + +def test_random_state_for_run_negative_run_is_fatal(xy_session_for_module, tmp_path): + """A negative run is a fatal error.""" + result = run_calls( + xy_session_for_module, + tmp_path, + "G_random_init_layout_exact 1337 5 7", + "G_random_state_for_run -1 0", + ) + assert_fatal(result, "run -1 is out of range") + + +def test_random_state_for_unit_is_run_zero(): + """G_random_state_for_unit() gives run 0 of G_random_state_for_run().""" + layout = bounded_layout(1337, 1000, 100) + state = struct_G_random_state() + G_random_state_for_unit(byref(state), byref(layout), 7) + assert draw_states(state, 10) == unit_states(layout, 7, 10) + + +def test_random_generate_seed_is_in_range_without_environment(monkeypatch): + """A seed from the time and process ID is one the generators accept.""" + monkeypatch.delenv("GRASS_RANDOM_SEED", raising=False) + monkeypatch.delenv("SOURCE_DATE_EPOCH", raising=False) + seed = G_random_generate_seed() + assert 0 <= seed < 2**32 + + +def test_random_generate_seed_reads_environment(monkeypatch): + """GRASS_RANDOM_SEED takes precedence over SOURCE_DATE_EPOCH.""" + monkeypatch.setenv("GRASS_RANDOM_SEED", "1337") + monkeypatch.setenv("SOURCE_DATE_EPOCH", "42") + assert G_random_generate_seed() == 1337 + monkeypatch.delenv("GRASS_RANDOM_SEED") + assert G_random_generate_seed() == 42 + + +def test_srand48_auto_uses_the_generated_seed(monkeypatch): + """The shared generator's automatic seed is the generated seed. + + The returned value seeds a generator of the program's own to the + sequence the shared generator then produces. + """ + monkeypatch.setenv("GRASS_RANDOM_SEED", "1337") + seed = G_srand48_auto() + assert seed == 1337 + shared = [int(G_drand48() * LCG_MODULUS) for _ in range(10)] + assert shared == seeded_states(seed, 10) + + +GENERATE_SEED_SCRIPT = """ +from grass.lib.gis import G_random_generate_seed, G_srand48_auto + +print(G_random_generate_seed()) +print(G_srand48_auto()) +""" + + +def generate_seed(session, tmp_path, **variables): + """Run G_random_generate_seed() and G_srand48_auto() in a subprocess. + + The variables replace GRASS_RANDOM_SEED and SOURCE_DATE_EPOCH, which + are unset unless given. A warning or fatal error is written by the + library to the process's standard error and needs a session to be + printed, and a fatal error exits the process. + """ + env = dict(session.env) + env.pop("GRASS_RANDOM_SEED", None) + env.pop("SOURCE_DATE_EPOCH", None) + env.update(variables) + script = tmp_path / "generate_seed.py" + script.write_text(GENERATE_SEED_SCRIPT, encoding="utf-8") + return subprocess.run( + [sys.executable, str(script)], + env=env, + capture_output=True, + text=True, + timeout=60, + check=False, + ) + + +@pytest.mark.parametrize("name", ["GRASS_RANDOM_SEED", "SOURCE_DATE_EPOCH"]) +def test_random_generate_seed_reduces_large_environment_value_with_a_warning( + xy_session_for_module, tmp_path, name +): + """A value past 32 bits is reduced to its low 32 bits with a warning. + + Both the generated seed and the shared generator's automatic seed are + the reduced value. + """ + result = generate_seed(xy_session_for_module, tmp_path, **{name: "5000000000"}) + assert result.returncode == 0, result.stderr + reduced = str(5000000000 % 2**32) + assert result.stdout.split() == [reduced, reduced] + message = " ".join(result.stderr.split()) + assert f"5000000000 from {name}" in message + assert "low 32 bits" in message + + +@pytest.mark.parametrize("value", ["-5", "-2147483648", "4294967295"]) +def test_random_generate_seed_keeps_value_in_range_without_warning( + xy_session_for_module, tmp_path, value +): + """A value from -2^31 to 2^32 - 1 is the seed as it is, also when negative.""" + result = generate_seed(xy_session_for_module, tmp_path, GRASS_RANDOM_SEED=value) + assert result.returncode == 0, result.stderr + assert result.stdout.split()[0] == value + assert "WARNING" not in result.stderr + + +def test_random_generate_seed_skips_empty_variable(xy_session_for_module, tmp_path): + """An empty GRASS_RANDOM_SEED counts as unset, so SOURCE_DATE_EPOCH is used.""" + result = generate_seed( + xy_session_for_module, tmp_path, GRASS_RANDOM_SEED="", SOURCE_DATE_EPOCH="42" + ) + assert result.returncode == 0, result.stderr + assert result.stdout.split() == ["42", "42"] + + +@pytest.mark.parametrize( + ("value", "message"), + [ + ("abc", "is not an integer"), + ("1e9", "is not an integer"), + ("12abc", "is not an integer"), + ("99999999999999999999", "is too large in magnitude"), + ("-99999999999999999999", "is too large in magnitude"), + ], +) +def test_random_generate_seed_rejects_invalid_value( + xy_session_for_module, tmp_path, value, message +): + """A value which is not an integer or overflows is a fatal error naming it.""" + result = generate_seed(xy_session_for_module, tmp_path, GRASS_RANDOM_SEED=value) + assert result.returncode != 0 + assert result.stdout == "" + stderr = " ".join(result.stderr.split()) + assert f"{value} from GRASS_RANDOM_SEED {message}" in stderr + + +def test_random_state_is_independent_of_shared_generator(): + """Drawing from the shared generator between draws changes nothing.""" + layout = spread_layout(1337, 4096) + expected = unit_states(layout, 5, 10) + state = struct_G_random_state() + G_random_state_for_unit(byref(state), byref(layout), 5) + G_srand48(99) + interleaved = [] + for _ in range(10): + G_lrand48() + interleaved.extend(draw_states(state, 1)) + assert interleaved == expected + + +def test_random_states_are_unaffected_by_threading(): + """Each thread drawing from its own state gets the serial values. + + ctypes releases the GIL around each call, so the threads do run the C + code concurrently. + """ + units = 4 + n = 5000 + layout = spread_layout(1337, units) + serial = [unit_states(layout, unit, n) for unit in range(units)] + threaded = [None] * units + + def worker(unit): + threaded[unit] = unit_states(layout, unit, n) + + threads = [threading.Thread(target=worker, args=(unit,)) for unit in range(units)] + for thread in threads: + thread.start() + for thread in threads: + thread.join() + assert threaded == serial + + +# First ten values of unit 3 with seed 1337, computed by the closed form at +# 3 * 17175675903 draws for a spread layout of 4096 units (4097 parts) and at 3 * 101 +# draws for a bounded layout of 1000 units with bound 100. A change of the +# engine or of the layouts changes these. +PINNED = { + "spread": [ + 0.39430735545757045, + 0.6290045178804426, + 0.714269874034521, + 0.5881407609385896, + 0.7378029873502499, + 0.712116402550361, + 0.027046777408891387, + 0.6296825457039468, + 0.13697898780445428, + 0.13722939092471265, + ], + "bounded": [ + 0.6813692756324983, + 0.069333342334577, + 0.21082745238814482, + 0.03296383528511271, + 0.4499312354256588, + 0.5150939051932433, + 0.6799364535715284, + 0.4718208970400859, + 0.8985155631446773, + 0.6221841681079674, + ], +} + + +@pytest.mark.parametrize( + ("kind", "layout", "offset"), + [ + ("spread", spread_layout(1337, 4096), 3 * 17175675903), + ("bounded", bounded_layout(1337, 1000, 100), 3 * 101), + ], + ids=["spread", "bounded"], +) +def test_layout_unit_values_are_pinned(kind, layout, offset): + """Unit 3's first ten values are fixed for a given seed and layout.""" + expected = PINNED[kind] + assert [ + state / LCG_MODULUS for state in reference_states(1337, offset, 10) + ] == expected + state = struct_G_random_state() + G_random_state_for_unit(byref(state), byref(layout), 3) + assert [G_random_double(byref(state)) for _ in range(10)] == expected diff --git a/lib/init/variables.html b/lib/init/variables.html index 15caf78bd31..0099bd29db5 100644 --- a/lib/init/variables.html +++ b/lib/init/variables.html @@ -281,6 +281,18 @@

List of selected (GRASS related) shell environment variables

On Mac OS X this should be the pythonw executable for the wxGUI to work. +
GRASS_RANDOM_SEED
+
[libgis]
+ seed for the random number generator of tools which let the user omit + the seed, such as r.mapcalc and + r.sim.water; the same value + gives the same result in every run. The value must be an integer; + anything else is an error. A value from -2147483648 to 4294967295 is + used as it is, and a value outside that range is reduced to its low 32 + bits with a warning. When the variable is not set or empty, + SOURCE_DATE_EPOCH is read the same way if set, and otherwise the seed + comes from the time and the process ID.
+
GRASS_VECTOR_LOWMEM
[vectorlib]
If the environment variable GRASS_VECTOR_LOWMEM exists, memory diff --git a/lib/init/variables.md b/lib/init/variables.md index 38f8c1969e7..88d5aaa830e 100644 --- a/lib/init/variables.md +++ b/lib/init/variables.md @@ -273,6 +273,18 @@ set to override Python executable. On Mac OS X this should be the `pythonw` executable for the wxGUI to work. +GRASS_RANDOM_SEED +\[libgis\] +seed for the random number generator of tools which let the user omit +the seed, such as *[r.mapcalc](r.mapcalc.md)* and +*[r.sim.water](r.sim.water.md)*; the same value gives the same result +in every run. The value must be an integer; anything else is an error. +A value from -2147483648 to 4294967295 is used as it is, and a value +outside that range is reduced to its low 32 bits with a warning. When +the variable is not set or empty, SOURCE_DATE_EPOCH is read the same +way if set, and otherwise the seed comes from the time and the process +ID. + GRASS_VECTOR_LOWMEM \[vectorlib\] If the environment variable GRASS_VECTOR_LOWMEM exists, memory From 7888ba01d5e73fb48174d71065ca568dc69850a8 Mon Sep 17 00:00:00 2001 From: Vaclav Petras Date: Wed, 30 Sep 2026 14:14:21 -0400 Subject: [PATCH 02/12] libgis: Test the environment seed in a subprocess The tests which set GRASS_RANDOM_SEED and SOURCE_DATE_EPOCH with pytest's monkeypatch in the test process failed on Windows, where the Python interpreter and the GRASS library use separate C runtimes with environments of their own, so a variable set from Python is not what the library's getenv() reads. Run these tests in a subprocess which starts with the variables set, as the other environment tests in the file do, and check the shared generator's sequence after G_srand48_auto() through the values the subprocess prints. Written with the help of Claude Code. --- lib/gis/tests/lib_gis_lrand48_test.py | 78 +++++++++++++++------------ 1 file changed, 43 insertions(+), 35 deletions(-) diff --git a/lib/gis/tests/lib_gis_lrand48_test.py b/lib/gis/tests/lib_gis_lrand48_test.py index 8fb3a02ffc0..eddffff89cb 100644 --- a/lib/gis/tests/lib_gis_lrand48_test.py +++ b/lib/gis/tests/lib_gis_lrand48_test.py @@ -27,7 +27,6 @@ G_mrand48, G_random_advance, G_random_double, - G_random_generate_seed, G_random_init_layout_bounded, G_random_init_layout_exact, G_random_init_layout_spread, @@ -37,7 +36,6 @@ G_random_state_for_run, G_random_state_for_unit, G_srand48, - G_srand48_auto, struct_G_random_layout, struct_G_random_state, ) @@ -951,41 +949,12 @@ def test_random_state_for_unit_is_run_zero(): assert draw_states(state, 10) == unit_states(layout, 7, 10) -def test_random_generate_seed_is_in_range_without_environment(monkeypatch): - """A seed from the time and process ID is one the generators accept.""" - monkeypatch.delenv("GRASS_RANDOM_SEED", raising=False) - monkeypatch.delenv("SOURCE_DATE_EPOCH", raising=False) - seed = G_random_generate_seed() - assert 0 <= seed < 2**32 - - -def test_random_generate_seed_reads_environment(monkeypatch): - """GRASS_RANDOM_SEED takes precedence over SOURCE_DATE_EPOCH.""" - monkeypatch.setenv("GRASS_RANDOM_SEED", "1337") - monkeypatch.setenv("SOURCE_DATE_EPOCH", "42") - assert G_random_generate_seed() == 1337 - monkeypatch.delenv("GRASS_RANDOM_SEED") - assert G_random_generate_seed() == 42 - - -def test_srand48_auto_uses_the_generated_seed(monkeypatch): - """The shared generator's automatic seed is the generated seed. - - The returned value seeds a generator of the program's own to the - sequence the shared generator then produces. - """ - monkeypatch.setenv("GRASS_RANDOM_SEED", "1337") - seed = G_srand48_auto() - assert seed == 1337 - shared = [int(G_drand48() * LCG_MODULUS) for _ in range(10)] - assert shared == seeded_states(seed, 10) - - GENERATE_SEED_SCRIPT = """ -from grass.lib.gis import G_random_generate_seed, G_srand48_auto +from grass.lib.gis import G_drand48, G_random_generate_seed, G_srand48_auto print(G_random_generate_seed()) print(G_srand48_auto()) +print(" ".join(str(int(G_drand48() * 2**48)) for _ in range(10))) """ @@ -1013,6 +982,45 @@ def generate_seed(session, tmp_path, **variables): ) +def test_random_generate_seed_is_in_range_without_environment( + xy_session_for_module, tmp_path +): + """A seed from the time and process ID is one the generators accept.""" + result = generate_seed(xy_session_for_module, tmp_path) + assert result.returncode == 0, result.stderr + generated, automatic = (int(value) for value in result.stdout.split()[:2]) + assert 0 <= generated < 2**32 + assert 0 <= automatic < 2**32 + + +def test_random_generate_seed_reads_environment(xy_session_for_module, tmp_path): + """GRASS_RANDOM_SEED takes precedence over SOURCE_DATE_EPOCH.""" + result = generate_seed( + xy_session_for_module, + tmp_path, + GRASS_RANDOM_SEED="1337", + SOURCE_DATE_EPOCH="42", + ) + assert result.returncode == 0, result.stderr + assert result.stdout.split()[0] == "1337" + result = generate_seed(xy_session_for_module, tmp_path, SOURCE_DATE_EPOCH="42") + assert result.returncode == 0, result.stderr + assert result.stdout.split()[0] == "42" + + +def test_srand48_auto_uses_the_generated_seed(xy_session_for_module, tmp_path): + """The shared generator's automatic seed is the generated seed. + + The returned value seeds a generator of the program's own to the + sequence the shared generator then produces. + """ + result = generate_seed(xy_session_for_module, tmp_path, GRASS_RANDOM_SEED="1337") + assert result.returncode == 0, result.stderr + lines = result.stdout.splitlines() + assert lines[:2] == ["1337", "1337"] + assert [int(state) for state in lines[2].split()] == seeded_states(1337, 10) + + @pytest.mark.parametrize("name", ["GRASS_RANDOM_SEED", "SOURCE_DATE_EPOCH"]) def test_random_generate_seed_reduces_large_environment_value_with_a_warning( xy_session_for_module, tmp_path, name @@ -1025,7 +1033,7 @@ def test_random_generate_seed_reduces_large_environment_value_with_a_warning( result = generate_seed(xy_session_for_module, tmp_path, **{name: "5000000000"}) assert result.returncode == 0, result.stderr reduced = str(5000000000 % 2**32) - assert result.stdout.split() == [reduced, reduced] + assert result.stdout.split()[:2] == [reduced, reduced] message = " ".join(result.stderr.split()) assert f"5000000000 from {name}" in message assert "low 32 bits" in message @@ -1048,7 +1056,7 @@ def test_random_generate_seed_skips_empty_variable(xy_session_for_module, tmp_pa xy_session_for_module, tmp_path, GRASS_RANDOM_SEED="", SOURCE_DATE_EPOCH="42" ) assert result.returncode == 0, result.stderr - assert result.stdout.split() == ["42", "42"] + assert result.stdout.split()[:2] == ["42", "42"] @pytest.mark.parametrize( From c8c87e21d4c58d872fa701af4cfd8c0eff667a29 Mon Sep 17 00:00:00 2001 From: Vaclav Petras Date: Wed, 30 Sep 2026 16:02:10 -0400 Subject: [PATCH 03/12] libgis: Add an author and license lines for the random number layouts Add Vaclav Petras to the authors of lrand48.c and drop the remark in parentheses after Maris Nartiss, since the header now describes the thread safety itself. Add the SPDX lines to the random number test file, whose header changes anyway. Written with the help of Claude Code. --- lib/gis/lrand48.c | 2 +- lib/gis/tests/lib_gis_lrand48_test.py | 3 +++ 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/lib/gis/lrand48.c b/lib/gis/lrand48.c index 0275f6441d4..1e48262dc3f 100644 --- a/lib/gis/lrand48.c +++ b/lib/gis/lrand48.c @@ -31,7 +31,7 @@ * SPDX-FileCopyrightText: 2014-2026 GRASS Development Team * SPDX-License-Identifier: GPL-2.0-or-later * - * \authors Glynn Clements, Maris Nartiss (thread safety) + * \authors Glynn Clements, Maris Nartiss, Vaclav Petras */ #include diff --git a/lib/gis/tests/lib_gis_lrand48_test.py b/lib/gis/tests/lib_gis_lrand48_test.py index eddffff89cb..e5a62b0e29f 100644 --- a/lib/gis/tests/lib_gis_lrand48_test.py +++ b/lib/gis/tests/lib_gis_lrand48_test.py @@ -1,3 +1,6 @@ +# SPDX-FileCopyrightText: 2026 GRASS Development Team +# SPDX-License-Identifier: GPL-2.0-or-later + """Regression tests for the library drand48-family PRNG These tests pin the exact output sequence of G_lrand48, G_mrand48 and From 01a58694f1afa9cd5918f5566eeb36a179428109 Mon Sep 17 00:00:00 2001 From: Vaclav Petras Date: Wed, 30 Sep 2026 16:04:13 -0400 Subject: [PATCH 04/12] libgis: State the condition for thread-safe program-owned generators The file header said the states a program owns are never shared, which the library cannot promise; whether a state is used by one thread at a time is up to the program. Say that instead, as the condition under which the functions need neither atomics nor locks and a unit's values do not depend on the thread which draws them. Written with the help of Claude Code. --- lib/gis/lrand48.c | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/lib/gis/lrand48.c b/lib/gis/lrand48.c index 1e48262dc3f..b1c518385e8 100644 --- a/lib/gis/lrand48.c +++ b/lib/gis/lrand48.c @@ -19,10 +19,11 @@ * The seeding functions are not thread-safe; see G_srand48(). * * The G_random_*() functions instead advance generators owned by the - * program, one struct G_random_state per unit of work, never shared. They - * need neither atomics nor locks, behave identically on every build, and - * give each unit a sequence fixed by the seed and the unit's number, not - * by the thread that draws it. The library places the units within the + * program, one struct G_random_state per unit of work. The program keeps + * a state to one thread at a time; the functions then need neither + * atomics nor locks, behave identically on every build, and give each + * unit a sequence fixed by the seed and the unit's number, not by the + * thread that draws it. The library places the units within the * span, the first 2^46 draws after the seed, a quarter of the generator's * cycle: a layout gives every unit of every run a stream of stride draws * in it, and a state is put at the start of its unit's stream. See \ref From 5f2cd1a7db4220e33f66c14593f421821dd08cdf Mon Sep 17 00:00:00 2001 From: Vaclav Petras Date: Wed, 30 Sep 2026 16:21:00 -0400 Subject: [PATCH 05/12] libgis: Share the seed range test between the seed checks The environment reader and the check of the seeds given to the layout and state functions each spelled out the range from -2^31 to 2^32 - 1. Test it in one place, with the reason for the range, and let both use it. Written with the help of Claude Code. --- lib/gis/lrand48.c | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/lib/gis/lrand48.c b/lib/gis/lrand48.c index b1c518385e8..518b7daecf7 100644 --- a/lib/gis/lrand48.c +++ b/lib/gis/lrand48.c @@ -128,6 +128,14 @@ void G_srand48(long seedval) seeded = 1; } +/* The seeds the generator tells apart: G_srand48() uses the low 32 bits of + * the seed, and this range names every such seed once, as a signed or an + * unsigned 32-bit value. */ +static int seed_in_range(long long seed) +{ + return seed >= -(1LL << 31) && seed <= (1LL << 32) - 1; +} + /* Read a seed from an environment variable. An unset or empty variable * gives no seed, so that the caller can try the next source. */ static int seed_from_environment(const char *name, long long *seed) @@ -148,7 +156,7 @@ static int seed_from_environment(const char *name, long long *seed) if (end == text || *end != '\0') G_fatal_error(_("Random number seed %s from %s is not an integer"), text, name); - if (value < -(1LL << 31) || value > (1LL << 32) - 1) { + if (!seed_in_range(value)) { long long reduced = (long long)((unsigned long long)value & 0xFFFFFFFF); G_warning(_("Random number seed %s from %s is used as %lld, " @@ -405,7 +413,7 @@ static unsigned long long lcg_jump(unsigned long long x, * rejected rather than silently losing its high bits. */ static void check_seed(long long seed) { - if (seed < -(1LL << 31) || seed > (1LL << 32) - 1) + if (!seed_in_range(seed)) G_fatal_error(_("Random number seed %lld is outside the range from " "-2147483648 to 4294967295 the generator can use"), seed); From 852652289d3957e573ec057991a422e3118ff46f Mon Sep 17 00:00:00 2001 From: Vaclav Petras Date: Wed, 30 Sep 2026 16:24:01 -0400 Subject: [PATCH 06/12] libgis: State the overflow check as a precondition of init_layout() The comment said the caller has checked that units * stride does not overflow, which reads as a description rather than as the requirement it is. Say that callers must check it beforehand. Written with the help of Claude Code. --- lib/gis/lrand48.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/gis/lrand48.c b/lib/gis/lrand48.c index 518b7daecf7..48a2270cad5 100644 --- a/lib/gis/lrand48.c +++ b/lib/gis/lrand48.c @@ -429,8 +429,8 @@ static void check_units(long long units) /* Fill in a layout of units of the given stride, all of one run first, * then those of the next run. The number of runs that fit is how many - * times units * stride draws fit into the span. The caller has checked - * that units * stride does not overflow. */ + * times units * stride draws fit into the span. Callers must check + * beforehand that units * stride does not overflow. */ static void init_layout(struct G_random_layout *layout, long long seed, long long units, long long stride) { From cf9f2fe1b45f0ac3dc2cfe3c6f53864bc25011d3 Mon Sep 17 00:00:00 2001 From: Vaclav Petras Date: Thu, 1 Oct 2026 15:44:25 -0400 Subject: [PATCH 07/12] libgis: Rename the random number seeding, whole-span and batch functions G_random_seed() becomes G_random_state_from_seed(): it places a state as G_random_state_for_unit() does, and the old name read as returning a seed next to G_random_generate_seed(). G_random_init_layout_spread() becomes G_random_init_layout(): the layout which gives the units the whole span is the basic one, of which a single seeded state is the case of one unit, while the exact and bounded layouts are those with a known stride. "Spread" also names a process in tools which would call the function, such as r.spread. A run becomes a batch, so G_random_layout_runs() and G_random_state_for_run() become G_random_layout_batches() and G_random_state_for_batch(). A batch of units says what the second index is whether a program uses it for the members of an ensemble, for repeated passes or for the processes of a model, and "run" keeps its ordinary meaning. The internal init_layout() becomes fill_layout(), the fields of the layout and the limit on the units of a layout of the whole span are renamed to match, and the messages say batch. The longer name of the limit realigns the neighboring definitions. Nothing else changes; the tests are renamed with the functions. Written with the help of Claude Code. --- include/grass/defs/gis.h | 13 +- include/grass/gis.h | 16 +- lib/gis/lrand48.c | 192 ++++++++++----------- lib/gis/tests/lib_gis_lrand48_test.py | 234 ++++++++++++++------------ 4 files changed, 236 insertions(+), 219 deletions(-) diff --git a/include/grass/defs/gis.h b/include/grass/defs/gis.h index bd9952d3202..529354c9596 100644 --- a/include/grass/defs/gis.h +++ b/include/grass/defs/gis.h @@ -518,20 +518,19 @@ long G_srand48_auto(void); long G_lrand48(void); long G_mrand48(void); double G_drand48(void); -void G_random_seed(struct G_random_state *, long long); +void G_random_state_from_seed(struct G_random_state *, long long); void G_random_init_layout_exact(struct G_random_layout *, long long, long long, long long); void G_random_init_layout_bounded(struct G_random_layout *, long long, long long, long long); -void G_random_init_layout_spread(struct G_random_layout *, long long, - long long); -long long G_random_layout_runs(const struct G_random_layout *); +void G_random_init_layout(struct G_random_layout *, long long, long long); +long long G_random_layout_batches(const struct G_random_layout *); long long G_random_layout_length(const struct G_random_layout *); void G_random_state_for_unit(struct G_random_state *, const struct G_random_layout *, long long); -void G_random_state_for_run(struct G_random_state *, - const struct G_random_layout *, long long, - long long); +void G_random_state_for_batch(struct G_random_state *, + const struct G_random_layout *, long long, + long long); void G_random_advance(struct G_random_state *, long long); double G_random_double(struct G_random_state *); long long G_random_generate_seed(void); diff --git a/include/grass/gis.h b/include/grass/gis.h index 9a1e8838454..f3cf925dc9c 100644 --- a/include/grass/gis.h +++ b/include/grass/gis.h @@ -629,23 +629,23 @@ struct Popen { }; /* A position of a pseudo-random number generator owned by the program; - * see G_random_seed() and G_random_state_for_unit(). The contents are - * private to the library and may change with the generator; read and - * advance the state only through the G_random_*() functions. It is a - * plain value: copying it copies the generator's position. */ + * see G_random_state_from_seed() and G_random_state_for_unit(). The + * contents are private to the library and may change with the generator; + * read and advance the state only through the G_random_*() functions. It + * is a plain value: copying it copies the generator's position. */ struct G_random_state { unsigned long long state; }; /* A layout: a seed and the way the generator's span is cut into streams - * of stride draws, one per unit of work and run; see - * G_random_init_layout_exact(). Private to the library. */ + * of stride draws, one per unit of work and batch; see + * G_random_init_layout(). Private to the library. */ struct G_random_layout { unsigned long long start; /* generator state at the seed */ long long units; long long stride; - long long runs; /* runs that fit into the span, 0 when none */ - int spread; /* a spread layout allows a single run */ + long long batches; /* batches that fit into the span, 0 when none */ + int whole_span; /* a layout of the whole span holds a single batch */ }; typedef int CELL; diff --git a/lib/gis/lrand48.c b/lib/gis/lrand48.c index 48a2270cad5..1673815ed25 100644 --- a/lib/gis/lrand48.c +++ b/lib/gis/lrand48.c @@ -25,7 +25,7 @@ * unit a sequence fixed by the seed and the unit's number, not by the * thread that draws it. The library places the units within the * span, the first 2^46 draws after the seed, a quarter of the generator's - * cycle: a layout gives every unit of every run a stream of stride draws + * cycle: a layout gives every unit of every batch a stream of stride draws * in it, and a state is put at the start of its unit's stream. See \ref * gislib_random_streams for the model. * @@ -66,19 +66,19 @@ typedef unsigned short uint16; typedef unsigned int uint32; typedef signed int int32; -#define LCG_A UINT64_C(0x5DEECE66D) -#define LCG_B UINT64_C(0xB) -#define MASK48 UINT64_C(0xFFFFFFFFFFFF) +#define LCG_A UINT64_C(0x5DEECE66D) +#define LCG_B UINT64_C(0xB) +#define MASK48 UINT64_C(0xFFFFFFFFFFFF) /* The multiplier has order 2^46 modulo 2^48, so two states 2^46, 2^47 or * 3 * 2^46 draws apart give values which differ by a constant for as long * as they run. The layouts therefore place all their units within the * first 2^46 draws after the seed, the span. */ -#define LCG_SPAN (UINT64_C(1) << 46) +#define LCG_SPAN (UINT64_C(1) << 46) -/* The most units a spread layout takes; see G_random_init_layout_spread() - * for why. */ -#define SPREAD_MAX_UNITS (1LL << 20) +/* The most units a layout of the whole span takes; see + * G_random_init_layout() for why. */ +#define WHOLE_SPAN_MAX_UNITS (1LL << 20) #if LRAND48_ATOMIC @@ -179,9 +179,9 @@ static int seed_from_environment(const char *name, long long *seed) * anything else is a fatal error naming the variable. A value from -2^31 * to 2^32 - 1 is returned as it is, and a value outside that range is * reduced to its low 32 bits, between 0 and 2^32 - 1, with a warning. The - * result is therefore a seed G_random_seed(), the layout functions and - * G_srand48() accept. Record it, for example in the history of the output - * map, so that the run can be repeated with it as the seed. + * result is therefore a seed G_random_state_from_seed(), the layout + * functions and G_srand48() accept. Record it, for example in the history of + * the output map, so that the run can be repeated with it as the seed. * * The hash of the time and process ID lies between 0 and 2^32 - 1. Two * calls in one process within the same microsecond, or within the same @@ -237,8 +237,8 @@ long long G_random_generate_seed(void) * * \return the seed passed to G_srand48(), between -2^31 and 2^32 - 1; * where long has 32 bits, a value of 2^31 or more comes back - * negative, which G_srand48() and G_random_seed() read as the same - * seed + * negative, which G_srand48() and G_random_state_from_seed() read + * as the same seed */ long G_srand48_auto(void) { @@ -427,19 +427,19 @@ static void check_units(long long units) units); } -/* Fill in a layout of units of the given stride, all of one run first, - * then those of the next run. The number of runs that fit is how many +/* Fill in a layout of units of the given stride, all of one batch first, + * then those of the next batch. The number of batches that fit is how many * times units * stride draws fit into the span. Callers must check * beforehand that units * stride does not overflow. */ -static void init_layout(struct G_random_layout *layout, long long seed, +static void fill_layout(struct G_random_layout *layout, long long seed, long long units, long long stride) { layout->start = lcg_seed(seed); layout->units = units; layout->stride = stride; - layout->runs = (long long)(LCG_SPAN / ((unsigned long long)units * - (unsigned long long)stride)); - layout->spread = 0; + layout->batches = (long long)(LCG_SPAN / ((unsigned long long)units * + (unsigned long long)stride)); + layout->whole_span = 0; } /*! @@ -448,8 +448,9 @@ static void init_layout(struct G_random_layout *layout, long long seed, * Puts the state at the seed: it then produces the sequence G_srand48() * followed by G_drand48() produces for the same seed, so code moving from * the shared generator to one of its own reproduces its existing results. - * Use it for a single sequence; for one sequence per unit of work, use a - * layout and G_random_state_for_unit(). + * Use it for a single sequence, which then has the whole span to itself; + * for one sequence per unit of work, use a layout and + * G_random_state_for_unit(). * * A seed outside -2^31 to 2^32 - 1 is a fatal error. A negative seed * means its two's complement 32-bit value, as for G_srand48(). @@ -459,7 +460,7 @@ static void init_layout(struct G_random_layout *layout, long long seed, * \param[out] state generator state to seed * \param[in] seed seed, from -2^31 to 2^32 - 1 */ -void G_random_seed(struct G_random_state *state, long long seed) +void G_random_state_from_seed(struct G_random_state *state, long long seed) { check_seed(seed); state->state = lcg_seed(seed); @@ -470,25 +471,25 @@ void G_random_seed(struct G_random_state *state, long long seed) * * The span is cut into streams of exactly \p draws_per_unit draws, one per * unit, placed one after another from the seed, and then one after - * another for the following runs. Unit u of run 0 therefore draws what a - * single sequence from G_random_seed() draws at positions u * - * \p draws_per_unit to (u + 1) * \p draws_per_unit - 1, so the units + * another for the following batches. Unit u of batch 0 therefore draws + * what a single sequence from G_random_state_from_seed() draws at positions + * u * \p draws_per_unit to (u + 1) * \p draws_per_unit - 1, so the units * together reproduce a serial run whatever order they are processed in. * A unit which draws more than \p draws_per_unit values runs into the * next unit's stream; when the number of draws is only bounded, use * G_random_init_layout_bounded(). * * The stride is \p draws_per_unit. G_random_layout_length() returns it, - * and G_random_layout_runs() returns the runs that fit, 2^46 divided by - * \p units * \p draws_per_unit, which is 0 when one run does not fit into - * the span. See \ref gislib_random_streams. + * and G_random_layout_batches() returns the batches that fit, 2^46 divided + * by \p units * \p draws_per_unit, which is 0 when one batch does not fit + * into the span. See \ref gislib_random_streams. * * A seed outside -2^31 to 2^32 - 1, a number of units or of draws which * is not positive, or a product of the two beyond the range of long long * is a fatal error. * * \param[out] layout layout to initialize - * \param[in] seed seed, see G_random_seed() + * \param[in] seed seed, see G_random_state_from_seed() * \param[in] units number of units of work, positive * \param[in] draws_per_unit number of values each unit draws, positive */ @@ -506,7 +507,7 @@ void G_random_init_layout_exact(struct G_random_layout *layout, long long seed, "each is too large (the number of draws must not " "exceed %lld)"), units, draws_per_unit, LLONG_MAX); - init_layout(layout, seed, units, draws_per_unit); + fill_layout(layout, seed, units, draws_per_unit); } /*! @@ -520,16 +521,16 @@ void G_random_init_layout_exact(struct G_random_layout *layout, long long seed, * even stride would add its own power of two to that count. * * The stride is \p max_draws rounded up to odd. G_random_layout_length() - * returns it, and G_random_layout_runs() returns the runs that fit, 2^46 - * divided by \p units times the stride, which is 0 when one run does not - * fit into the span. See \ref gislib_random_streams. + * returns it, and G_random_layout_batches() returns the batches that fit, + * 2^46 divided by \p units times the stride, which is 0 when one batch does + * not fit into the span. See \ref gislib_random_streams. * * A seed outside -2^31 to 2^32 - 1, a number of units or a bound which is * not positive, or a product of the units and the stride beyond the range * of long long is a fatal error. * * \param[out] layout layout to initialize - * \param[in] seed seed, see G_random_seed() + * \param[in] seed seed, see G_random_state_from_seed() * \param[in] units number of units of work, positive * \param[in] max_draws most values any unit draws, positive */ @@ -553,11 +554,11 @@ void G_random_init_layout_bounded(struct G_random_layout *layout, "times the bound rounded up to odd must not exceed " "%lld)"), units, max_draws, LLONG_MAX); - init_layout(layout, seed, units, stride); + fill_layout(layout, seed, units, stride); } /*! - * \brief Initialize a layout which spreads the units over the whole span + * \brief Initialize a layout which gives the units the whole span * * For units whose number of draws is not known in advance. The span is * divided into parts, as many as there are units, or one more when that @@ -571,7 +572,8 @@ void G_random_init_layout_bounded(struct G_random_layout *layout, * 0.24, so a unit meets such a relation only after drawing about a * quarter of its stride. The odd stride puts units 2^i apart in number an * odd multiple of 2^i draws apart, as in a bounded layout. A single unit - * keeps the whole span of 2^46 draws. The layout holds a single run. + * keeps the whole span of 2^46 draws, as a state from + * G_random_state_from_seed() does. The layout holds a single batch. * * The number of units is limited to 2^20: the rounded stride loses up to * two draws per unit, which accumulate along the span, so beyond that @@ -579,27 +581,27 @@ void G_random_init_layout_bounded(struct G_random_layout *layout, * bounded layout serves any number of units. * * The stride is the length of a part. G_random_layout_length() returns - * it, and G_random_layout_runs() returns 1. + * it, and G_random_layout_batches() returns 1. * * A seed outside -2^31 to 2^32 - 1, or a number of units which is not * positive or above 2^20, is a fatal error. * * \param[out] layout layout to initialize - * \param[in] seed seed, see G_random_seed() + * \param[in] seed seed, see G_random_state_from_seed() * \param[in] units number of units of work, from 1 to 2^20 */ -void G_random_init_layout_spread(struct G_random_layout *layout, long long seed, - long long units) +void G_random_init_layout(struct G_random_layout *layout, long long seed, + long long units) { unsigned long long parts, stride; check_seed(seed); check_units(units); - if (units > SPREAD_MAX_UNITS) - G_fatal_error(_("A spread random number layout takes at most %lld " - "units, not %lld; a bounded layout serves any number " - "of units"), - SPREAD_MAX_UNITS, units); + if (units > WHOLE_SPAN_MAX_UNITS) + G_fatal_error(_("A random number layout of the whole span takes at " + "most %lld units, not %lld; a bounded layout serves " + "any number of units"), + WHOLE_SPAN_MAX_UNITS, units); parts = (unsigned long long)units | 1; stride = LCG_SPAN / parts; if (parts > 1 && stride % 2 == 0) @@ -607,21 +609,21 @@ void G_random_init_layout_spread(struct G_random_layout *layout, long long seed, layout->start = lcg_seed(seed); layout->units = units; layout->stride = (long long)stride; - layout->runs = 1; - layout->spread = 1; + layout->batches = 1; + layout->whole_span = 1; } /*! - * \brief Return the number of runs that fit into the span + * \brief Return the number of batches that fit into the span * * \param[in] layout an initialized layout * - * \return the runs that fit, 1 for a spread layout, 0 when one run of the - * layout does not fit into the span + * \return the batches that fit, 1 for a layout of the whole span, 0 when + * one batch of the layout does not fit into the span */ -long long G_random_layout_runs(const struct G_random_layout *layout) +long long G_random_layout_batches(const struct G_random_layout *layout) { - return layout->runs; + return layout->batches; } /*! @@ -638,30 +640,31 @@ long long G_random_layout_length(const struct G_random_layout *layout) } /*! - * \brief Put a generator state at the start of a unit's stream in a run + * \brief Put a generator state at the start of a unit's stream in a batch * - * The stream of \p unit in \p run starts (\p run * units + \p unit) * stride - * draws after the seed. The state is set whatever it held before, so the - * same call always restarts the same sequence. See - * \ref gislib_random_streams. + * A batch is one stream for every unit of the layout, and the batches + * follow one another along the span: the stream of \p unit in \p batch + * starts (\p batch * units + \p unit) * stride draws after the seed. The + * state is set whatever it held before, so the same call always restarts + * the same sequence. See \ref gislib_random_streams. * - * A \p unit outside 0 to units - 1, a negative \p run, a \p run other - * than 0 of a spread layout, or a \p run beyond the runs that fit is a - * fatal error. When no run fits, run 0 is still allowed, and only run 0: - * a tool which warned that its layout does not fit into the span may use - * it as earlier versions did. + * A \p unit outside 0 to units - 1, a negative \p batch, a \p batch other + * than 0 of a layout of the whole span, or a \p batch beyond the batches + * that fit is a fatal error. When no batch fits, batch 0 is still allowed, + * and only batch 0: a tool which warned that its layout does not fit into + * the span may use it as earlier versions did. * * Thread-safe as long as no two threads use the same state; the layout * is only read. * * \param[out] state generator state to set * \param[in] layout an initialized layout - * \param[in] run number of the run, from 0 + * \param[in] batch number of the batch, from 0 * \param[in] unit number of the unit, from 0 to units - 1 */ -void G_random_state_for_run(struct G_random_state *state, - const struct G_random_layout *layout, long long run, - long long unit) +void G_random_state_for_batch(struct G_random_state *state, + const struct G_random_layout *layout, + long long batch, long long unit) { unsigned long long offset; @@ -669,31 +672,31 @@ void G_random_state_for_run(struct G_random_state *state, G_fatal_error(_("Random number unit %lld is out of range " "(must be between 0 and %lld)"), unit, layout->units - 1); - if (run < 0) - G_fatal_error(_("Random number run %lld is out of range " + if (batch < 0) + G_fatal_error(_("Random number batch %lld is out of range " "(must not be negative)"), - run); - if (run > 0 && layout->spread) - G_fatal_error(_("Random number run %lld is out of range " - "(a spread layout holds a single run)"), - run); - if (layout->runs >= 1 && run >= layout->runs) - G_fatal_error(n_("Random number run %lld is out of range " - "(%lld run fits into the generator's span)", - "Random number run %lld is out of range " - "(%lld runs fit into the generator's span)", - (unsigned long)layout->runs), - run, layout->runs); - if (layout->runs == 0 && run > 0) - G_fatal_error(_("Random number run %lld is out of range " - "(no run fits into the generator's span, only run 0 " - "can be used)"), - run); - - /* With run below the runs that fit, the offset is below 2^46; with - * run 0 of a layout which does not fit, it may reach past the span. + batch); + if (batch > 0 && layout->whole_span) + G_fatal_error(_("Random number batch %lld is out of range " + "(a layout of the whole span holds a single batch)"), + batch); + if (layout->batches >= 1 && batch >= layout->batches) + G_fatal_error(n_("Random number batch %lld is out of range " + "(%lld batch fits into the generator's span)", + "Random number batch %lld is out of range " + "(%lld batches fit into the generator's span)", + (unsigned long)layout->batches), + batch, layout->batches); + if (layout->batches == 0 && batch > 0) + G_fatal_error(_("Random number batch %lld is out of range " + "(no batch fits into the generator's span, only batch " + "0 can be used)"), + batch); + + /* With batch below the batches that fit, the offset is below 2^46; with + * batch 0 of a layout which does not fit, it may reach past the span. * It is below 2^63 either way, since units times stride is. */ - offset = ((unsigned long long)run * (unsigned long long)layout->units + + offset = ((unsigned long long)batch * (unsigned long long)layout->units + (unsigned long long)unit) * (unsigned long long)layout->stride; state->state = lcg_jump(layout->start, offset); @@ -702,7 +705,7 @@ void G_random_state_for_run(struct G_random_state *state, /*! * \brief Put a generator state at the start of a unit's stream * - * The same as G_random_state_for_run() with run 0. + * The same as G_random_state_for_batch() with batch 0. * * \param[out] state generator state to set * \param[in] layout an initialized layout @@ -712,7 +715,7 @@ void G_random_state_for_unit(struct G_random_state *state, const struct G_random_layout *layout, long long unit) { - G_random_state_for_run(state, layout, 0, unit); + G_random_state_for_batch(state, layout, 0, unit); } /*! @@ -744,8 +747,9 @@ void G_random_advance(struct G_random_state *state, long long draws) * G_drand48(), this needs no atomics and so behaves identically on every * build. * - * \param[in,out] state generator state, set with G_random_seed(), - * G_random_state_for_unit() or G_random_state_for_run() + * \param[in,out] state generator state, set with + * G_random_state_from_seed(), G_random_state_for_unit() or + * G_random_state_for_batch() * * \return the generated value */ diff --git a/lib/gis/tests/lib_gis_lrand48_test.py b/lib/gis/tests/lib_gis_lrand48_test.py index e5a62b0e29f..892ecf65e11 100644 --- a/lib/gis/tests/lib_gis_lrand48_test.py +++ b/lib/gis/tests/lib_gis_lrand48_test.py @@ -30,14 +30,14 @@ G_mrand48, G_random_advance, G_random_double, + G_random_init_layout, G_random_init_layout_bounded, G_random_init_layout_exact, - G_random_init_layout_spread, + G_random_layout_batches, G_random_layout_length, - G_random_layout_runs, - G_random_seed, - G_random_state_for_run, + G_random_state_for_batch, G_random_state_for_unit, + G_random_state_from_seed, G_srand48, struct_G_random_layout, struct_G_random_state, @@ -291,13 +291,23 @@ def test_srand48_is_reproducible(): # has order 2^46, so states 2^46 apart give values differing by a constant. SPAN = 2**46 -# The most units a spread layout takes. -SPREAD_MAX_UNITS = 2**20 +# The most units a layout of the whole span takes. +WHOLE_SPAN_MAX_UNITS = 2**20 -# Unit counts for the spread layout: small counts, powers of two, which +# Unit counts for the layout of the whole span: small counts, powers of two, which # give an even stride before rounding, a hundred thousand, and the two # largest counts allowed. -SPREAD_UNITS = [2, 3, 4, 6, 64, 4096, 100000, SPREAD_MAX_UNITS - 1, SPREAD_MAX_UNITS] +WHOLE_SPAN_UNITS = [ + 2, + 3, + 4, + 6, + 64, + 4096, + 100000, + WHOLE_SPAN_MAX_UNITS - 1, + WHOLE_SPAN_MAX_UNITS, +] def lcg_jump_reference(state, steps): @@ -337,7 +347,7 @@ def draw_states(state, n): def seeded_states(seed, n): """States of the first n draws of the single sequence of a seed""" state = struct_G_random_state() - G_random_seed(byref(state), seed) + G_random_state_from_seed(byref(state), seed) return draw_states(state, n) @@ -353,21 +363,21 @@ def bounded_layout(seed, units, max_draws): return layout -def spread_layout(seed, units): +def whole_span_layout(seed, units): layout = struct_G_random_layout() - G_random_init_layout_spread(byref(layout), seed, units) + G_random_init_layout(byref(layout), seed, units) return layout -def unit_states(layout, unit, n, run=0): - """States of the first n draws of a unit in a run of a layout""" +def unit_states(layout, unit, n, batch=0): + """States of the first n draws of a unit in a batch of a layout""" state = struct_G_random_state() - G_random_state_for_run(byref(state), byref(layout), run, unit) + G_random_state_for_batch(byref(state), byref(layout), batch, unit) return draw_states(state, n) -def spread_stride(units): - """The spread stride: the span divided by the parts, rounded down to odd. +def whole_span_stride(units): + """The whole-span stride: the span divided by the parts, rounded down to odd. The parts are the units, or one more when the units are even, so that the unit halfway along does not start 2^45 draws after unit 0. @@ -402,7 +412,7 @@ def constant_shift(first, second, max_lag=4): state = gis.struct_G_random_state() layout = gis.struct_G_random_layout() -gis.G_random_seed(byref(state), 1337) +gis.G_random_state_from_seed(byref(state), 1337) for call in sys.argv[1:]: name, *numbers = call.split() function = getattr(gis, name) @@ -463,7 +473,7 @@ def test_lcg_jump_reference_matches_stepping(): @pytest.mark.parametrize("seed", [0, 1, 42, 1337, 2147483647, -1, 4294967295]) def test_random_seed_matches_shared_generator(seed): - """G_random_seed() gives the sequence G_srand48() and G_drand48() give. + """G_random_state_from_seed() gives the sequence G_srand48() and G_drand48() give. This is what lets code switch from the shared generator to one of its own without changing its single-threaded results. @@ -471,7 +481,7 @@ def test_random_seed_matches_shared_generator(seed): G_srand48(seed) shared = [G_drand48() for _ in range(100)] state = struct_G_random_state() - G_random_seed(byref(state), seed) + G_random_state_from_seed(byref(state), seed) assert [G_random_double(byref(state)) for _ in range(100)] == shared if seed in REFERENCE: assert shared[:10] == REFERENCE[seed]["drand48"] @@ -495,12 +505,12 @@ def test_random_seed_accepts_the_range_boundaries(seed): @pytest.mark.parametrize( "call", [ - "G_random_seed {}", + "G_random_state_from_seed {}", "G_random_init_layout_exact {} 5 7", "G_random_init_layout_bounded {} 5 7", - "G_random_init_layout_spread {} 5", + "G_random_init_layout {} 5", ], - ids=["seed", "exact", "bounded", "spread"], + ids=["seed", "exact", "bounded", "whole_span"], ) def test_random_seed_outside_range_is_fatal( xy_session_for_module, tmp_path, seed, call @@ -530,68 +540,68 @@ def test_exact_layout_units_follow_the_single_sequence(seed): ("units", "draws"), [(5, 7), (1, 1), (1000, 100), (4096, 2**34), (2**23, 2**23), (2**23, 2**23 + 1)], ) -def test_exact_layout_runs_and_length(units, draws): - """The runs that fit are the span divided by units times draws. +def test_exact_layout_batches_and_length(units, draws): + """The batches that fit are the span divided by units times draws. - The last pair's product exceeds 2^46, so no run fits and runs() is 0. + The last pair's product exceeds 2^46, so no batch fits and batches() is 0. The length is the number of draws exactly. """ layout = exact_layout(1337, units, draws) - assert G_random_layout_runs(byref(layout)) == SPAN // (units * draws) + assert G_random_layout_batches(byref(layout)) == SPAN // (units * draws) assert G_random_layout_length(byref(layout)) == draws -def test_exact_layout_runs_follow_each_other(): - """Run r's unit u starts where the single sequence has drawn (r * units + u) * draws.""" +def test_exact_layout_batches_follow_each_other(): + """Unit u of batch b starts (b * units + u) * draws into the single sequence.""" layout = exact_layout(1337, 5, 7) - assert unit_states(layout, 2, 3, run=4) == reference_states( + assert unit_states(layout, 2, 3, batch=4) == reference_states( 1337, (4 * 5 + 2) * 7, 3 ) -def test_exact_layout_which_does_not_fit_allows_run_zero( +def test_exact_layout_which_does_not_fit_allows_batch_zero( xy_session_for_module, tmp_path ): - """With no run fitting, run 0 is still allowed but run 1 is fatal. + """With no batch fitting, batch 0 is still allowed but batch 1 is fatal. A tool which warned that its layout does not fit into the span may still use it, as earlier versions did. """ units, draws = 2**23, 2**23 + 1 layout = exact_layout(1337, units, draws) - assert G_random_layout_runs(byref(layout)) == 0 + assert G_random_layout_batches(byref(layout)) == 0 unit = units - 1 assert unit_states(layout, unit, 3) == reference_states(1337, unit * draws, 3) result = run_calls( xy_session_for_module, tmp_path, f"G_random_init_layout_exact 1337 {units} {draws}", - "G_random_state_for_run 1 0", + "G_random_state_for_batch 1 0", ) - assert_fatal(result, "run 1 is out of range") + assert_fatal(result, "batch 1 is out of range") -def test_exact_layout_run_past_the_runs_that_fit_is_fatal( +def test_exact_layout_batch_past_the_batches_that_fit_is_fatal( xy_session_for_module, tmp_path ): - """A run beyond those that fit is a fatal error naming both numbers.""" + """A batch beyond those that fit is a fatal error naming both numbers.""" units, draws = 4096, 2**34 - runs = SPAN // (units * draws) - assert runs == 1 + batches = SPAN // (units * draws) + assert batches == 1 ok = run_calls( xy_session_for_module, tmp_path, f"G_random_init_layout_exact 1337 {units} {draws}", - f"G_random_state_for_run {runs - 1} {units - 1}", + f"G_random_state_for_batch {batches - 1} {units - 1}", ) assert ok.returncode == 0, ok.stderr result = run_calls( xy_session_for_module, tmp_path, f"G_random_init_layout_exact 1337 {units} {draws}", - f"G_random_state_for_run {runs} 0", + f"G_random_state_for_batch {batches} 0", ) - assert_fatal(result, f"run {runs} is out of range ({runs} run fits") + assert_fatal(result, f"batch {batches} is out of range ({batches} batch fits") @pytest.mark.parametrize( @@ -605,7 +615,7 @@ def test_exact_layout_run_past_the_runs_that_fit_is_fatal( ("G_random_init_layout_bounded 1337 5 0", "must be positive, not 0"), ("G_random_init_layout_bounded 1337 0 7", "must be positive, not 0"), ("G_random_init_layout_bounded 1337 4294967296 4294967296", "too large"), - ("G_random_init_layout_spread 1337 0", "must be positive, not 0"), + ("G_random_init_layout 1337 0", "must be positive, not 0"), ], ) def test_layout_with_impossible_counts_is_fatal( @@ -620,93 +630,93 @@ def test_bounded_layout_rounds_an_even_bound_up_to_odd(bound, stride): """An even bound becomes the next odd number, an odd bound stays.""" layout = bounded_layout(1337, 1000, bound) assert G_random_layout_length(byref(layout)) == stride - assert G_random_layout_runs(byref(layout)) == SPAN // (1000 * stride) + assert G_random_layout_batches(byref(layout)) == SPAN // (1000 * stride) @pytest.mark.parametrize("seed", [0, 1337, -1]) def test_bounded_layout_unit_starts(seed): - """Unit u starts u * stride draws after the seed, run 1's (units + u) * stride.""" + """Unit u starts u * stride draws after the seed, batch 1's (units + u) * stride.""" units, stride = 1000, 101 layout = bounded_layout(seed, units, 100) for unit in (0, 1, 2, 999): assert unit_states(layout, unit, 3) == reference_states( seed, unit * stride, 3 ), unit - assert unit_states(layout, unit, 3, run=1) == reference_states( + assert unit_states(layout, unit, 3, batch=1) == reference_states( seed, (units + unit) * stride, 3 ), unit -@pytest.mark.parametrize("units", SPREAD_UNITS) -def test_spread_layout_stride(units): +@pytest.mark.parametrize("units", WHOLE_SPAN_UNITS) +def test_whole_span_layout_stride(units): """The stride is the span divided by the parts, rounded down to odd.""" - layout = spread_layout(1337, units) + layout = whole_span_layout(1337, units) stride = G_random_layout_length(byref(layout)) - assert stride == spread_stride(units) + assert stride == whole_span_stride(units) assert stride % 2 == 1 - assert G_random_layout_runs(byref(layout)) == 1 + assert G_random_layout_batches(byref(layout)) == 1 -def test_spread_layout_of_one_unit_keeps_the_span(): +def test_whole_span_layout_of_one_unit_keeps_the_span(): """A single unit keeps the whole span of 2^46 draws. Its stream is the single sequence of the seed. """ - layout = spread_layout(1337, 1) + layout = whole_span_layout(1337, 1) assert G_random_layout_length(byref(layout)) == SPAN - assert G_random_layout_runs(byref(layout)) == 1 + assert G_random_layout_batches(byref(layout)) == 1 assert unit_states(layout, 0, 20) == seeded_states(1337, 20) @pytest.mark.parametrize("seed", [0, 1337, -1]) -@pytest.mark.parametrize("units", [4, 4096, 100000, SPREAD_MAX_UNITS]) -def test_spread_layout_unit_starts(seed, units): +@pytest.mark.parametrize("units", [4, 4096, 100000, WHOLE_SPAN_MAX_UNITS]) +def test_whole_span_layout_unit_starts(seed, units): """Unit u starts u * stride draws after the seed. The first, a middle and the last unit are checked; the last one also shows that the highest unit is accepted. """ - layout = spread_layout(seed, units) - stride = spread_stride(units) + layout = whole_span_layout(seed, units) + stride = whole_span_stride(units) for unit in sorted({0, 1, units // 2, units - 1}): assert unit_states(layout, unit, 2) == reference_states( seed, unit * stride, 2 ), unit -def test_spread_layout_run_one_is_fatal(xy_session_for_module, tmp_path): - """A spread layout holds a single run.""" +def test_whole_span_layout_batch_one_is_fatal(xy_session_for_module, tmp_path): + """A layout of the whole span holds a single batch.""" result = run_calls( xy_session_for_module, tmp_path, - "G_random_init_layout_spread 1337 4", - "G_random_state_for_run 1 0", + "G_random_init_layout 1337 4", + "G_random_state_for_batch 1 0", ) - assert_fatal(result, "a spread layout holds a single run") + assert_fatal(result, "a layout of the whole span holds a single batch") -@pytest.mark.parametrize("units", [SPREAD_MAX_UNITS + 1, 2**46 + 1, 2**62]) -def test_spread_layout_with_too_many_units_is_fatal( +@pytest.mark.parametrize("units", [WHOLE_SPAN_MAX_UNITS + 1, 2**46 + 1, 2**62]) +def test_whole_span_layout_with_too_many_units_is_fatal( xy_session_for_module, tmp_path, units ): - """A spread layout takes at most 2^20 units.""" + """A layout of the whole span takes at most 2^20 units.""" result = run_calls( - xy_session_for_module, tmp_path, f"G_random_init_layout_spread 1337 {units}" + xy_session_for_module, tmp_path, f"G_random_init_layout 1337 {units}" ) - assert_fatal(result, f"at most {SPREAD_MAX_UNITS} units, not {units}") + assert_fatal(result, f"at most {WHOLE_SPAN_MAX_UNITS} units, not {units}") -def test_spread_layout_of_the_most_units(): +def test_whole_span_layout_of_the_most_units(): """2^20 units, the most allowed, get a stride of 67,108,799 draws. The span divided by 2^20 + 1 parts is just over 67,108,800, which is even, so the stride is rounded down to the odd 67,108,799. The last unit starts where the formula says. """ - layout = spread_layout(1337, SPREAD_MAX_UNITS) + layout = whole_span_layout(1337, WHOLE_SPAN_MAX_UNITS) assert G_random_layout_length(byref(layout)) == 67108799 - assert spread_stride(SPREAD_MAX_UNITS) == 67108799 - last = SPREAD_MAX_UNITS - 1 + assert whole_span_stride(WHOLE_SPAN_MAX_UNITS) == 67108799 + last = WHOLE_SPAN_MAX_UNITS - 1 assert unit_states(layout, last, 3) == reference_states(1337, last * 67108799, 3) @@ -716,7 +726,7 @@ def nearest_unit_distance(units, distance): Units start a multiple of the stride apart, so this is also the nearest any two units come to being that distance apart. """ - stride = spread_stride(units) + stride = whole_span_stride(units) nearest = distance // stride return min( abs(k * stride - distance) / stride @@ -725,7 +735,7 @@ def nearest_unit_distance(units, distance): ) -def test_spread_layout_keeps_units_off_the_shifted_distances(): +def test_whole_span_layout_keeps_units_off_the_shifted_distances(): """No unit starts near 2^45 or 2^44 draws after another unit. At those distances the values of two positions differ by two or four @@ -736,28 +746,30 @@ def test_spread_layout_keeps_units_off_the_shifted_distances(): those distances, less the draws lost by rounding the stride down, which accumulate along the span. For every allowed number of units, the nearest unit is at least 0.48 of a stride from 2^45 and 0.24 of a - stride from 2^44. This checks the rule; test_spread_layout_stride ties + stride from 2^44. This checks the rule; test_whole_span_layout_stride ties the rule to the library. """ closest_to_half = min( - nearest_unit_distance(units, 2**45) for units in range(2, SPREAD_MAX_UNITS + 1) + nearest_unit_distance(units, 2**45) + for units in range(2, WHOLE_SPAN_MAX_UNITS + 1) ) closest_to_quarter = min( - nearest_unit_distance(units, 2**44) for units in range(4, SPREAD_MAX_UNITS + 1) + nearest_unit_distance(units, 2**44) + for units in range(4, WHOLE_SPAN_MAX_UNITS + 1) ) assert closest_to_half >= 0.48 assert closest_to_quarter >= 0.24 def layouts_for_span_check(): - """Layouts with at least one run fitting, as (layout, units, stride) + """Layouts with at least one batch fitting, as (layout, units, stride) The stride is computed here from the documented rules, not read from the layout. """ cases = [ - (spread_layout(1337, units), units, spread_stride(units)) - for units in [1, *SPREAD_UNITS] + (whole_span_layout(1337, units), units, whole_span_stride(units)) + for units in [1, *WHOLE_SPAN_UNITS] ] cases += [ (exact_layout(1337, units, draws), units, draws) @@ -772,30 +784,30 @@ def layouts_for_span_check(): @pytest.mark.parametrize(("layout", "units", "stride"), layouts_for_span_check()) def test_layout_lies_within_the_span(layout, units, stride): - """The last stream of the last run that fits ends within the span. + """The last stream of the last batch that fits ends within the span. Every position in a stream is then less than 2^46 draws from every other, so no two streams can be 2^46, 2^47 or 3 * 2^46 draws apart, the distances at which this generator repeats its values plus a constant. """ - runs = G_random_layout_runs(byref(layout)) + batches = G_random_layout_batches(byref(layout)) length = G_random_layout_length(byref(layout)) assert length == stride - assert runs >= 1 - last_start = (runs * units - 1) * stride + assert batches >= 1 + last_start = (batches * units - 1) * stride assert last_start + length <= SPAN -def test_spread_layout_units_have_no_constant_shift(): - """Two units of a spread layout are not one sequence shifted in value. +def test_whole_span_layout_units_have_no_constant_shift(): + """Two units of a layout of the whole span are not one sequence shifted in value. - The first and the last unit of the largest spread layout start less + The first and the last unit of the largest layout of the whole span start less than 2^46 draws apart, so no lag relates their values by a constant. """ - layout = spread_layout(1337, SPREAD_MAX_UNITS) + layout = whole_span_layout(1337, WHOLE_SPAN_MAX_UNITS) first = unit_states(layout, 0, 100) - second = unit_states(layout, SPREAD_MAX_UNITS - 1, 100) + second = unit_states(layout, WHOLE_SPAN_MAX_UNITS - 1, 100) assert constant_shift(first, second) is None @@ -820,7 +832,7 @@ def test_states_half_a_span_apart_alternate_between_two_shifts(): assert set(shifts[0::2]) == {0.625} assert set(shifts[1::2]) == {0.125} state = struct_G_random_state() - G_random_seed(byref(state), 42) + G_random_state_from_seed(byref(state), 42) G_random_advance(byref(state), 2**45) assert draw_states(state, 50) == advanced @@ -834,15 +846,15 @@ def low_bits_relation(first, second): return bits -def test_spread_layout_units_two_strides_apart_relate_three_low_bits(): +def test_whole_span_layout_units_two_strides_apart_relate_three_low_bits(): """Units 2 * stride apart, stride odd, share only the low 3 bits of a shift. A distance of 2^k times an odd number makes the differences of the two units' states constant in their low k + 2 bits; above those bits the difference varies from draw to draw. """ - layout = spread_layout(1337, 4) - assert spread_stride(4) % 2 == 1 + layout = whole_span_layout(1337, 4) + assert whole_span_stride(4) % 2 == 1 first = unit_states(layout, 0, 200) second = unit_states(layout, 2, 200) assert low_bits_relation(first, second) == 3 @@ -864,7 +876,7 @@ def test_bounded_layout_even_bound_relates_three_low_bits_too(): def advanced_states(seed, draws, n): state = struct_G_random_state() - G_random_seed(byref(state), seed) + G_random_state_from_seed(byref(state), seed) G_random_advance(byref(state), draws) return draw_states(state, n) @@ -878,7 +890,7 @@ def test_random_advance_equals_drawing(): def test_random_advance_composes(): """Two advances are one advance by their sum.""" state = struct_G_random_state() - G_random_seed(byref(state), 1337) + G_random_state_from_seed(byref(state), 1337) G_random_advance(byref(state), 25) G_random_advance(byref(state), 35) assert draw_states(state, 10) == advanced_states(1337, 60, 10) @@ -914,8 +926,8 @@ def test_random_advance_negative_is_fatal(xy_session_for_module, tmp_path): @pytest.mark.parametrize( ("layout_call", "unit", "message"), [ - ("G_random_init_layout_spread 1337 4", 4, "unit 4 is out of range"), - ("G_random_init_layout_spread 1337 4", -1, "unit -1 is out of range"), + ("G_random_init_layout 1337 4", 4, "unit 4 is out of range"), + ("G_random_init_layout 1337 4", -1, "unit -1 is out of range"), ("G_random_init_layout_exact 1337 5 7", 5, "unit 5 is out of range"), ("G_random_init_layout_bounded 1337 5 7", -1, "unit -1 is out of range"), ], @@ -933,19 +945,21 @@ def test_random_state_for_unit_out_of_range_is_fatal( assert_fatal(result, message) -def test_random_state_for_run_negative_run_is_fatal(xy_session_for_module, tmp_path): - """A negative run is a fatal error.""" +def test_random_state_for_batch_negative_batch_is_fatal( + xy_session_for_module, tmp_path +): + """A negative batch is a fatal error.""" result = run_calls( xy_session_for_module, tmp_path, "G_random_init_layout_exact 1337 5 7", - "G_random_state_for_run -1 0", + "G_random_state_for_batch -1 0", ) - assert_fatal(result, "run -1 is out of range") + assert_fatal(result, "batch -1 is out of range") -def test_random_state_for_unit_is_run_zero(): - """G_random_state_for_unit() gives run 0 of G_random_state_for_run().""" +def test_random_state_for_unit_is_batch_zero(): + """G_random_state_for_unit() gives batch 0 of G_random_state_for_batch().""" layout = bounded_layout(1337, 1000, 100) state = struct_G_random_state() G_random_state_for_unit(byref(state), byref(layout), 7) @@ -1085,7 +1099,7 @@ def test_random_generate_seed_rejects_invalid_value( def test_random_state_is_independent_of_shared_generator(): """Drawing from the shared generator between draws changes nothing.""" - layout = spread_layout(1337, 4096) + layout = whole_span_layout(1337, 4096) expected = unit_states(layout, 5, 10) state = struct_G_random_state() G_random_state_for_unit(byref(state), byref(layout), 5) @@ -1105,7 +1119,7 @@ def test_random_states_are_unaffected_by_threading(): """ units = 4 n = 5000 - layout = spread_layout(1337, units) + layout = whole_span_layout(1337, units) serial = [unit_states(layout, unit, n) for unit in range(units)] threaded = [None] * units @@ -1121,11 +1135,11 @@ def worker(unit): # First ten values of unit 3 with seed 1337, computed by the closed form at -# 3 * 17175675903 draws for a spread layout of 4096 units (4097 parts) and at 3 * 101 -# draws for a bounded layout of 1000 units with bound 100. A change of the -# engine or of the layouts changes these. +# 3 * 17175675903 draws for a layout of the whole span of 4096 units (4097 +# parts) and at 3 * 101 draws for a bounded layout of 1000 units with bound +# 100. A change of the engine or of the layouts changes these. PINNED = { - "spread": [ + "whole_span": [ 0.39430735545757045, 0.6290045178804426, 0.714269874034521, @@ -1155,10 +1169,10 @@ def worker(unit): @pytest.mark.parametrize( ("kind", "layout", "offset"), [ - ("spread", spread_layout(1337, 4096), 3 * 17175675903), + ("whole_span", whole_span_layout(1337, 4096), 3 * 17175675903), ("bounded", bounded_layout(1337, 1000, 100), 3 * 101), ], - ids=["spread", "bounded"], + ids=["whole_span", "bounded"], ) def test_layout_unit_values_are_pinned(kind, layout, offset): """Unit 3's first ten values are fixed for a given seed and layout.""" From dfd0c6b16eb313a07dc73c34f1b521c856a0ba6e Mon Sep 17 00:00:00 2001 From: Vaclav Petras Date: Fri, 2 Oct 2026 13:56:02 -0400 Subject: [PATCH 08/12] libgis: Start the batches of random numbers an odd number of draws apart Batches started units * stride draws apart. When that number is even, the values which different batches give at the same draw are related, the more simply the more often 2 divides the distance: with a distance divisible by 2^20, the value of every batch is given by a quadratic through those of the three batches before it, much as consecutive seeds give values on a straight line. Make the distance units * stride, or one more when that is even, as the strides of bounded layouts and of layouts of the whole span are already odd for the same reason between units. Batch 0 does not change, a batch costs at most one draw more, and no tool uses further batches yet. This does not remove the relation. It moves it from the draws which carry the same number in every batch to draws whose numbers differ by one from batch to batch. A layout gives streams which do not overlap, not independent ones. The batches that fit are now those whose last stream ends within the span, so a batch of exactly 2^46 draws still fits once. Written with the help of Claude Code. --- lib/gis/lrand48.c | 65 ++++++++++++++++++--------- lib/gis/tests/lib_gis_lrand48_test.py | 51 ++++++++++++++------- 2 files changed, 79 insertions(+), 37 deletions(-) diff --git a/lib/gis/lrand48.c b/lib/gis/lrand48.c index 1673815ed25..21f39d0cfd3 100644 --- a/lib/gis/lrand48.c +++ b/lib/gis/lrand48.c @@ -427,18 +427,36 @@ static void check_units(long long units) units); } +/* The draws from the start of one batch to the start of the next: those of + * one batch, units * stride, or one more when that number is even. With an + * odd distance, as with an odd stride between units, the same draw of + * different batches is as unrelated as the generator allows; an even + * distance would relate those draws, the more simply the more often 2 + * divides it. */ +static unsigned long long batch_distance(const struct G_random_layout *layout) +{ + return ((unsigned long long)layout->units * + (unsigned long long)layout->stride) | + 1; +} + /* Fill in a layout of units of the given stride, all of one batch first, - * then those of the next batch. The number of batches that fit is how many - * times units * stride draws fit into the span. Callers must check - * beforehand that units * stride does not overflow. */ + * then those of the next batch. The batches that fit are those whose last + * stream ends within the span. Callers must check beforehand that + * units * stride does not overflow. */ static void fill_layout(struct G_random_layout *layout, long long seed, long long units, long long stride) { + unsigned long long draws = + (unsigned long long)units * (unsigned long long)stride; + layout->start = lcg_seed(seed); layout->units = units; layout->stride = stride; - layout->batches = (long long)(LCG_SPAN / ((unsigned long long)units * - (unsigned long long)stride)); + layout->batches = 0; + if (draws <= LCG_SPAN) + layout->batches = + (long long)((LCG_SPAN - draws) / batch_distance(layout) + 1); layout->whole_span = 0; } @@ -470,19 +488,19 @@ void G_random_state_from_seed(struct G_random_state *state, long long seed) * \brief Initialize a layout whose units draw an exact number of values * * The span is cut into streams of exactly \p draws_per_unit draws, one per - * unit, placed one after another from the seed, and then one after - * another for the following batches. Unit u of batch 0 therefore draws - * what a single sequence from G_random_state_from_seed() draws at positions - * u * \p draws_per_unit to (u + 1) * \p draws_per_unit - 1, so the units - * together reproduce a serial run whatever order they are processed in. + * unit, placed one after another from the seed. Unit u of batch 0 therefore + * draws what a single sequence from G_random_state_from_seed() draws at + * positions u * \p draws_per_unit to (u + 1) * \p draws_per_unit - 1, so + * the units together reproduce a serial run whatever order they are + * processed in. Further batches follow, see G_random_state_for_batch(). * A unit which draws more than \p draws_per_unit values runs into the * next unit's stream; when the number of draws is only bounded, use * G_random_init_layout_bounded(). * * The stride is \p draws_per_unit. G_random_layout_length() returns it, - * and G_random_layout_batches() returns the batches that fit, 2^46 divided - * by \p units * \p draws_per_unit, which is 0 when one batch does not fit - * into the span. See \ref gislib_random_streams. + * and G_random_layout_batches() returns the batches that fit, which is 0 + * when \p units * \p draws_per_unit draws do not fit into the span. See \ref + * gislib_random_streams. * * A seed outside -2^31 to 2^32 - 1, a number of units or of draws which * is not positive, or a product of the two beyond the range of long long @@ -522,8 +540,8 @@ void G_random_init_layout_exact(struct G_random_layout *layout, long long seed, * * The stride is \p max_draws rounded up to odd. G_random_layout_length() * returns it, and G_random_layout_batches() returns the batches that fit, - * 2^46 divided by \p units times the stride, which is 0 when one batch does - * not fit into the span. See \ref gislib_random_streams. + * which is 0 when \p units times the stride draws do not fit into the span. See + * \ref gislib_random_streams. * * A seed outside -2^31 to 2^32 - 1, a number of units or a bound which is * not positive, or a product of the units and the stride beyond the range @@ -643,10 +661,14 @@ long long G_random_layout_length(const struct G_random_layout *layout) * \brief Put a generator state at the start of a unit's stream in a batch * * A batch is one stream for every unit of the layout, and the batches - * follow one another along the span: the stream of \p unit in \p batch - * starts (\p batch * units + \p unit) * stride draws after the seed. The - * state is set whatever it held before, so the same call always restarts - * the same sequence. See \ref gislib_random_streams. + * follow one another along the span. They start units * stride draws + * apart, or one more when that number is even, and the stream of \p unit + * in \p batch starts \p unit * stride draws after the start of the batch. + * With the odd distance, as with an odd stride, the same draw of different + * batches is as unrelated as the generator allows; what an even distance + * would relate there is related between draws whose numbers differ + * instead. The state is set whatever it held before, so the same call + * always restarts the same sequence. See \ref gislib_random_streams. * * A \p unit outside 0 to units - 1, a negative \p batch, a \p batch other * than 0 of a layout of the whole span, or a \p batch beyond the batches @@ -696,9 +718,8 @@ void G_random_state_for_batch(struct G_random_state *state, /* With batch below the batches that fit, the offset is below 2^46; with * batch 0 of a layout which does not fit, it may reach past the span. * It is below 2^63 either way, since units times stride is. */ - offset = ((unsigned long long)batch * (unsigned long long)layout->units + - (unsigned long long)unit) * - (unsigned long long)layout->stride; + offset = (unsigned long long)batch * batch_distance(layout) + + (unsigned long long)unit * (unsigned long long)layout->stride; state->state = lcg_jump(layout->start, offset); } diff --git a/lib/gis/tests/lib_gis_lrand48_test.py b/lib/gis/tests/lib_gis_lrand48_test.py index 892ecf65e11..baecdf6c453 100644 --- a/lib/gis/tests/lib_gis_lrand48_test.py +++ b/lib/gis/tests/lib_gis_lrand48_test.py @@ -369,6 +369,19 @@ def whole_span_layout(seed, units): return layout +def batch_distance(units, stride): + """Draws from the start of one batch to the next: units * stride, made odd""" + return (units * stride) | 1 + + +def batches_that_fit(units, stride): + """Batches whose last stream ends within the span""" + draws = units * stride + if draws > SPAN: + return 0 + return (SPAN - draws) // batch_distance(units, stride) + 1 + + def unit_states(layout, unit, n, batch=0): """States of the first n draws of a unit in a batch of a layout""" state = struct_G_random_state() @@ -541,22 +554,30 @@ def test_exact_layout_units_follow_the_single_sequence(seed): [(5, 7), (1, 1), (1000, 100), (4096, 2**34), (2**23, 2**23), (2**23, 2**23 + 1)], ) def test_exact_layout_batches_and_length(units, draws): - """The batches that fit are the span divided by units times draws. + """The batches that fit are those whose last stream ends within the span. - The last pair's product exceeds 2^46, so no batch fits and batches() is 0. - The length is the number of draws exactly. + A batch of exactly 2^46 draws fits once. The last pair's product exceeds + 2^46, so no batch fits and batches() is 0. The length is the number of + draws exactly. """ layout = exact_layout(1337, units, draws) - assert G_random_layout_batches(byref(layout)) == SPAN // (units * draws) + assert G_random_layout_batches(byref(layout)) == batches_that_fit(units, draws) assert G_random_layout_length(byref(layout)) == draws -def test_exact_layout_batches_follow_each_other(): - """Unit u of batch b starts (b * units + u) * draws into the single sequence.""" - layout = exact_layout(1337, 5, 7) - assert unit_states(layout, 2, 3, batch=4) == reference_states( - 1337, (4 * 5 + 2) * 7, 3 - ) +@pytest.mark.parametrize(("units", "draws", "distance"), [(5, 7, 35), (4, 6, 25)]) +def test_exact_layout_batches_start_an_odd_distance_apart(units, draws, distance): + """Batches start units * draws apart, or one more when that is even. + + Unit u of batch b starts b * distance + u * draws into the single + sequence. + """ + assert batch_distance(units, draws) == distance + layout = exact_layout(1337, units, draws) + for batch, unit in [(0, 2), (1, 0), (4, 2)]: + assert unit_states(layout, unit, 3, batch=batch) == reference_states( + 1337, batch * distance + unit * draws, 3 + ), (batch, unit) def test_exact_layout_which_does_not_fit_allows_batch_zero( @@ -586,7 +607,7 @@ def test_exact_layout_batch_past_the_batches_that_fit_is_fatal( ): """A batch beyond those that fit is a fatal error naming both numbers.""" units, draws = 4096, 2**34 - batches = SPAN // (units * draws) + batches = batches_that_fit(units, draws) assert batches == 1 ok = run_calls( xy_session_for_module, @@ -630,12 +651,12 @@ def test_bounded_layout_rounds_an_even_bound_up_to_odd(bound, stride): """An even bound becomes the next odd number, an odd bound stays.""" layout = bounded_layout(1337, 1000, bound) assert G_random_layout_length(byref(layout)) == stride - assert G_random_layout_batches(byref(layout)) == SPAN // (1000 * stride) + assert G_random_layout_batches(byref(layout)) == batches_that_fit(1000, stride) @pytest.mark.parametrize("seed", [0, 1337, -1]) def test_bounded_layout_unit_starts(seed): - """Unit u starts u * stride draws after the seed, batch 1's (units + u) * stride.""" + """Unit u starts u * stride draws after the seed or after the start of batch 1.""" units, stride = 1000, 101 layout = bounded_layout(seed, units, 100) for unit in (0, 1, 2, 999): @@ -643,7 +664,7 @@ def test_bounded_layout_unit_starts(seed): seed, unit * stride, 3 ), unit assert unit_states(layout, unit, 3, batch=1) == reference_states( - seed, (units + unit) * stride, 3 + seed, batch_distance(units, stride) + unit * stride, 3 ), unit @@ -795,7 +816,7 @@ def test_layout_lies_within_the_span(layout, units, stride): length = G_random_layout_length(byref(layout)) assert length == stride assert batches >= 1 - last_start = (batches * units - 1) * stride + last_start = (batches - 1) * batch_distance(units, stride) + (units - 1) * stride assert last_start + length <= SPAN From 5fd1ddffa80196d1dd2cb68fecd346f94a378e3c Mon Sep 17 00:00:00 2001 From: Vaclav Petras Date: Fri, 2 Oct 2026 15:42:27 -0400 Subject: [PATCH 09/12] libgis: Use uint_least64_t for the state of the program-owned generators The shared generator keeps its state in an atomic_uint_least64_t and advances it in uint_least64_t, with UINT64_C constants. The generators a program owns used unsigned long long for the same 48-bit state and the arithmetic on it. Use uint_least64_t there too, in the helpers and in the two private structs, so that the new code has the type which the file already uses for a generator state. gis.h includes for the two struct fields. It is a C99 header, which the C11 the build requires includes, and sources across the libraries include it already. Seeds, units, draws and batches stay long long: they are range-checked counts which callers parse with strtoll() and messages print with %lld. The typedefs of the 2014 code are left as they are. The values drawn do not change. Written with the help of Claude Code. --- include/grass/gis.h | 5 +++-- lib/gis/lrand48.c | 46 +++++++++++++++++++++------------------------ 2 files changed, 24 insertions(+), 27 deletions(-) diff --git a/include/grass/gis.h b/include/grass/gis.h index f3cf925dc9c..b24c8a0d2d7 100644 --- a/include/grass/gis.h +++ b/include/grass/gis.h @@ -17,6 +17,7 @@ /*============================= Include Files ==============================*/ /* System include files */ +#include #include #include @@ -634,14 +635,14 @@ struct Popen { * read and advance the state only through the G_random_*() functions. It * is a plain value: copying it copies the generator's position. */ struct G_random_state { - unsigned long long state; + uint_least64_t state; }; /* A layout: a seed and the way the generator's span is cut into streams * of stride draws, one per unit of work and batch; see * G_random_init_layout(). Private to the library. */ struct G_random_layout { - unsigned long long start; /* generator state at the seed */ + uint_least64_t start; /* generator state at the seed */ long long units; long long stride; long long batches; /* batches that fit into the span, 0 when none */ diff --git a/lib/gis/lrand48.c b/lib/gis/lrand48.c index 21f39d0cfd3..d39e98ff68b 100644 --- a/lib/gis/lrand48.c +++ b/lib/gis/lrand48.c @@ -157,7 +157,7 @@ static int seed_from_environment(const char *name, long long *seed) G_fatal_error(_("Random number seed %s from %s is not an integer"), text, name); if (!seed_in_range(value)) { - long long reduced = (long long)((unsigned long long)value & 0xFFFFFFFF); + long long reduced = (long long)((uint_least64_t)value & 0xFFFFFFFF); G_warning(_("Random number seed %s from %s is used as %lld, " "its low 32 bits"), @@ -196,13 +196,13 @@ static int seed_from_environment(const char *name, long long *seed) long long G_random_generate_seed(void) { long long given; - unsigned long long seed; + uint_least64_t seed; if (seed_from_environment("GRASS_RANDOM_SEED", &given) || seed_from_environment("SOURCE_DATE_EPOCH", &given)) return given; - seed = (unsigned long long)getpid(); + seed = (uint_least64_t)getpid(); #ifdef HAVE_GETTIMEOFDAY { @@ -210,14 +210,14 @@ long long G_random_generate_seed(void) if (gettimeofday(&tv, NULL) < 0) G_fatal_error(_("gettimeofday failed: %s"), strerror(errno)); - seed += (unsigned long long)tv.tv_sec; - seed += (unsigned long long)tv.tv_usec; + seed += (uint_least64_t)tv.tv_sec; + seed += (uint_least64_t)tv.tv_usec; } #else { time_t t = time(NULL); - seed += (unsigned long long)t; + seed += (uint_least64_t)t; } #endif @@ -370,7 +370,7 @@ double G_drand48(void) /* Advance a generator of the program's own by one step. The * multiplication may wrap around at 2^64; that does not change the result * modulo 2^48. */ -static unsigned long long lcg_step(unsigned long long x) +static uint_least64_t lcg_step(uint_least64_t x) { return (LCG_A * x + LCG_B) & MASK48; } @@ -378,20 +378,19 @@ static unsigned long long lcg_step(unsigned long long x) /* Turn a seed into a generator state the way G_srand48() does: only the * low 32 bits are used, so a negative seed gives the state of its two's * complement 32-bit value. */ -static unsigned long long lcg_seed(long long seed) +static uint_least64_t lcg_seed(long long seed) { - return (((unsigned long long)seed & 0xFFFFFFFF) << 16) | 0x330E; + return (((uint_least64_t)seed & 0xFFFFFFFF) << 16) | 0x330E; } /* Advance the generator by an arbitrary number of steps without taking * them one at a time. One step is the affine map x -> a * x + c, and * composing two such maps gives another, so the map for `steps` steps is * built by repeated squaring, as an integer power would be. */ -static unsigned long long lcg_jump(unsigned long long x, - unsigned long long steps) +static uint_least64_t lcg_jump(uint_least64_t x, uint_least64_t steps) { - unsigned long long a_total = 1, c_total = 0; /* the identity map */ - unsigned long long a = LCG_A, c = LCG_B; /* one step */ + uint_least64_t a_total = 1, c_total = 0; /* the identity map */ + uint_least64_t a = LCG_A, c = LCG_B; /* one step */ while (steps) { if (steps & 1) { @@ -433,11 +432,9 @@ static void check_units(long long units) * different batches is as unrelated as the generator allows; an even * distance would relate those draws, the more simply the more often 2 * divides it. */ -static unsigned long long batch_distance(const struct G_random_layout *layout) +static uint_least64_t batch_distance(const struct G_random_layout *layout) { - return ((unsigned long long)layout->units * - (unsigned long long)layout->stride) | - 1; + return ((uint_least64_t)layout->units * (uint_least64_t)layout->stride) | 1; } /* Fill in a layout of units of the given stride, all of one batch first, @@ -447,8 +444,7 @@ static unsigned long long batch_distance(const struct G_random_layout *layout) static void fill_layout(struct G_random_layout *layout, long long seed, long long units, long long stride) { - unsigned long long draws = - (unsigned long long)units * (unsigned long long)stride; + uint_least64_t draws = (uint_least64_t)units * (uint_least64_t)stride; layout->start = lcg_seed(seed); layout->units = units; @@ -611,7 +607,7 @@ void G_random_init_layout_bounded(struct G_random_layout *layout, void G_random_init_layout(struct G_random_layout *layout, long long seed, long long units) { - unsigned long long parts, stride; + uint_least64_t parts, stride; check_seed(seed); check_units(units); @@ -620,7 +616,7 @@ void G_random_init_layout(struct G_random_layout *layout, long long seed, "most %lld units, not %lld; a bounded layout serves " "any number of units"), WHOLE_SPAN_MAX_UNITS, units); - parts = (unsigned long long)units | 1; + parts = (uint_least64_t)units | 1; stride = LCG_SPAN / parts; if (parts > 1 && stride % 2 == 0) stride--; @@ -688,7 +684,7 @@ void G_random_state_for_batch(struct G_random_state *state, const struct G_random_layout *layout, long long batch, long long unit) { - unsigned long long offset; + uint_least64_t offset; if (unit < 0 || unit >= layout->units) G_fatal_error(_("Random number unit %lld is out of range " @@ -718,8 +714,8 @@ void G_random_state_for_batch(struct G_random_state *state, /* With batch below the batches that fit, the offset is below 2^46; with * batch 0 of a layout which does not fit, it may reach past the span. * It is below 2^63 either way, since units times stride is. */ - offset = (unsigned long long)batch * batch_distance(layout) + - (unsigned long long)unit * (unsigned long long)layout->stride; + offset = (uint_least64_t)batch * batch_distance(layout) + + (uint_least64_t)unit * (uint_least64_t)layout->stride; state->state = lcg_jump(layout->start, offset); } @@ -757,7 +753,7 @@ void G_random_advance(struct G_random_state *state, long long draws) G_fatal_error(_("Cannot advance a random number generator by %lld " "draws (the number must not be negative)"), draws); - state->state = lcg_jump(state->state, (unsigned long long)draws); + state->state = lcg_jump(state->state, (uint_least64_t)draws); } /*! From eebd877f29ba1331c194339d158c8b49722d9029 Mon Sep 17 00:00:00 2001 From: Vaclav Petras Date: Fri, 2 Oct 2026 16:24:51 -0400 Subject: [PATCH 10/12] libgis: Use exact-width integers inside the program-owned generators In the new code, a generator state is a uint64_t and the counts of a layout, its units, stride and batches, are int64_t, in the structs and in the helpers. The previous commit gave the state uint_least64_t, as the shared generator has, which left the layout struct with two spellings of a 64-bit integer side by side. The shared generator keeps its least-width types, those of its atomic state, and its 2014 typedefs. The public functions still take and return long long, which callers parse with strtoll() and messages print with %lld; the conversion at that boundary changes no value. All these names compile to the same code on the platforms the library builds on, so this is a matter of reading, and the values drawn do not change. Written with the help of Claude Code. --- include/grass/gis.h | 12 +++++------ lib/gis/lrand48.c | 50 ++++++++++++++++++++++----------------------- 2 files changed, 31 insertions(+), 31 deletions(-) diff --git a/include/grass/gis.h b/include/grass/gis.h index b24c8a0d2d7..bc1a4c7e6c9 100644 --- a/include/grass/gis.h +++ b/include/grass/gis.h @@ -635,18 +635,18 @@ struct Popen { * read and advance the state only through the G_random_*() functions. It * is a plain value: copying it copies the generator's position. */ struct G_random_state { - uint_least64_t state; + uint64_t state; }; /* A layout: a seed and the way the generator's span is cut into streams * of stride draws, one per unit of work and batch; see * G_random_init_layout(). Private to the library. */ struct G_random_layout { - uint_least64_t start; /* generator state at the seed */ - long long units; - long long stride; - long long batches; /* batches that fit into the span, 0 when none */ - int whole_span; /* a layout of the whole span holds a single batch */ + uint64_t start; /* generator state at the seed */ + int64_t units; + int64_t stride; + int64_t batches; /* batches that fit into the span, 0 when none */ + int whole_span; /* a layout of the whole span holds a single batch */ }; typedef int CELL; diff --git a/lib/gis/lrand48.c b/lib/gis/lrand48.c index d39e98ff68b..28298727ed4 100644 --- a/lib/gis/lrand48.c +++ b/lib/gis/lrand48.c @@ -157,7 +157,7 @@ static int seed_from_environment(const char *name, long long *seed) G_fatal_error(_("Random number seed %s from %s is not an integer"), text, name); if (!seed_in_range(value)) { - long long reduced = (long long)((uint_least64_t)value & 0xFFFFFFFF); + long long reduced = (long long)((uint64_t)value & 0xFFFFFFFF); G_warning(_("Random number seed %s from %s is used as %lld, " "its low 32 bits"), @@ -196,13 +196,13 @@ static int seed_from_environment(const char *name, long long *seed) long long G_random_generate_seed(void) { long long given; - uint_least64_t seed; + uint64_t seed; if (seed_from_environment("GRASS_RANDOM_SEED", &given) || seed_from_environment("SOURCE_DATE_EPOCH", &given)) return given; - seed = (uint_least64_t)getpid(); + seed = (uint64_t)getpid(); #ifdef HAVE_GETTIMEOFDAY { @@ -210,14 +210,14 @@ long long G_random_generate_seed(void) if (gettimeofday(&tv, NULL) < 0) G_fatal_error(_("gettimeofday failed: %s"), strerror(errno)); - seed += (uint_least64_t)tv.tv_sec; - seed += (uint_least64_t)tv.tv_usec; + seed += (uint64_t)tv.tv_sec; + seed += (uint64_t)tv.tv_usec; } #else { time_t t = time(NULL); - seed += (uint_least64_t)t; + seed += (uint64_t)t; } #endif @@ -370,7 +370,7 @@ double G_drand48(void) /* Advance a generator of the program's own by one step. The * multiplication may wrap around at 2^64; that does not change the result * modulo 2^48. */ -static uint_least64_t lcg_step(uint_least64_t x) +static uint64_t lcg_step(uint64_t x) { return (LCG_A * x + LCG_B) & MASK48; } @@ -378,19 +378,19 @@ static uint_least64_t lcg_step(uint_least64_t x) /* Turn a seed into a generator state the way G_srand48() does: only the * low 32 bits are used, so a negative seed gives the state of its two's * complement 32-bit value. */ -static uint_least64_t lcg_seed(long long seed) +static uint64_t lcg_seed(long long seed) { - return (((uint_least64_t)seed & 0xFFFFFFFF) << 16) | 0x330E; + return (((uint64_t)seed & 0xFFFFFFFF) << 16) | 0x330E; } /* Advance the generator by an arbitrary number of steps without taking * them one at a time. One step is the affine map x -> a * x + c, and * composing two such maps gives another, so the map for `steps` steps is * built by repeated squaring, as an integer power would be. */ -static uint_least64_t lcg_jump(uint_least64_t x, uint_least64_t steps) +static uint64_t lcg_jump(uint64_t x, uint64_t steps) { - uint_least64_t a_total = 1, c_total = 0; /* the identity map */ - uint_least64_t a = LCG_A, c = LCG_B; /* one step */ + uint64_t a_total = 1, c_total = 0; /* the identity map */ + uint64_t a = LCG_A, c = LCG_B; /* one step */ while (steps) { if (steps & 1) { @@ -432,9 +432,9 @@ static void check_units(long long units) * different batches is as unrelated as the generator allows; an even * distance would relate those draws, the more simply the more often 2 * divides it. */ -static uint_least64_t batch_distance(const struct G_random_layout *layout) +static uint64_t batch_distance(const struct G_random_layout *layout) { - return ((uint_least64_t)layout->units * (uint_least64_t)layout->stride) | 1; + return ((uint64_t)layout->units * (uint64_t)layout->stride) | 1; } /* Fill in a layout of units of the given stride, all of one batch first, @@ -444,7 +444,7 @@ static uint_least64_t batch_distance(const struct G_random_layout *layout) static void fill_layout(struct G_random_layout *layout, long long seed, long long units, long long stride) { - uint_least64_t draws = (uint_least64_t)units * (uint_least64_t)stride; + uint64_t draws = (uint64_t)units * (uint64_t)stride; layout->start = lcg_seed(seed); layout->units = units; @@ -452,7 +452,7 @@ static void fill_layout(struct G_random_layout *layout, long long seed, layout->batches = 0; if (draws <= LCG_SPAN) layout->batches = - (long long)((LCG_SPAN - draws) / batch_distance(layout) + 1); + (int64_t)((LCG_SPAN - draws) / batch_distance(layout) + 1); layout->whole_span = 0; } @@ -607,7 +607,7 @@ void G_random_init_layout_bounded(struct G_random_layout *layout, void G_random_init_layout(struct G_random_layout *layout, long long seed, long long units) { - uint_least64_t parts, stride; + uint64_t parts, stride; check_seed(seed); check_units(units); @@ -616,13 +616,13 @@ void G_random_init_layout(struct G_random_layout *layout, long long seed, "most %lld units, not %lld; a bounded layout serves " "any number of units"), WHOLE_SPAN_MAX_UNITS, units); - parts = (uint_least64_t)units | 1; + parts = (uint64_t)units | 1; stride = LCG_SPAN / parts; if (parts > 1 && stride % 2 == 0) stride--; layout->start = lcg_seed(seed); layout->units = units; - layout->stride = (long long)stride; + layout->stride = (int64_t)stride; layout->batches = 1; layout->whole_span = 1; } @@ -684,12 +684,12 @@ void G_random_state_for_batch(struct G_random_state *state, const struct G_random_layout *layout, long long batch, long long unit) { - uint_least64_t offset; + uint64_t offset; if (unit < 0 || unit >= layout->units) G_fatal_error(_("Random number unit %lld is out of range " "(must be between 0 and %lld)"), - unit, layout->units - 1); + unit, (long long)(layout->units - 1)); if (batch < 0) G_fatal_error(_("Random number batch %lld is out of range " "(must not be negative)"), @@ -704,7 +704,7 @@ void G_random_state_for_batch(struct G_random_state *state, "Random number batch %lld is out of range " "(%lld batches fit into the generator's span)", (unsigned long)layout->batches), - batch, layout->batches); + batch, (long long)layout->batches); if (layout->batches == 0 && batch > 0) G_fatal_error(_("Random number batch %lld is out of range " "(no batch fits into the generator's span, only batch " @@ -714,8 +714,8 @@ void G_random_state_for_batch(struct G_random_state *state, /* With batch below the batches that fit, the offset is below 2^46; with * batch 0 of a layout which does not fit, it may reach past the span. * It is below 2^63 either way, since units times stride is. */ - offset = (uint_least64_t)batch * batch_distance(layout) + - (uint_least64_t)unit * (uint_least64_t)layout->stride; + offset = (uint64_t)batch * batch_distance(layout) + + (uint64_t)unit * (uint64_t)layout->stride; state->state = lcg_jump(layout->start, offset); } @@ -753,7 +753,7 @@ void G_random_advance(struct G_random_state *state, long long draws) G_fatal_error(_("Cannot advance a random number generator by %lld " "draws (the number must not be negative)"), draws); - state->state = lcg_jump(state->state, (uint_least64_t)draws); + state->state = lcg_jump(state->state, (uint64_t)draws); } /*! From 9bb7cf7564888f0c369a2b78efc5b133731006d5 Mon Sep 17 00:00:00 2001 From: Vaclav Petras Date: Fri, 2 Oct 2026 17:10:19 -0400 Subject: [PATCH 11/12] libgis: Make the whole-span flag of a random number layout a bool The field says whether a layout is one of the whole span, which holds a single batch and gets an error message of its own for a further batch. It is a truth value, for which gis.h asks for the C99 type bool, and it stood as an int beside the exact-width fields of the struct. gis.h now includes itself, as the raster and vector headers do. No source in the tree defines bool, true or false on its own. Written with the help of Claude Code. --- include/grass/gis.h | 3 ++- lib/gis/lrand48.c | 4 ++-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/include/grass/gis.h b/include/grass/gis.h index bc1a4c7e6c9..5a58dfc1b17 100644 --- a/include/grass/gis.h +++ b/include/grass/gis.h @@ -17,6 +17,7 @@ /*============================= Include Files ==============================*/ /* System include files */ +#include #include #include #include @@ -646,7 +647,7 @@ struct G_random_layout { int64_t units; int64_t stride; int64_t batches; /* batches that fit into the span, 0 when none */ - int whole_span; /* a layout of the whole span holds a single batch */ + bool whole_span; /* a layout of the whole span holds a single batch */ }; typedef int CELL; diff --git a/lib/gis/lrand48.c b/lib/gis/lrand48.c index 28298727ed4..9d7941d8cf4 100644 --- a/lib/gis/lrand48.c +++ b/lib/gis/lrand48.c @@ -453,7 +453,7 @@ static void fill_layout(struct G_random_layout *layout, long long seed, if (draws <= LCG_SPAN) layout->batches = (int64_t)((LCG_SPAN - draws) / batch_distance(layout) + 1); - layout->whole_span = 0; + layout->whole_span = false; } /*! @@ -624,7 +624,7 @@ void G_random_init_layout(struct G_random_layout *layout, long long seed, layout->units = units; layout->stride = (int64_t)stride; layout->batches = 1; - layout->whole_span = 1; + layout->whole_span = true; } /*! From fe5f83ffdcec6a87876cc2202e1b7cdeeebb0a73 Mon Sep 17 00:00:00 2001 From: Vaclav Petras Date: Fri, 2 Oct 2026 23:32:50 -0400 Subject: [PATCH 12/12] libgis: Take and return int64_t in the random number layout functions The functions took and returned long long while the layout struct holds int64_t, so the same values were named by two types at the boundary of the API. Use int64_t throughout, with INT64_MAX in the overflow checks and PRId64 in the messages, as the vector library does. Tools which hold a seed or a count as long long call the functions unchanged. Written with the help of Claude Code. --- include/grass/defs/gis.h | 25 ++++---- lib/gis/lrand48.c | 132 +++++++++++++++++++-------------------- 2 files changed, 77 insertions(+), 80 deletions(-) diff --git a/include/grass/defs/gis.h b/include/grass/defs/gis.h index 529354c9596..106c9801ef4 100644 --- a/include/grass/defs/gis.h +++ b/include/grass/defs/gis.h @@ -518,22 +518,21 @@ long G_srand48_auto(void); long G_lrand48(void); long G_mrand48(void); double G_drand48(void); -void G_random_state_from_seed(struct G_random_state *, long long); -void G_random_init_layout_exact(struct G_random_layout *, long long, long long, - long long); -void G_random_init_layout_bounded(struct G_random_layout *, long long, - long long, long long); -void G_random_init_layout(struct G_random_layout *, long long, long long); -long long G_random_layout_batches(const struct G_random_layout *); -long long G_random_layout_length(const struct G_random_layout *); +void G_random_state_from_seed(struct G_random_state *, int64_t); +void G_random_init_layout_exact(struct G_random_layout *, int64_t, int64_t, + int64_t); +void G_random_init_layout_bounded(struct G_random_layout *, int64_t, int64_t, + int64_t); +void G_random_init_layout(struct G_random_layout *, int64_t, int64_t); +int64_t G_random_layout_batches(const struct G_random_layout *); +int64_t G_random_layout_length(const struct G_random_layout *); void G_random_state_for_unit(struct G_random_state *, - const struct G_random_layout *, long long); + const struct G_random_layout *, int64_t); void G_random_state_for_batch(struct G_random_state *, - const struct G_random_layout *, long long, - long long); -void G_random_advance(struct G_random_state *, long long); + const struct G_random_layout *, int64_t, int64_t); +void G_random_advance(struct G_random_state *, int64_t); double G_random_double(struct G_random_state *); -long long G_random_generate_seed(void); +int64_t G_random_generate_seed(void); /* ls.c */ void G_set_ls_filter(int (*)(const char *, void *), void *); diff --git a/lib/gis/lrand48.c b/lib/gis/lrand48.c index 9d7941d8cf4..7c2f192a835 100644 --- a/lib/gis/lrand48.c +++ b/lib/gis/lrand48.c @@ -36,7 +36,7 @@ */ #include -#include +#include #include #include #include @@ -78,7 +78,7 @@ typedef signed int int32; /* The most units a layout of the whole span takes; see * G_random_init_layout() for why. */ -#define WHOLE_SPAN_MAX_UNITS (1LL << 20) +#define WHOLE_SPAN_MAX_UNITS (INT64_C(1) << 20) #if LRAND48_ATOMIC @@ -131,18 +131,18 @@ void G_srand48(long seedval) /* The seeds the generator tells apart: G_srand48() uses the low 32 bits of * the seed, and this range names every such seed once, as a signed or an * unsigned 32-bit value. */ -static int seed_in_range(long long seed) +static int seed_in_range(int64_t seed) { - return seed >= -(1LL << 31) && seed <= (1LL << 32) - 1; + return seed >= -(INT64_C(1) << 31) && seed <= (INT64_C(1) << 32) - 1; } /* Read a seed from an environment variable. An unset or empty variable * gives no seed, so that the caller can try the next source. */ -static int seed_from_environment(const char *name, long long *seed) +static int seed_from_environment(const char *name, int64_t *seed) { const char *text = getenv(name); char *end; - long long value; + int64_t value; if (!text || !*text) return 0; @@ -157,9 +157,9 @@ static int seed_from_environment(const char *name, long long *seed) G_fatal_error(_("Random number seed %s from %s is not an integer"), text, name); if (!seed_in_range(value)) { - long long reduced = (long long)((uint64_t)value & 0xFFFFFFFF); + int64_t reduced = (int64_t)((uint64_t)value & 0xFFFFFFFF); - G_warning(_("Random number seed %s from %s is used as %lld, " + G_warning(_("Random number seed %s from %s is used as %" PRId64 ", " "its low 32 bits"), text, name, reduced); value = reduced; @@ -174,7 +174,7 @@ static int seed_from_environment(const char *name, long long *seed) * The seed is the value of the environment variable GRASS_RANDOM_SEED, or * of SOURCE_DATE_EPOCH when GRASS_RANDOM_SEED is not set or empty, and * otherwise a weak hash of the current time and process ID. A value from - * the environment must be a decimal integer within the range of long long, + * the environment must be a decimal integer within the range of int64_t, * with nothing after it (leading white space and a sign are allowed); * anything else is a fatal error naming the variable. A value from -2^31 * to 2^32 - 1 is returned as it is, and a value outside that range is @@ -193,9 +193,9 @@ static int seed_from_environment(const char *name, long long *seed) * * \return the seed */ -long long G_random_generate_seed(void) +int64_t G_random_generate_seed(void) { - long long given; + int64_t given; uint64_t seed; if (seed_from_environment("GRASS_RANDOM_SEED", &given) || @@ -221,7 +221,7 @@ long long G_random_generate_seed(void) } #endif - return (long long)(seed & 0xFFFFFFFF); + return (int64_t)(seed & 0xFFFFFFFF); } /*! @@ -242,7 +242,7 @@ long long G_random_generate_seed(void) */ long G_srand48_auto(void) { - long long seed = G_random_generate_seed(); + int64_t seed = G_random_generate_seed(); G_srand48((long)seed); return (long)seed; @@ -378,7 +378,7 @@ static uint64_t lcg_step(uint64_t x) /* Turn a seed into a generator state the way G_srand48() does: only the * low 32 bits are used, so a negative seed gives the state of its two's * complement 32-bit value. */ -static uint64_t lcg_seed(long long seed) +static uint64_t lcg_seed(int64_t seed) { return (((uint64_t)seed & 0xFFFFFFFF) << 16) | 0x330E; } @@ -410,19 +410,19 @@ static uint64_t lcg_jump(uint64_t x, uint64_t steps) * as signed or as unsigned 32-bit values are accepted, so -1 and * 4294967295 are both accepted and are the same seed; anything else is * rejected rather than silently losing its high bits. */ -static void check_seed(long long seed) +static void check_seed(int64_t seed) { if (!seed_in_range(seed)) - G_fatal_error(_("Random number seed %lld is outside the range from " - "-2147483648 to 4294967295 the generator can use"), + G_fatal_error(_("Random number seed %" PRId64 " is outside the range " + "from -2147483648 to 4294967295 the generator can use"), seed); } -static void check_units(long long units) +static void check_units(int64_t units) { if (units <= 0) G_fatal_error(_("The number of units of a random number layout must " - "be positive, not %lld"), + "be positive, not %" PRId64), units); } @@ -441,8 +441,8 @@ static uint64_t batch_distance(const struct G_random_layout *layout) * then those of the next batch. The batches that fit are those whose last * stream ends within the span. Callers must check beforehand that * units * stride does not overflow. */ -static void fill_layout(struct G_random_layout *layout, long long seed, - long long units, long long stride) +static void fill_layout(struct G_random_layout *layout, int64_t seed, + int64_t units, int64_t stride) { uint64_t draws = (uint64_t)units * (uint64_t)stride; @@ -474,7 +474,7 @@ static void fill_layout(struct G_random_layout *layout, long long seed, * \param[out] state generator state to seed * \param[in] seed seed, from -2^31 to 2^32 - 1 */ -void G_random_state_from_seed(struct G_random_state *state, long long seed) +void G_random_state_from_seed(struct G_random_state *state, int64_t seed) { check_seed(seed); state->state = lcg_seed(seed); @@ -499,7 +499,7 @@ void G_random_state_from_seed(struct G_random_state *state, long long seed) * gislib_random_streams. * * A seed outside -2^31 to 2^32 - 1, a number of units or of draws which - * is not positive, or a product of the two beyond the range of long long + * is not positive, or a product of the two beyond the range of int64_t * is a fatal error. * * \param[out] layout layout to initialize @@ -507,20 +507,20 @@ void G_random_state_from_seed(struct G_random_state *state, long long seed) * \param[in] units number of units of work, positive * \param[in] draws_per_unit number of values each unit draws, positive */ -void G_random_init_layout_exact(struct G_random_layout *layout, long long seed, - long long units, long long draws_per_unit) +void G_random_init_layout_exact(struct G_random_layout *layout, int64_t seed, + int64_t units, int64_t draws_per_unit) { check_seed(seed); check_units(units); if (draws_per_unit <= 0) G_fatal_error(_("The number of random numbers a unit draws must be " - "positive, not %lld"), + "positive, not %" PRId64), draws_per_unit); - if (units > LLONG_MAX / draws_per_unit) - G_fatal_error(_("A random number layout of %lld units of %lld draws " - "each is too large (the number of draws must not " - "exceed %lld)"), - units, draws_per_unit, LLONG_MAX); + if (units > INT64_MAX / draws_per_unit) + G_fatal_error(_("A random number layout of %" PRId64 " units of " + "%" PRId64 " draws each is too large (the number of " + "draws must not exceed %" PRId64 ")"), + units, draws_per_unit, INT64_MAX); fill_layout(layout, seed, units, draws_per_unit); } @@ -541,33 +541,32 @@ void G_random_init_layout_exact(struct G_random_layout *layout, long long seed, * * A seed outside -2^31 to 2^32 - 1, a number of units or a bound which is * not positive, or a product of the units and the stride beyond the range - * of long long is a fatal error. + * of int64_t is a fatal error. * * \param[out] layout layout to initialize * \param[in] seed seed, see G_random_state_from_seed() * \param[in] units number of units of work, positive * \param[in] max_draws most values any unit draws, positive */ -void G_random_init_layout_bounded(struct G_random_layout *layout, - long long seed, long long units, - long long max_draws) +void G_random_init_layout_bounded(struct G_random_layout *layout, int64_t seed, + int64_t units, int64_t max_draws) { - long long stride; + int64_t stride; check_seed(seed); check_units(units); if (max_draws <= 0) G_fatal_error(_("The most random numbers a unit draws must be " - "positive, not %lld"), + "positive, not %" PRId64), max_draws); - /* LLONG_MAX is odd, so an even bound can be rounded up. */ + /* INT64_MAX is odd, so an even bound can be rounded up. */ stride = max_draws % 2 == 0 ? max_draws + 1 : max_draws; - if (units > LLONG_MAX / stride) - G_fatal_error(_("A random number layout of %lld units of at most " - "%lld draws each is too large (the number of units " - "times the bound rounded up to odd must not exceed " - "%lld)"), - units, max_draws, LLONG_MAX); + if (units > INT64_MAX / stride) + G_fatal_error(_("A random number layout of %" PRId64 " units of at " + "most %" PRId64 " draws each is too large (the number " + "of units times the bound rounded up to odd must not " + "exceed %" PRId64 ")"), + units, max_draws, INT64_MAX); fill_layout(layout, seed, units, stride); } @@ -604,8 +603,8 @@ void G_random_init_layout_bounded(struct G_random_layout *layout, * \param[in] seed seed, see G_random_state_from_seed() * \param[in] units number of units of work, from 1 to 2^20 */ -void G_random_init_layout(struct G_random_layout *layout, long long seed, - long long units) +void G_random_init_layout(struct G_random_layout *layout, int64_t seed, + int64_t units) { uint64_t parts, stride; @@ -613,8 +612,8 @@ void G_random_init_layout(struct G_random_layout *layout, long long seed, check_units(units); if (units > WHOLE_SPAN_MAX_UNITS) G_fatal_error(_("A random number layout of the whole span takes at " - "most %lld units, not %lld; a bounded layout serves " - "any number of units"), + "most %" PRId64 " units, not %" PRId64 "; a bounded " + "layout serves any number of units"), WHOLE_SPAN_MAX_UNITS, units); parts = (uint64_t)units | 1; stride = LCG_SPAN / parts; @@ -635,7 +634,7 @@ void G_random_init_layout(struct G_random_layout *layout, long long seed, * \return the batches that fit, 1 for a layout of the whole span, 0 when * one batch of the layout does not fit into the span */ -long long G_random_layout_batches(const struct G_random_layout *layout) +int64_t G_random_layout_batches(const struct G_random_layout *layout) { return layout->batches; } @@ -648,7 +647,7 @@ long long G_random_layout_batches(const struct G_random_layout *layout) * \return the stride of the layout, the number of values every unit may * draw without running into the next unit's stream */ -long long G_random_layout_length(const struct G_random_layout *layout) +int64_t G_random_layout_length(const struct G_random_layout *layout) { return layout->stride; } @@ -682,31 +681,31 @@ long long G_random_layout_length(const struct G_random_layout *layout) */ void G_random_state_for_batch(struct G_random_state *state, const struct G_random_layout *layout, - long long batch, long long unit) + int64_t batch, int64_t unit) { uint64_t offset; if (unit < 0 || unit >= layout->units) - G_fatal_error(_("Random number unit %lld is out of range " - "(must be between 0 and %lld)"), - unit, (long long)(layout->units - 1)); + G_fatal_error(_("Random number unit %" PRId64 " is out of range " + "(must be between 0 and %" PRId64 ")"), + unit, layout->units - 1); if (batch < 0) - G_fatal_error(_("Random number batch %lld is out of range " + G_fatal_error(_("Random number batch %" PRId64 " is out of range " "(must not be negative)"), batch); if (batch > 0 && layout->whole_span) - G_fatal_error(_("Random number batch %lld is out of range " + G_fatal_error(_("Random number batch %" PRId64 " is out of range " "(a layout of the whole span holds a single batch)"), batch); if (layout->batches >= 1 && batch >= layout->batches) - G_fatal_error(n_("Random number batch %lld is out of range " - "(%lld batch fits into the generator's span)", - "Random number batch %lld is out of range " - "(%lld batches fit into the generator's span)", + G_fatal_error(n_("Random number batch %" PRId64 " is out of range " + "(%" PRId64 " batch fits into the generator's span)", + "Random number batch %" PRId64 " is out of range " + "(%" PRId64 " batches fit into the generator's span)", (unsigned long)layout->batches), - batch, (long long)layout->batches); + batch, layout->batches); if (layout->batches == 0 && batch > 0) - G_fatal_error(_("Random number batch %lld is out of range " + G_fatal_error(_("Random number batch %" PRId64 " is out of range " "(no batch fits into the generator's span, only batch " "0 can be used)"), batch); @@ -729,8 +728,7 @@ void G_random_state_for_batch(struct G_random_state *state, * \param[in] unit number of the unit, from 0 to units - 1 */ void G_random_state_for_unit(struct G_random_state *state, - const struct G_random_layout *layout, - long long unit) + const struct G_random_layout *layout, int64_t unit) { G_random_state_for_batch(state, layout, 0, unit); } @@ -747,11 +745,11 @@ void G_random_state_for_unit(struct G_random_state *state, * \param[in,out] state a seeded generator state * \param[in] draws number of values to skip, not negative */ -void G_random_advance(struct G_random_state *state, long long draws) +void G_random_advance(struct G_random_state *state, int64_t draws) { if (draws < 0) - G_fatal_error(_("Cannot advance a random number generator by %lld " - "draws (the number must not be negative)"), + G_fatal_error(_("Cannot advance a random number generator by " + "%" PRId64 " draws (the number must not be negative)"), draws); state->state = lcg_jump(state->state, (uint64_t)draws); }