埃拉托斯特尼筛法有巨大的“透支”——到底桑达拉姆筛法更好吗？

Question

标准埃拉托斯特尼筛法多次剔除大多数复合材料；事实上，唯一不会被标记多次的就是那些恰好是两个素数的乘积。当然，随着筛子变大，透支也会增加。

对于奇数筛（即没有偶数），透支达到 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）。

当且仅当循环限制是智能计算的时， Sundaram 筛法（在maths.org上有另一种解释）可能会更聪明一些。然而，我所看到的消息来源似乎将其掩盖为“优化”，而且似乎未优化的 Sundaram 每次都会被仅赔率的埃拉托色尼击败。

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

有人拥有竞技性、调校 Sundaram 的实践经验吗？它能打败古希腊语吗？

我精简了小因子筛的代码（2^32，仅赔率的希腊语），并将段大小调整为 256 KB，这对于具有 256 KB L2 的旧 Nehalem 和较新的 CPU 来说都是最佳的（甚至尽管后者对更大的细分市场更加宽容）。但现在我已经碰壁了，该死的筛子仍然需要 8.5 秒来初始化。从硬盘加载筛子并不是一个非常有吸引力的选择，并且多线程很难以可移植的方式完成（因为像 boost 这样的库往往会增加可移植性）...

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

PS：透支本身不是问题，会被二级缓存吸收。关键是，标准的埃拉托斯特尼似乎做了比必要的工作多一倍的工作，这表明可能有可能更快地完成更少的工作。

Answer 1

由于没有人接受“桑达拉姆与埃拉托色尼”问题，我坐下来分析了它。结果：经典的 Sundaram 的透支率严格高于仅赔率的埃拉托斯特尼 (Eratosthenes)；如果你应用一个明显的、小的优化，那么透支是完全相同的——原因将变得显而易见。如果你修复 Sundaram 来完全避免透支，那么你会得到类似Pritchard 筛子的东西，它要复杂得多。

《幸运笔记》中对桑达拉姆筛子的阐述可能是迄今为止最好的；稍微重写以使用假设的（即非标准且此处未提供）类型，bitmap_t它看起来有点像这样。为了测量过度绘制，位图类型需要与BTS（位测试和设置）CPU 指令相对应的操作，该指令可通过 Wintel 编译器和 MinGW 版本的 gcc 的 _bittestandset() 内在函数获得。内在函数对性能非常不利，但对于计算透支非常方便。

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

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

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

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