mercredi 24 juin 2015

Fast modular multiplication modulo prime for linear congruential generator in C

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:

  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 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