6

我正在尝试实现一个以梅森素数 (2 31 -1) 作为模数的随机数生成器。以下工作代码基于几个相关帖子：

  1. 如何在 C 中提取 32 位无符号整数的特定“n”位？
  2. 以素数为模的快速乘法和减法
  3. 快速乘法模 2^16 + 1

然而，

它不适用于uint32_t hi, lo;，这意味着我不了解问题的有符号与无符号方面。

基于上面的＃2，我期待答案是（hi + lo）。这意味着，我不明白为什么需要以下语句。

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

  • 有人可以澄清我的困惑的根源吗？

  • 代码本身可以改进吗？

  • 生成器应该避免 0 或 2 31 -1 作为种子吗？

  • 素数（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);
}
4

2 回答 2

2

下面的 if 语句

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

应该写成

if (x1 > r)
    x1 -= r;

两个结果都是相同的模 2^31：

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

第一个解决方案溢出 anint32_t并假设从uint64_tto的转换int32_t是模 2^31。虽然许多 C 编译器以这种方式处理转换，但 C 标准并未强制要求这样做。实际结果是实现定义的。

第二种解决方案避免了溢出并与int32_t和一起使用uint32_t。

您还可以将整数常量用于r：

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

或者干脆

uint64_t r = INT32_MAX;

编辑：对于 2^pk 形式的素数，您必须使用带有 p 位的掩码并计算结果

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

如果k * hi + lo可以溢出 a uint32_t（即(k + 1) * (2^p - 1) >= 2^32），则必须使用 64 位算术：

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

根据平台的不同，后者可能会更快。

于 2015-06-25T11:07:32.323 回答
1

苏提供了这个作为解决方案：

通过一些实验（底部的新代码），我能够使用 uint32_t，这进一步表明我不了解有符号整数如何与位运算一起使用。

以下代码uint32_t用于输入以及hi和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);
  }

问题是我假设减少将在一次通过后完成。如果 (x > 2 31 -1)，则根据定义，还原没有发生，需要第二次通过。在这种情况下，减去 2 31 -1 就可以了。在上面的第二次尝试中，r = 2^31-1因此是模数。x -= r达到最终的还原。

也许具有随机数或模约简专业知识的人可以更好地解释它。

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);
}
于 2015-06-25T10:25:56.940 回答