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

Faunie

Look at the assembler output of the sgined operations. You will generally see Logical Shifts to the Left, but Arithmetic Shifts to the Right with signed values. Other bitwise operators should behave the same. The ASL operator will preserve the sign bits when shifting downwards. – Bosson

Thanks. The first attempt was signed and therefore added to the confusion. The second solution used unsigned and revealed the real problem. – Faunie

uint64_t x1 = (uint64_t ) ((hi + lo) ); Hmmm maybe uint64_t x1 = (uint64_t ) hi + lo;? – Gnosticize

Minor: r |= 1 << i; --> r |= (uint32_t)1 << i; (Cope with 16-bit int. and signed overflow.) – Gnosticize

@chux. Thanks for catching these errors! – Faunie

r |=1 part just creates the modulus r = m = 2147483647. Hard coding this eliminates the issue. – Faunie

It would be nicer to remove the solution from the question, post it as an answer (what it is !) and accept it. Currently, the question appears as a not answered one. – Outgroup

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.

Voorhees

Thanks. Even though, the code 'worked' I did not really understand what was going on. Your answer addresses all points of confusion. – Faunie

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

Sequestration