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í como hi y lo .

 #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 ³¹ -1), entonces, por definición, la reducción no ha ocurrido y es necesario un segundo paso. Restar 2 ³¹ -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); }