Fast modular multiplication modulo simple for linear congruent generator in C

Question

I am trying to implement a random number generator with a simple Mersenne (2 ³¹ -1) as a module. The following working code was based on several related posts:

• How to extract specific "n" bits of a 32-bit unsigned integer in C?
• Fast multiplication and subtraction modulo
• Fast multiplication modulo 2 ^ 16 + 1

but

It does not work with uint32_t hi, lo; , which means that I do not understand the signed or unsigned aspect of the problem.

Based on No. 2 above, I expected an answer (hi + lo). This means that I do not understand why the following statement is required.

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

• Can someone clarify the source of my confusion?

• Is it possible to improve the code?

• Should the generator avoid 0 or 2 ³¹ -1 as a seed?

• How would the code change for a simple (2 ^p -k)?

Source

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

Answer 1

Next if statement

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

should be written as

 if (x1 > r) x1 -= r; 

Both results coincide 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 int32_t and assumes that the conversion from uint64_t to int32_t is modulo 2 ^ 31. Although many C compilers handle the conversion in this way, this is not provided by the C standard. The actual result is determined by the implementation.

The second solution avoids overflow and works with both int32_t and uint32_t .

You can also use integer constant for r :

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

Or simply

 uint64_t r = INT32_MAX; 

EDIT: For prime numbers of the form 2 ^ pk you need 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 uint32_t (i.e. (k + 1) * (2^p - 1) >= 2^32 ), you should use 64-bit arithmetic:

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

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

Answer 2

Sue provided this as a solution:

With some experiments (new code below), I was able to use uint32_t , which also indicates that I don’t understand how 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 problem was that my assumption was that the cut would be completed after one pass. If (x> 2 ³¹ -1), then by definition the reduction did not occur, and a second pass is necessary. Subtracting 2 ³¹ -1, in this case the trick. In the second attempt above, r = 2^31-1 and therefore the modulus. x -= r reaches the final contraction.

Perhaps someone with experience in random numbers or modular reduction can explain this better.

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