I am trying to implement a random-number generator with Mersenne prime (2^{31}-1) as the modulus. The 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 `uint32_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);
}
```

·
Santiago Trujillo

The following if statement

```
if (x1 > r)
x1 += r + 2;
```

should be written as

```
if (x1 > r)
x1 -= r;
```

Both results are the same modulo 2^31:

```
x1 + r + 2 = x1 + 2^31 - 1 + 2 = x1 + 2^31 + 1
x1 - r = x1 - (2^31 - 1) = x1 - 2^31 + 1
```

The first solution overflows an `int32_t`

and assumes that conversion from `uint64_t`

to `int32_t`

is modulo 2^31. While many C compilers handle the conversion this way, this is not mandated by the C standard. The actual result is implementation-defined.

The second solution avoids the overflow and works with both `int32_t`

and `uint32_t`

.

You can also use an integer constant for `r`

:

```
uint64_t r = 0x7FFFFFFF; // 2^31 - 1
```

Or simply

```
uint64_t r = INT32_MAX;
```

**EDIT:** For primes of the form 2^p-k, you have to use masks with p bits and calculate the result with

```
uint32_t x1 = (k * hi + lo) % ((1 << p) - k)
```

If `k * hi + lo`

can overflow a `uint32_t`

(that is `(k + 1) * (2^p - 1) >= 2^32`

), you have to use 64-bit arithmetic:

```
uint32_t x1 = ((uint64_t)a * x) % ((1 << p) - k)
```

Depending on the platform, the latter might be faster anyway.

·
Santiago Trujillo
Report

Sue provided this as a solution:

With some experimentation (new code at the bottom), I was able to use

`uint32_t`

, which further suggests that I do not understand how the signed integers work with bit operations.The following code uses

`uint32_t`

for input as well as`hi`

and`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); }`

The issue was that my assumption that the reduction will be complete after one pass. If (x > 2

^{31}-1), then by definition the reduction has not occurred and a second pass is necessary. Subtracting 2^{31}-1, in that case does the trick. In the second attempt above,`r = 2^31-1`

and is therefore the modulus.`x -= r`

achieves the final reduction.Perhaps someone with expertise in random numbers or modular reduction could explain it better.

Cleaned function without

`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); }`

·
Santiago Trujillo
Report