如何有效地筛选素数的选定范围？

Question

我一直在解决Project Euler 和 Sphere Online Judge 问题。在这个特定的问题中，我必须找到两个给定数字中的所有素数。我有一个看起来很有前途的功能（基于Sieve of Eratosthenes），但它太慢了。有人能发现是什么让我的功能如此缓慢，并暗示我该如何解决吗？此外，一些关于如何进行一般优化的 cmets（或指向此类 cmets/书籍/文章等的链接）将不胜感激。

代码：

def ranged_sieve(l, b)
  primes = (l..b).to_a
  primes[0]=nil if primes[0] < 2
  (2..Math.sqrt(b).to_i).each do |counter|
    step_from = l / counter
    step_from = step_from * counter
    l > 3 ? j = step_from : j = counter + counter
    (j..b).step(counter) do |stepped|
      index = primes.index(stepped)
      primes[index] = nil if index
    end
  end
  primes.compact
end

Answer 1

我还没有完全看清楚，但一个因素是，您将primes 中的某个值替换为nil，然后compact-ing 将其删除。这是一种浪费。只需直接使用delete_at 执行此操作，速度就会提高两倍以上：

def ranged_sieve2(l, b)
  primes = (l..b).to_a
  primes.delete_at(0) if primes[0] < 2
  (2..Math.sqrt(b).to_i).each do |counter|
    step_from = l / counter
    step_from = step_from * counter
    l > 3 ? j = step_from : j = counter + counter
    (j..b).step(counter) do |stepped|
      index = primes.index(stepped)
      primes.delete_at(index) if index
    end
  end
  primes
end

ranged_sieve(1, 100) # => Took approx 8e-4 seconds on my computer
ranged_sieve2(1, 100) # => Took approx 3e-4 seconds on my computer

另一个需要改进的地方是，当相关大小变大时，使用散列比数组快得多。用哈希替换你的数组实现，你可以得到这个：

def ranged_sieve3(l, b)
  primes = (l..b).inject({}){|h, i| h[i] = true; h}
  primes.delete(0)
  primes.delete(1)
  (2..Math.sqrt(b).to_i).each do |counter|
    step_from = l / counter
    step_from = step_from * counter
    l > 3 ? j = step_from : j = counter + counter
    (j..b).step(counter) do |stepped|
      primes.delete(stepped)
    end
  end
  primes.keys
end

当你这样做range_sieve3(1, 100) 时，由于开销，它比range_sieve2(1, 100) 慢。但是当你把数字变大时，优势就会变得突出。例如，我在我的电脑上得到了这个结果：

ranged_sieve(1, 1000) # => Took 1e-01 secs
ranged_sieve2(1, 1000) # => Took 3e-02 secs
ranged_sieve3(1, 1000) # => Took 8e-04 secs

Answer 2

SPOJ（Sphere Online Judges）的PRIME1 problem 旨在让你不能简单地筛选到上限，因为在这种情况下你会被超时击中。

一种可能的方法是提高速度；通过在标准筛子上添加一些花里胡哨的东西，它可以运行得足够快，以保持远低于超时限制。速度优化包括：

• 仅表示筛子中的奇数（节省 50% 的空间）
• 筛选适合缓存的小段，适合 L1 缓存 (32 KB)
• 通过小素数进行预测（即在筛段上爆破预先计算的模式）
• 记住跨段的每个素数的上一个（或下一个）工作偏移，而不是使用慢除法重新计算它们

将所有这些放在一起将筛选整个 2^32 范围的时间从 20 秒缩短到 2 秒，远低于 SPOI 超时。 My pastebin 有三个简单的 C++ 演示程序，您可以在其中检查每个优化的实际效果并查看它们的效果。

一个更简单的方法是只做必要的工作：筛选到目标范围的最后一个数字的平方根以获得所有潜在的质因子，然后只筛选目标范围本身。这样，您可以在不到两打代码和几毫秒的运行时间中解决 SPOJ 问题。我刚刚完成了一个demo .cpp for this type of segmented sieving（困难的部分不是筛子而是舒适测试的测试框架，由于几乎没有任何参考数据，因此验证正确操作高达2^64-1）。

筛子本身是这样的；筛子是仅赔率的打包位图，并且筛子范围以位为单位指定以确保稳健性（所有内容都在 .cpp 中进行了解释），因此您将通过 (range_start / 2) 获取offset：

unsigned char odd_composites32[UINT32_MAX / (2 * CHAR_BIT) + 1];   // the small factor sieve
uintxx_t sieved_bits = 0;                                          // how far it's been initialised

void extend_factor_sieve_to_cover (uintxx_t max_factor_bit);       // bit, not number!

void sieve32 (unsigned char *target_segment, uint64_t offset, uintxx_t bit_count)
{
   assert( bit_count > 0 && bit_count <= UINT32_MAX / 2 + 1 );

   uintxx_t max_bit = bit_count - 1;
   uint64_t max_num = 2 * (offset + max_bit) + 1;
   uintxx_t max_factor_bit = (max_factor32(max_num) - 1) / 2;

   if (target_segment != odd_composites32)
   {
      extend_factor_sieve_to_cover(max_factor_bit);
   }

   std::memset(target_segment, 0, std::size_t((max_bit + CHAR_BIT) / CHAR_BIT));

   for (uintxx_t i = 3u >> 1; i <= max_factor_bit; ++i)
   {
      if (bit(odd_composites32, i))  continue;

      uintxx_t n = (i << 1) + 1;   // the actual prime represented by bit i (< 2^32)

      uintxx_t stride = n;         // == (n * 2) / 2
      uint64_t start = (uint64_t(n) * n) >> 1;
      uintxx_t k;

      if (start >= offset)
      {
         k = uintxx_t(start - offset);
      }
      else // start < offset
      {
         uintxx_t before_the_segment = (offset - start) % stride;

         k = before_the_segment == 0 ? 0 : stride - before_the_segment;
      }

      while (k <= max_bit)
      {
         set_bit(target_segment, k);

         // k can wrap since strides go up to almost 2^32
         if ((k += stride) < stride)
         {
            break;
         }
      }
   }
}

对于 SPOJ 问题 - 小于 2^32 的数字 - 对于所有变量（即 uint32_t 而不是 uintxx_t 和 uint64_t），无符号整数就足够了，而且有些事情可以进一步简化。此外，对于这些小范围，您可以使用 sqrt() 而不是 max_factor32()。

在演示代码中，extend_factor_sieve_to_cover() 在道德上与sieve32(odd_composites32, 0, max_factor_bit + 1) 等效，只是采用了缓存友好的小步骤。对于 SPOJ 问题，您可以简单地使用单个 sieve32() 调用，因为在小于 2^32 的数字中只有 6541 个小的奇质因数，您可以立即对其进行筛分。

因此，解决这个 SPOJ 问题的诀窍是使用分段筛分，它将总运行时间缩短到几毫秒。