为什么我的Atkin sieve实现比Eratosthenes慢？

Question

我正在做Ruby中的Project Euler中的问题，并实现了Atkin的sieve来寻找质数，但它比Eratosthenes的sieve运行得慢。有什么问题吗？

def atkin_sieve(n)
  primes = [2,3,5]

  sieve = Array.new(n+1, false)

  y_upper = n-4 > 0 ? Math.sqrt(n-4).truncate : 1
  for x in (1..Math.sqrt(n/4).truncate) 
    for y in (1..y_upper)
      k = 4*x**2 + y**2
      sieve[k] = !sieve[k] if k%12 == 1 or k%12 == 5
    end
  end

  y_upper = n-3 > 0 ? Math.sqrt(n-3).truncate : 1
  for x in (1..Math.sqrt(n/3).truncate)
    for y in (1..y_upper)
      k = 3*x**2 + y**2
      sieve[k] = !sieve[k] if k%12 == 7
    end
  end

  for x in (1..Math.sqrt(n).truncate)
    for y in (1..x)
      k = 3*x**2 - y**2
      if k < n and k%12 == 11
    sieve[k] = !sieve[k]
      end
    end
  end

  for j in (5...n)
    if sieve[j]
      prime = true
      for i in (0...primes.length)
        if j % (primes[i]**2) == 0
          prime = false
          break
        end
      end
      primes << j if prime
    end
  end
  primes
end

def erato_sieve(n)
  primes = []
  for i in (2..n)
    if primes.all?{|x| i % x != 0}
      primes << i
    end
  end
  primes
end

Answer 1

正如Wikipedia says所说，“阿特金的现代筛子更复杂，但如果适当优化，速度会更快”(我的重点)。

在第一组循环中节省一些时间的第一个明显的地方是，当4*x**2 + y**2大于n时，停止迭代y。例如，如果n是1,000,000，x是450，那么当y大于435时，您应该停止迭代(而不是像现在那样继续到999 )。因此您可以将第一个循环重写为：

for x in (1..Math.sqrt(n/4).truncate)
    X = 4 * x ** 2
    for y in (1..Math.sqrt(n - X).truncate)
        k = X + y ** 2
        sieve[k] = !sieve[k] if k%12 == 1 or k%12 == 5
    end
end

(这也避免了每次循环都重新计算4*x**2，尽管这可能是一个非常小的改进(如果有的话)。

当然，类似的说明也适用于y上的其他循环。

第二个可以提高速度的地方是在y上循环的策略。您循环遍历范围内y的所有值，然后检查哪些值导致k的值具有模数为12的正确余数。相反，您可以只循环遍历y的正确值，并避免完全测试余数。

如果4*x**2是4模12，则y**2必须是1或9模12，因此y必须是1、3、5、7或11模12。如果4*x**2是8模12，则y**2必须是5或9模12，因此y必须是3或9模12。最后，如果4*x**2是0模12，则y**2必须是1或5模12，因此y必须是1、5、7、9或11模12。

我还注意到你的Eratosthenes筛子通过测试所有低于i的素数的除性来做无用的工作。一旦你测试了所有小于或等于i平方根的素数的可除性，你就可以停止迭代。

Answer 2

如果你首先正确地实现了Eratosthenes的筛子，这将会有很大的帮助。

该筛子的关键特性是，每次素数除以一个数时，只需执行一次操作。相比之下，你是在为每个小于这个数字的素数做功。差异是微妙的，但对性能的影响是巨大的。

以下是您未能实现的实际筛选器：

def eratosthenes_primes(n)
    primes = []
    could_be_prime = (0..n).map{|i| true}
    could_be_prime[0] = false
    could_be_prime[1] = false
    i = 0
    while i*i <= n
        if could_be_prime[i]
            j = i*i
            while j <= n
                could_be_prime[j] = false
                j += i
            end
        end
        i += 1
    end
    return (2..n).find_all{|i| could_be_prime[i]}
end

将它与您的代码进行比较，以找到50,000个以内的所有素数。还要注意的是，通过特殊地封装偶数的逻辑，可以很容易地将速度提高2倍。通过这个调整，对于每个需要计算大量素数的Project Euler问题，这个算法应该足够快了。

Answer 3

@Gareth提到了一些关于4x^2+y^2的冗余计算。无论是在这里还是在其他你在循环中进行计算的地方，你都可以利用你已经执行的计算，并将其简化为简单的加法。

您可以依赖于X已经具有4 * (x-1) ** 2值这一事实，而不是X=4 * x ** 2。由于4x^2 = 4(x-1)^2 + 4(2x - 1)，您需要做的就是将8 * x - 4添加到X中。您可以对k和其他需要重复计算的地方(如3x^2 + y^2)使用相同的技巧。