Fast modular multiplication modulo prime for linear congruential generator in C

Question

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

1. How do I extract specific 'n' bits of a 32-bit unsigned integer in C?
2. Fast multiplication and subtraction modulo a prime
3. 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³¹-1 as a seed?

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

Original code

#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

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.

Answer 2

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³¹-1), then by definition the reduction has not occurred and a second pass is necessary. Subtracting 2³¹-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);
}