diff --git a/include/grass/defs/gis.h b/include/grass/defs/gis.h index 498418fb01e..106c9801ef4 100644 --- a/include/grass/defs/gis.h +++ b/include/grass/defs/gis.h @@ -518,6 +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 *, 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 *, int64_t); +void G_random_state_for_batch(struct G_random_state *, + 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 *); +int64_t 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..5a58dfc1b17 100644 --- a/include/grass/gis.h +++ b/include/grass/gis.h @@ -17,6 +17,8 @@ /*============================= Include Files ==============================*/ /* System include files */ +#include +#include #include #include @@ -628,6 +630,26 @@ struct Popen { int pid; }; +/* A position of a pseudo-random number generator owned by the program; + * 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 { + 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 { + 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 */ + bool whole_span; /* a layout of the whole span holds a single batch */ +}; + typedef int CELL; typedef double DCELL; typedef float FCELL; diff --git a/lib/gis/lrand48.c b/lib/gis/lrand48.c index 1491dfe63a7..7c2f192a835 100644 --- a/lib/gis/lrand48.c +++ b/lib/gis/lrand48.c @@ -18,21 +18,33 @@ * * 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. 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 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. + * * 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 +#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 +66,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 layout of the whole span takes; see + * G_random_init_layout() for why. */ +#define WHOLE_SPAN_MAX_UNITS (INT64_C(1) << 20) + #if LRAND48_ATOMIC /* The whole 48-bit state is kept in one atomic integer so that it can be @@ -63,10 +89,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 +128,121 @@ 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(int64_t seed) +{ + 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, int64_t *seed) +{ + const char *text = getenv(name); + char *end; + int64_t 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 (!seed_in_range(value)) { + int64_t reduced = (int64_t)((uint64_t)value & 0xFFFFFFFF); + + G_warning(_("Random number seed %s from %s is used as %" PRId64 ", " + "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 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 + * reduced to its low 32 bits, between 0 and 2^32 - 1, with a warning. The + * 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. * - * 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) +int64_t G_random_generate_seed(void) { - unsigned long seed; - char *grass_random_seed = getenv("GRASS_RANDOM_SEED"); + int64_t given; + uint64_t 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 = (uint64_t)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 += (uint64_t)tv.tv_sec; + seed += (uint64_t)tv.tv_usec; + } #else - { - time_t t = time(NULL); + { + time_t t = time(NULL); - seed += (unsigned long)t; - } -#endif + seed += (uint64_t)t; } +#endif + + return (int64_t)(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_state_from_seed() read + * as the same seed + */ +long G_srand48_auto(void) +{ + int64_t seed = G_random_generate_seed(); G_srand48((long)seed); return (long)seed; @@ -272,6 +367,414 @@ 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 uint64_t lcg_step(uint64_t 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 uint64_t lcg_seed(int64_t seed) +{ + 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 uint64_t lcg_jump(uint64_t x, uint64_t steps) +{ + 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) { + 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(int64_t seed) +{ + if (!seed_in_range(seed)) + G_fatal_error(_("Random number seed %" PRId64 " is outside the range " + "from -2147483648 to 4294967295 the generator can use"), + seed); +} + +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 %" PRId64), + 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 uint64_t batch_distance(const struct G_random_layout *layout) +{ + 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, + * 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, int64_t seed, + int64_t units, int64_t stride) +{ + uint64_t draws = (uint64_t)units * (uint64_t)stride; + + layout->start = lcg_seed(seed); + layout->units = units; + layout->stride = stride; + layout->batches = 0; + if (draws <= LCG_SPAN) + layout->batches = + (int64_t)((LCG_SPAN - draws) / batch_distance(layout) + 1); + layout->whole_span = false; +} + +/*! + * \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, 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(). + * + * 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_state_from_seed(struct G_random_state *state, int64_t 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. 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, 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 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] draws_per_unit number of values each unit draws, positive + */ +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 %" PRId64), + draws_per_unit); + 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); +} + +/*! + * \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_batches() returns the batches that fit, + * 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 + * 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, int64_t seed, + int64_t units, int64_t max_draws) +{ + 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 %" PRId64), + max_draws); + /* INT64_MAX is odd, so an even bound can be rounded up. */ + stride = max_draws % 2 == 0 ? max_draws + 1 : max_draws; + 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); +} + +/*! + * \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 + * 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, 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 + * 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_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_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, int64_t seed, + int64_t units) +{ + uint64_t parts, stride; + + check_seed(seed); + check_units(units); + if (units > WHOLE_SPAN_MAX_UNITS) + G_fatal_error(_("A random number layout of the whole span takes at " + "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; + if (parts > 1 && stride % 2 == 0) + stride--; + layout->start = lcg_seed(seed); + layout->units = units; + layout->stride = (int64_t)stride; + layout->batches = 1; + layout->whole_span = true; +} + +/*! + * \brief Return the number of batches that fit into the span + * + * \param[in] layout an initialized layout + * + * \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 + */ +int64_t G_random_layout_batches(const struct G_random_layout *layout) +{ + return layout->batches; +} + +/*! + * \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 + */ +int64_t 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 batch + * + * A batch is one stream for every unit of the layout, and the batches + * 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 + * 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] batch number of the batch, from 0 + * \param[in] unit number of the unit, from 0 to units - 1 + */ +void G_random_state_for_batch(struct G_random_state *state, + const struct G_random_layout *layout, + int64_t batch, int64_t unit) +{ + uint64_t offset; + + if (unit < 0 || unit >= layout->units) + 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 %" PRId64 " is out of range " + "(must not be negative)"), + batch); + if (batch > 0 && layout->whole_span) + 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 %" 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, layout->batches); + if (layout->batches == 0 && batch > 0) + 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); + + /* 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 = (uint64_t)batch * batch_distance(layout) + + (uint64_t)unit * (uint64_t)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_batch() with batch 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, int64_t unit) +{ + G_random_state_for_batch(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, int64_t draws) +{ + if (draws < 0) + 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); +} + +/*! + * \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_state_from_seed(), G_random_state_for_unit() or + * G_random_state_for_batch() + * + * \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..baecdf6c453 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 @@ -6,11 +9,39 @@ 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_init_layout, + G_random_init_layout_bounded, + G_random_init_layout_exact, + G_random_layout_batches, + G_random_layout_length, + 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, +) # 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 +278,929 @@ 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 layout of the whole span takes. +WHOLE_SPAN_MAX_UNITS = 2**20 + +# 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. +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): + """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_state_from_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 whole_span_layout(seed, units): + layout = struct_G_random_layout() + G_random_init_layout(byref(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() + G_random_state_for_batch(byref(state), byref(layout), batch, unit) + return draw_states(state, n) + + +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. + """ + 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_state_from_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_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. + """ + G_srand48(seed) + shared = [G_drand48() for _ in range(100)] + state = struct_G_random_state() + 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"] + + +@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_state_from_seed {}", + "G_random_init_layout_exact {} 5 7", + "G_random_init_layout_bounded {} 5 7", + "G_random_init_layout {} 5", + ], + ids=["seed", "exact", "bounded", "whole_span"], +) +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_batches_and_length(units, draws): + """The batches that fit are those whose last stream ends within the span. + + 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)) == batches_that_fit(units, draws) + assert G_random_layout_length(byref(layout)) == draws + + +@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( + xy_session_for_module, tmp_path +): + """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_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_batch 1 0", + ) + assert_fatal(result, "batch 1 is out of range") + + +def test_exact_layout_batch_past_the_batches_that_fit_is_fatal( + xy_session_for_module, tmp_path +): + """A batch beyond those that fit is a fatal error naming both numbers.""" + units, draws = 4096, 2**34 + batches = batches_that_fit(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_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_batch {batches} 0", + ) + assert_fatal(result, f"batch {batches} is out of range ({batches} batch 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 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_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 or after the start of batch 1.""" + 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, batch=1) == reference_states( + seed, batch_distance(units, stride) + unit * stride, 3 + ), unit + + +@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 = whole_span_layout(1337, units) + stride = G_random_layout_length(byref(layout)) + assert stride == whole_span_stride(units) + assert stride % 2 == 1 + assert G_random_layout_batches(byref(layout)) == 1 + + +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 = whole_span_layout(1337, 1) + assert G_random_layout_length(byref(layout)) == SPAN + 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, 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 = 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_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 1337 4", + "G_random_state_for_batch 1 0", + ) + assert_fatal(result, "a layout of the whole span holds a single batch") + + +@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 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 1337 {units}" + ) + assert_fatal(result, f"at most {WHOLE_SPAN_MAX_UNITS} units, not {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 = whole_span_layout(1337, WHOLE_SPAN_MAX_UNITS) + assert G_random_layout_length(byref(layout)) == 67108799 + 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) + + +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 = whole_span_stride(units) + nearest = distance // stride + return min( + abs(k * stride - distance) / stride + for k in (nearest, nearest + 1) + if 1 <= k < units + ) + + +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 + 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_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, WHOLE_SPAN_MAX_UNITS + 1) + ) + closest_to_quarter = min( + 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 batch fitting, as (layout, units, stride) + + The stride is computed here from the documented rules, not read from the + layout. + """ + cases = [ + (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) + 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 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. + """ + batches = G_random_layout_batches(byref(layout)) + length = G_random_layout_length(byref(layout)) + assert length == stride + assert batches >= 1 + last_start = (batches - 1) * batch_distance(units, stride) + (units - 1) * stride + assert last_start + length <= SPAN + + +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 layout of the whole span start less + than 2^46 draws apart, so no lag relates their values by a constant. + """ + layout = whole_span_layout(1337, WHOLE_SPAN_MAX_UNITS) + first = unit_states(layout, 0, 100) + second = unit_states(layout, WHOLE_SPAN_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_state_from_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_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 = 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 + + +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_state_from_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_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) + + +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 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"), + ], +) +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_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_batch -1 0", + ) + assert_fatal(result, "batch -1 is out of range") + + +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) + assert draw_states(state, 10) == unit_states(layout, 7, 10) + + +GENERATE_SEED_SCRIPT = """ +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))) +""" + + +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, + ) + + +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 +): + """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()[:2] == [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()[:2] == ["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 = 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) + 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 = whole_span_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 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 = { + "whole_span": [ + 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"), + [ + ("whole_span", whole_span_layout(1337, 4096), 3 * 17175675903), + ("bounded", bounded_layout(1337, 1000, 100), 3 * 101), + ], + 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.""" + 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