Eratosthenes 的筛子有巨大的“透支”——到底是 Sundaram 更好吗？

Question

Eratosthenes 的标准筛子多次删除了大多数复合材料；事实上，唯一不会被多次标记的是那些恰好是两个素数的乘积。自然，随着筛子变大，透支也会增加。

对于奇数筛（即没有偶数），透支达到 100%，n = 3,509,227，有 1,503,868 个复合数和 1,503,868 个已划掉数字的划线。对于 n = 2^32，透支上升到 134.25%（透支 2,610,022,328 与弹出计数 1,944,203,427 = (2^32 / 2) - 203,280,221）。

Sieve of Sundaram - 在maths.org 有另一种解释 - 可能会更聪明一点，如果 - 且仅当 - 智能计算循环限制。然而，我所看到的消息来源似乎掩盖了这一点“优化”，而且似乎未优化的 Sundaram 每次都会被只有赔率的 Eratosthenes 击败。

有趣的是，两者都创建了完全相同的最终位图，即位 k 对应于数字 (2 * k + 1) 的位图。所以两种算法最终都必须设置完全相同的位，它们只是有不同的处理方式。

有人亲身体验过具有竞争力的、经过调整的 Sundaram 吗？能打败古希腊吗？

我已经精简了我的小因子筛子的代码（2^32，一个只有赔率的希腊语），并将片段大小调整为 256 KB，这对于具有 256 KB L2 的旧 Nehalem 和新版 Nehalem 来说是最佳的CPU（即使后者对更大的段更宽容）。但是现在我撞到了一堵砖墙，血筛仍然需要 8.5 秒来初始化。从硬盘加载筛子不是一个很有吸引力的选择，而且多线程很难以可移植的方式进行（因为像 boost 这样的库往往会给可移植性带来麻烦）......

Sundaram 能否将启动时间缩短几秒钟？

P.S.：透支本身不是问题，会被 L2 缓存吸收。关键是标准的埃拉托色尼似乎比必要的工作量增加了一倍以上，这表明可能有可能做的工作更少，速度更快。

Answer 1

由于“Sundaram vs. Eratosthenes”的问题没有任何接受者，所以我坐下来分析了它。结果：经典的 Sundaram's 的透支率高于仅赔率的 Eratosthenes；如果你应用一个明显的、小的优化，那么透支是完全一样的——原因很明显。如果您修复 Sundaram 以完全避免过度绘制，那么您会得到类似 Pritchard's Sieve 的东西，这要复杂得多。

exposition of Sundaram's Sieve in Lucky's Notes 可能是迄今为止最好的；稍微重写以使用假设的（即非标准且未在此处提供）类型bitmap_t 它看起来有点像这样。为了测量过度绘制，位图类型需要对应于BTS（位测试和设置）CPU 指令的操作，该指令可通过 Wintel 编译器和 gcc 的 MinGW 版本的 _bittestandset() 内部函数获得。内在函数对性能非常不利，但对于计算过度绘制非常方便。

注意：对于筛选所有素数最多为 N 的人将调用具有 max_bit = N/2 的筛子；如果结果位图的位 i 被设置，则数字 (2 * i + 1) 是复合数。 该函数的名称中包含“31”，因为对于大于 2^31 的位图，索引数学会中断；因此这段代码只能筛选最多 2^32-1 的数字（对应于 max_bit

uint64_t Sundaram31_a (bitmap_t &bm, uint32_t max_bit)
{
   assert( max_bit <= UINT32_MAX / 2 );

   uint32_t m = max_bit;
   uint64_t overdraw = 0;

   bm.set_all(0);

   for (uint32_t i = 1; i < m / 2; ++i)
   {
      for (uint32_t j = i; j <= (m - i) / (2 * i + 1); ++j)
      {
         uint32_t k = i + j + 2 * i * j;

         overdraw += bm.bts(k);
      }
   }

   return overdraw;
}

