I am trying to implement a random-number generator with Mersenne prime (2^31-1) as the modulus. Following working code was based on several related posts:
- How do I extract specific 'n' bits of a 32-bit unsigned integer in C?
- Fast multiplication and subtraction modulo a prime
- Fast multiplication modulo 2^16 + 1
However,
It does not work with unint32_t hi, lo;
. Which means, I do not understand signed vs. unsigned aspect of the problem.
Based on #2 above, I was expecting the answer to be (hi+lo). Which means, I do not understand why the following statement is needed.
if (x1 > r)
x1 += r + 2;
-
Can someone please clarify the source of my confusion?
-
Can the code itself be improved?
-
Should the generator avoid 0 or 2^31-1 as a seed?
-
How would the code change for a prime (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);
}
Aucun commentaire:
Enregistrer un commentaire