Estoy tratando de implementar un generador de números aleatorios con Mersenne prime (2 31 -1) como módulo. El siguiente código de trabajo se basó en varias publicaciones relacionadas:
Sin embargo,
No funciona con uint32_t hi, lo;
, lo que significa que no entiendo el aspecto firmado frente al no firmado del problema.
Según el n. ° 2 anterior, esperaba que la respuesta fuera (hola + bajo). Lo que significa que no entiendo por qué se necesita la siguiente declaración.
if (x1 > r) x1 += r + 2;
¿Puede alguien aclarar el origen de mi confusión?
¿Se puede mejorar el código en sí?
¿Debe el generador evitar 0 o 2 31 -1 como semilla?
¿Cómo cambiaría el código para un número primo (2 p -k)?
#include <inttypes.h> // x1 = a*x0 (mod 2^31-1) int32_t lgc_m(int32_t a, int32_t x) { printf("x %"PRId32"\n", x); if (x == 2147483647){ printf("x1 %"PRId64"\n", 0); return (0); } uint64_t c, r = 1; c = (uint64_t)a * (uint64_t)x; if (c < 2147483647){ printf("x1 %"PRId64"\n", c); return (c); } int32_t hi=0, lo=0; int i, p = 31;//2^31-1 for (i = 1; i < p; ++i){ r |= 1 << i; } lo = (c & r) ; hi = (c & ~r) >> p; uint64_t x1 = (uint64_t ) (hi + lo); // NOT SURE ABOUT THE NEXT STATEMENT if (x1 > r) x1 += r + 2; printf("c %"PRId64"\n", c); printf("r %"PRId64"\n", r); printf("\tlo %"PRId32"\n", lo); printf("\thi %"PRId32"\n", hi); printf("x1 %"PRId64"\n", x1); printf("\n" ); return((int32_t) x1); } int main(void) { int32_t r; r = lgc_m(1583458089, 1); r = lgc_m(1583458089, 2000000000); r = lgc_m(1583458089, 2147483646); r = lgc_m(1583458089, 2147483647); return(0); }
La siguiente sentencia if
if (x1 > r) x1 += r + 2;
debe escribirse como
if (x1 > r) x1 -= r;
Ambos resultados son el mismo módulo 2^31:
x1 + r + 2 = x1 + 2^31 - 1 + 2 = x1 + 2^31 + 1 x1 - r = x1 - (2^31 - 1) = x1 - 2^31 + 1
La primera solución desborda un int32_t
y asume que la conversión de uint64_t
a int32_t
es módulo 2^31. Si bien muchos compiladores de C manejan la conversión de esta manera, el estándar C no lo exige. El resultado real está definido por la implementación.
La segunda solución evita el desbordamiento y funciona tanto con int32_t
como con uint32_t
.
También puede usar una constante entera para r
:
uint64_t r = 0x7FFFFFFF; // 2^31 - 1
O simplemente
uint64_t r = INT32_MAX;
EDITAR: para números primos de la forma 2^pk, debe usar máscaras con p bits y calcular el resultado con
uint32_t x1 = (k * hi + lo) % ((1 << p) - k)
Si k * hi + lo
puede desbordar un uint32_t
(es decir, (k + 1) * (2^p - 1) >= 2^32
), debe usar aritmética de 64 bits:
uint32_t x1 = ((uint64_t)a * x) % ((1 << p) - k)
Dependiendo de la plataforma, este último podría ser más rápido de todos modos.
Sue proporcionó esto como una solución:
Con un poco de experimentación (código nuevo en la parte inferior), pude usar
uint32_t
, lo que sugiere además que no entiendo cómo funcionan los enteros con signo con las operaciones de bits.El siguiente código usa
uint32_t
para la entrada, así comohi
ylo
.#include <inttypes.h> // x1 = a*x0 (mod 2^31-1) uint32_t lgc_m(uint32_t a, uint32_t x) { printf("x %"PRId32"\n", x); if (x == 2147483647){ printf("x1 %"PRId64"\n", 0); return (0); } uint64_t c, r = 1; c = (uint64_t)a * (uint64_t)x; if (c < 2147483647){ printf("x1 %"PRId64"\n", c); return (c); } uint32_t hi=0, lo=0; int i, p = 31;//2^31-1 for (i = 1; i < p; ++i){ r |= 1 << i; } hi = c >> p; lo = (c & r) ; uint64_t x1 = (uint64_t ) ((hi + lo) ); // NOT SURE ABOUT THE NEXT STATEMENT if (x1 > r){ printf("x1 - r = %"PRId64"\n", x1- r); x1 -= r; } printf("c %"PRId64"\n", c); printf("r %"PRId64"\n", r); printf("\tlo %"PRId32"\n", lo); printf("\thi %"PRId32"\n", hi); printf("x1 %"PRId64"\n", x1); printf("\n" ); return((uint32_t) x1); } int main(void) { uint32_t r; r = lgc_m(1583458089, 1583458089); r = lgc_m(1583458089, 2147483645); return(0); }
El problema era que mi suposición de que la reducción se completará después de un pase. Si (x > 2 31 -1), entonces, por definición, la reducción no ha ocurrido y es necesario un segundo paso. Restar 2 31 -1, en ese caso funciona. En el segundo intento anterior,
r = 2^31-1
y por lo tanto es el módulo.x -= r
logra la reducción final.Quizás alguien con experiencia en números aleatorios o reducción modular podría explicarlo mejor.
Función limpia sin
printf()
s.uint32_t lgc_m(uint32_t a, uint32_t x){ uint64_t c, x1, m = 2147483647; //modulus: m = 2^31-1 if (x == m) return (0); c = (uint64_t)a * (uint64_t)x; if (c < m)//no reduction necessary return (c); uint32_t hi, lo, p = 31;//2^p-1, p = 31 hi = c >> p; lo = c & m; x1 = (uint64_t)(hi + lo); if (x1 > m){//one more pass needed //this block can be replaced by x1 -= m; hi = x1 >> p; lo = (x1 & m); x1 = (uint64_t)(hi + lo); } return((uint32_t) x1); }