Lucky 对j 的限制是准确的，但对i 的限制非常宽松。收紧它并丢失我添加的m 别名，以使代码看起来更像网络上的常见说明，我们得到：

uint64_t Sundaram31_b (bitmap_t &bm, uint32_t max_bit)
{
   uint32_t i_max = uint32_t(std::sqrt(double(2 * max_bit + 1)) - 1) / 2;
   uint64_t overdraw = 0;

   bm.set_all(0);

   for (uint32_t i = 1; i <= i_max; ++i)
   {
      for (uint32_t j = i; j <= (max_bit - i) / (2 * i + 1); ++j)
      {
         uint32_t k = i + j + 2 * i * j;

         overdraw += bm.bts(k);
      }
   }

   return overdraw;
}

assert 被丢弃是为了减少噪音，但它实际上仍然有效且必要。现在是时候进行一点强度降低，将乘法转换为迭代加法：

uint64_t Sundaram31_c (bitmap_t &bm, uint32_t max_bit)
{
   uint32_t i_max = uint32_t(std::sqrt(double(2 * max_bit + 1)) - 1) / 2;
   uint64_t overdraw = 0;

   bm.set_all(0);

   for (uint32_t i = 1; i <= i_max; ++i)
   {
      uint32_t n = 2 * i + 1;
      uint32_t k = n * i + i;   // <= max_bit because that's how we computed i_max
      uint32_t j_max = (max_bit - i) / n;

      for (uint32_t j = i; j <= j_max; ++j, k += n)
      {
         overdraw += bm.bts(k);
      }
   }

   return overdraw;
}

将循环条件转换为使用k可以让我们失去j；事情现在应该看起来非常熟悉......

uint64_t Sundaram31_d (bitmap_t &bm, uint32_t max_bit)
{
   uint32_t i_max = uint32_t(std::sqrt(double(2 * max_bit + 1)) - 1) / 2;
   uint64_t overdraw = 0;

   bm.set_all(0);

   for (uint32_t i = 1; i <= i_max; ++i)
   {
      uint32_t n = 2 * i + 1;     
      uint32_t k = n * i + i;   // <= max_bit because that's how we computed i_max

      for ( ; k <= max_bit; k += n)
      {        
         overdraw += bm.bts(k);
      }
   }

   return overdraw;
}

随着事情的发展，是时候分析一下数学上是否有理由证明某个明显的小变化。证明留给读者作为练习......

uint64_t Sundaram31_e (bitmap_t &bm, uint32_t max_bit)
{
   uint32_t i_max = unsigned(std::sqrt(double(2 * max_bit + 1)) - 1) / 2;
   uint64_t overdraw = 0;

   bm.set_all(0);

   for (uint32_t i = 1; i <= i_max; ++i)
   {
      if (bm.bt(i))  continue;

      uint32_t n = 2 * i + 1;     
      uint32_t k = n * i + i;   // <= m because we computed i_max to get this bound

      for ( ; k <= max_bit; k += n)
      {        
         overdraw += bm.bts(k);
      }
   }

   return overdraw;
}

与经典的仅赔率 Eratosthenes（除了名称）唯一不同的是k 的初始值，对于古希腊语，通常为(n * n) / 2。但是，将2 * i + 1 替换为n 差异结果为 1/2，四舍五入为 0。因此，Sundaram 的 仅赔率 Eratosthenes 没有跳过复合材料的“优化”以避免至少有一些已经划掉的数字划掉了。 i_max 的值与希腊语的 max_factor_bit 相同，只是使用完全不同的逻辑步骤得出并使用略有不同的公式计算得出。

PS：在代码中看到overdraw 这么多次之后，人们可能会想知道它实际上是什么...筛选到 2^32-1 的数字（即完整的 2^31 位图） Sundaram 的透支为 8,643,678,027（大约 2 * 2^32）或 444.6%；通过将其变为仅赔率的小修复 Eratosthenes，透支变为 2,610,022,328 或 134.2%。