为什么我的Atkin筛的实施比Eratosthenes慢？

Question

我正在用Ruby中的Project Euler做问题,并用Atkin的筛子找到素数,但它的运行速度比Eratosthenes筛的慢.问题是什么？

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

正如维基百科所说,"Atkin的现代筛选更复杂,但在适当优化时更快"(我的重点).

在第一组循环中节省一些时间的第一个显而易见的地方是在大于y时停止迭代.例如,如果是1,000,000并且是450,那么当大于435 时你应该停止迭代(而不是像你现在那样继续到999).所以你可以重写第一个循环:4*x**2 + y**2nnxy

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**20是模12,那么y**2必须是1或5模12,所以y必须是1,5,7,9或11模12.

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

Answer 2

如果您实际上首先正确实施了Eratosthenes的Sieve,那将会有很大帮助.

该筛子的关键特征是每次只能对一个数字进行一次操作.相比之下,你正在为每个素数工作而不是数量.差异很微妙,但性能影响巨大.

以下是您未能实施的实际筛选:

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问题都应该足够快.