生成一个不超过一定数量的素数列表答案

【问题标题】：Generate a list of primes up to a certain number生成一个不超过一定数量的素数列表
【发布时间】：2011-04-16 22:34:05
【问题描述】：

我正在尝试生成一个低于 10 亿的素数列表。我正在尝试这个，但这种结构非常糟糕。有什么建议吗？

a <- 1:1000000000
d <- 0
b <- for (i in a) {for (j in 1:i) {if (i %% j !=0) {d <- c(d,i)}}}

【问题讨论】：

你需要一次所有的素数吗？如果需要，您似乎可以从这样的地方下载或获取它们：primes.utm.edu/lists/small/millions
这里可以得到素数列表prime-numbers.org
如果您只是对计算进行矢量化，由于 R 的性质，它们的运行速度会比筛子快得多。（请参阅我的答案）

标签： algorithm r primes

【解决方案1】：

没有建议，但请允许我进行一些扩展评论。我用以下代码运行了这个实验：

get_primes <- function(n_min, n_max){
  options(scipen=999)
    result = vector()
      for (x in seq(max(n_min,2), n_max)){
        has_factor <- F
        for (p in seq(2, ceiling(sqrt(x)))){
          if(x %% p == 0) has_factor <- T
          if(has_factor == T) break
          }
        if(has_factor==F) result <- c(result,x)
        }
    result
}

经过近 24 小时不间断的计算机操作，我得到了一个 5,245,897 素数列表。 π(1,000,000,000) = 50,847,534，因此完成此计算需要 10 天。

这些前~ 5 百万素数中的Here is the file。

【讨论】：

为了获得大约 500 万个素数，您的试运行上限是多少？我猜是一亿，但我不清楚。其次，您的pi(100,000,000) 语句中似乎缺少一个零（应该是pi(1,000,000,000)）。
@JosephWood 是的，就像在 OP 中一样。感谢您指出 pi 转录错误。
只是好奇，你看到我上面的回答了吗？它表明，你不仅可以在不到一秒的时间内生成 10 亿以下的素数，而且可以在大约 5 秒内得到 100 亿以下的素数。它还展示了其他能够在合理时间内生成十亿以下素数的基本 R 方法。我之所以这么说，是因为您似乎试图表明，这种生成 10 亿以下素数的练习不是一项可行的任务。如果我误解了，请纠正我。
@JosephWood 我没有读过它，因为它开始谈论包。但事实证明包裹是你的！这是非常令人印象深刻的！

【解决方案2】：

R 中的质数

OP 要求生成所有低于 10 亿的素数。到目前为止提供的所有答案要么无法做到这一点，要么执行需要很长时间，要么目前在 R 中不可用（参见@Charles 的answer）。包RcppAlgos（我是作者）能够仅使用一个线程在1 second 上生成请求的输出。它基于 Kim Walisch 的 Eratosthenes 分段筛。

RcppAlgos

library(RcppAlgos)
system.time(primeSieve(1e9))  ## using 1 thread
  user  system elapsed 
 1.099   0.077   1.176

使用多线程

在最近的版本中（即>= 2.3.0），我们可以利用多个线程来更快地生成。例如，现在我们可以在不到半秒的时间内生成高达 10 亿个素数！

system.time(primeSieve(10^9, nThreads = 8))
  user  system elapsed 
 2.046   0.048   0.375

R 中用于生成素数的可用包摘要

library(schoolmath)
library(primefactr)
library(sfsmisc)
library(primes)
library(numbers)
library(spuRs)
library(randtoolbox)
library(matlab)
## and 'sieve' from @John

在开始之前，我们注意到@Henrik 在schoolmath 中指出的问题仍然存在。观察：

## 1 is NOT a prime number
schoolmath::primes(start = 1, end = 20)
[1]  1  2  3  5  7 11 13 17 19   

## This should return 1, however it is saying that 52
##  "prime" numbers less than 10^4 are divisible by 7!!
sum(schoolmath::primes(start = 1, end = 10^4) %% 7L == 0)
[1] 52

关键是，此时不要使用schoolmath 来生成素数（作者无意冒犯……事实上，我已经向维护者提出了问题）。

让我们看看randtoolbox，因为它看起来非常高效。观察：

library(microbenchmark)
## the argument for get.primes is for how many prime numbers you need
## whereas most packages get all primes less than a certain number
microbenchmark(priRandtoolbox = get.primes(78498),
              priRcppAlgos = RcppAlgos::primeSieve(10^6), unit = "relative")
Unit: relative
          expr      min       lq     mean   median       uq       max neval
priRandtoolbox  1.00000  1.00000 1.000000 1.000000 1.000000 1.0000000   100
  priRcppAlgos 12.79832 12.55065 6.493295 7.355044 7.363331 0.3490306   100

仔细观察会发现它本质上是一个查找表（在源代码的文件randtoolbox.c 中找到）。

#include "primes.h"

void reconstruct_primes()
{
    int i;
    if (primeNumber[2] == 1)
        for (i = 2; i < 100000; i++)
            primeNumber[i] = primeNumber[i-1] + 2*primeNumber[i];
}

其中primes.h 是一个头文件，其中包含一个“素数差的一半”数组。因此，您将受到该数组中用于生成素数的元素数量（即前十万个素数）的限制。如果您只使用较小的素数（小于1,299,709（即第 100,000 个素数）），并且您正在处理需要nth 素数的项目，那么randtoolbox 是可行的方法。

下面，我们对其余包执行基准测试。

质数高达一百万

microbenchmark(priRcppAlgos = RcppAlgos::primeSieve(10^6),
               priNumbers = numbers::Primes(10^6),
               priSpuRs = spuRs::primesieve(c(), 2:10^6),
               priPrimes = primes::generate_primes(1, 10^6),
               priPrimefactr = primefactr::AllPrimesUpTo(10^6),
               priSfsmisc = sfsmisc::primes(10^6),
               priMatlab = matlab::primes(10^6),
               priJohnSieve = sieve(10^6),
               unit = "relative")
Unit: relative
          expr        min        lq      mean     median        uq       max neval
  priRcppAlgos   1.000000   1.00000   1.00000   1.000000   1.00000   1.00000   100
    priNumbers  21.550402  23.19917  26.67230  23.140031  24.56783  53.58169   100
      priSpuRs 232.682764 223.35847 233.65760 235.924538 236.09220 212.17140   100
     priPrimes  46.591868  43.64566  40.72524  39.106107  39.60530  36.47959   100
 priPrimefactr  39.609560  40.58511  42.64926  37.835497  38.89907  65.00466   100
    priSfsmisc   9.271614  10.68997  12.38100   9.761438  11.97680  38.12275   100
     priMatlab  21.756936  24.39900  27.08800  23.433433  24.85569  49.80532   100
  priJohnSieve  10.630835  11.46217  12.55619  10.792553  13.30264  38.99460   100

质数高达一千万

microbenchmark(priRcppAlgos = RcppAlgos::primeSieve(10^7),
               priNumbers = numbers::Primes(10^7),
               priSpuRs = spuRs::primesieve(c(), 2:10^7),
               priPrimes = primes::generate_primes(1, 10^7),
               priPrimefactr = primefactr::AllPrimesUpTo(10^7),
               priSfsmisc = sfsmisc::primes(10^7),
               priMatlab = matlab::primes(10^7),
               priJohnSieve = sieve(10^7),
               unit = "relative", times = 20)
Unit: relative
          expr       min        lq      mean    median        uq       max neval
  priRcppAlgos   1.00000   1.00000   1.00000   1.00000   1.00000   1.00000    20
    priNumbers  30.57896  28.91780  31.26486  30.47751  29.81762  40.43611    20
      priSpuRs 533.99400 497.20484 490.39989 494.89262 473.16314 470.87654    20
     priPrimes 125.04440 114.71349 112.30075 113.54464 107.92360 103.74659    20
 priPrimefactr  52.03477  50.32676  52.28153  51.72503  52.32880  59.55558    20
    priSfsmisc  16.89114  16.44673  17.48093  16.64139  18.07987  22.88660    20
     priMatlab  30.13476  28.30881  31.70260  30.73251  32.92625  41.21350    20
  priJohnSieve  18.25245  17.95183  19.08338  17.92877  18.35414  32.57675    20

质数高达一亿

对于接下来的两个基准测试，我们只考虑RcppAlgos、numbers、sfsmisc、matlab 和@John 的sieve 函数。

microbenchmark(priRcppAlgos = RcppAlgos::primeSieve(10^8),
               priNumbers = numbers::Primes(10^8),
               priSfsmisc = sfsmisc::primes(10^8),
               priMatlab = matlab::primes(10^8),
               priJohnSieve = sieve(10^8),
               unit = "relative", times = 20)
Unit: relative
         expr      min       lq     mean   median       uq      max neval
 priRcppAlgos  1.00000  1.00000  1.00000  1.00000  1.00000  1.00000    20
   priNumbers 35.64097 33.75777 32.83526 32.25151 31.74193 31.95457    20
   priSfsmisc 21.68673 20.47128 20.01984 19.65887 19.43016 19.51961    20
    priMatlab 35.34738 33.55789 32.67803 32.21343 31.56551 31.65399    20
 priJohnSieve 23.28720 22.19674 21.64982 21.27136 20.95323 21.31737    20

质数高达十亿

注意我们必须删除sieve函数中的条件if(n > 1e8) stop("n too large")。

## See top section
## system.time(primeSieve(10^9))
##  user  system elapsed 
## 1.099   0.077   1.176      ## RcppAlgos single-threaded

## gc()
system.time(matlab::primes(10^9))
   user  system elapsed 
 31.780  12.456  45.549        ## ~39x slower than RcppAlgos

## gc()
system.time(numbers::Primes(10^9))
   user  system elapsed 
 32.252   9.257  41.441        ## ~35x slower than RcppAlgos

## gc()
system.time(sieve(10^9))
  user  system elapsed 
26.266   3.906  30.201         ## ~26x slower than RcppAlgos

## gc()
system.time(sfsmisc::primes(10^9))
  user  system elapsed 
24.292   3.389  27.710         ## ~24x slower than RcppAlgos

从这些比较中，我们看到RcppAlgos 随着n 的增大而扩展得更好。

 _________________________________________________________
|            |   1e6   |   1e7    |   1e8     |    1e9    |
|            |---------|----------|-----------|-----------
| RcppAlgos  |   1.00  |   1.00   |    1.00   |    1.00   |
|   sfsmisc  |   9.76  |  16.64   |   19.66   |   23.56   |
| JohnSieve  |  10.79  |  17.93   |   21.27   |   25.68   |
|   numbers  |  23.14  |  30.48   |   32.25   |   34.86   |
|    matlab  |  23.43  |  30.73   |   32.21   |   38.73   |
 ---------------------------------------------------------

当我们使用多线程时，差异会更加显着：

microbenchmark(ser = primeSieve(1e6),
               par = primeSieve(1e6, nThreads = 8), unit = "relative")
Unit: relative
expr      min       lq     mean   median       uq      max neval
 ser 1.741342 1.492707 1.481546 1.512804 1.432601 1.275733   100
 par 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000   100

microbenchmark(ser = primeSieve(1e7),
               par = primeSieve(1e7, nThreads = 8), unit = "relative")
Unit: relative
 expr      min      lq     mean   median       uq      max neval
  ser 2.632054 2.50671 2.405262 2.418097 2.306008 2.246153   100
  par 1.000000 1.00000 1.000000 1.000000 1.000000 1.000000   100

microbenchmark(ser = primeSieve(1e8),
               par = primeSieve(1e8, nThreads = 8), unit = "relative", times = 20)
Unit: relative
 expr      min       lq     mean   median       uq      max neval
  ser 2.914836 2.850347 2.761313 2.709214 2.755683 2.438048    20
  par 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000    20

microbenchmark(ser = primeSieve(1e9),
               par = primeSieve(1e9, nThreads = 8), unit = "relative", times = 10)
Unit: relative
 expr      min       lq     mean   median       uq      max neval
  ser 3.081841 2.999521 2.980076 2.987556 2.961563 2.841023    10
  par 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000    10

然后将上表乘以系列结果的相应中值时间：

 _____________________________________________________________
|                |   1e6   |   1e7    |   1e8     |    1e9    |
|                |---------|----------|-----------|-----------
| RcppAlgos-Par  |   1.00  |   1.00   |    1.00   |    1.00   |
| RcppAlgos-Ser  |   1.51  |   2.42   |    2.71   |    2.99   |
|     sfsmisc    |  14.76  |  40.24   |   53.26   |   70.39   |
|   JohnSieve    |  16.32  |  43.36   |   57.62   |   76.72   |
|     numbers    |  35.01  |  73.70   |   87.37   |  104.15   |
|      matlab    |  35.44  |  74.31   |   87.26   |  115.71   |
 -------------------------------------------------------------

某个范围内的素数

microbenchmark(priRcppAlgos = RcppAlgos::primeSieve(10^9, 10^9 + 10^6),
               priNumbers = numbers::Primes(10^9, 10^9 + 10^6),
               priPrimes = primes::generate_primes(10^9, 10^9 + 10^6),
               unit = "relative", times = 20)
Unit: relative
         expr      min       lq    mean   median       uq      max neval
 priRcppAlgos   1.0000   1.0000   1.000   1.0000   1.0000   1.0000    20
   priNumbers 115.3000 112.1195 106.295 110.3327 104.9106  81.6943    20
    priPrimes 983.7902 948.4493 890.243 919.4345 867.5775 708.9603    20

在 6 秒内启动多达 100 亿次

##  primes less than 10 billion
system.time(tenBillion <- RcppAlgos::primeSieve(10^10, nThreads = 8))
  user  system elapsed 
26.077   2.063   5.602

length(tenBillion)
[1] 455052511

## Warning!!!... Large object created
tenBillionSize <- object.size(tenBillion)
print(tenBillionSize, units = "Gb")
3.4 Gb

非常大数范围内的质数：

在2.3.0 版本之前，我们只是对每个数量级的数字使用相同的算法。当大多数筛选素数在每个段中至少有一个倍数时，这对于较小的数字是可以的（通常，段大小受L1 Cache ~32KiB 的大小限制）。但是，当我们处理较大的数字时，筛选素数将包含许多数字，每个段的倍数少于一个。这种情况会产生很多开销，因为我们正在执行许多污染缓存的毫无价值的检查。因此，当数字非常大时，我们观察到素数的生成要慢得多。观察版本2.2.0（见Installing older version of R package）：

## Install version 2.2.0
## packageurl <- "http://cran.r-project.org/src/contrib/Archive/RcppAlgos/RcppAlgos_2.2.0.tar.gz"
## install.packages(packageurl, repos=NULL, type="source")

system.time(old <- RcppAlgos::primeSieve(1e15, 1e15 + 1e9))
 user  system elapsed 
7.932   0.134   8.067

现在使用最初由Tomás Oliveira 开发的缓存友好改进，我们看到了巨大的改进：

## Reinstall current version from CRAN
## install.packages("RcppAlgos"); library(RcppAlgos)
system.time(cacheFriendly <- primeSieve(1e15, 1e15 + 1e9))
 user  system elapsed 
2.258   0.166   2.424   ## Over 3x faster than older versions

system.time(primeSieve(1e15, 1e15 + 1e9, nThreads = 8))
 user  system elapsed 
4.852   0.780   0.911   ##  Over 8x faster using multiple threads

带走

有许多很棒的包可用于生成素数
如果您一般都在寻找速度，那么RcppAlgos::primeSieve 是不匹配的，尤其是对于更大的数字。
如果您正在使用小质数，请查看randtoolbox::get.primes。
如果您需要某个范围内的素数，则可以使用 numbers、primes 和 RcppAlgos 包。
良好编程实践的重要性怎么强调都不过分（例如矢量化、使用正确的数据类型等）。 @John 提供的纯基础 R 解决方案最恰当地证明了这一点。它简洁、清晰且非常高效。

【讨论】：

除了返回大量找到的素数之外，您的代码还不错；您不能在每个段页面的非复合材料上而不是巨大的结果数组上返回一个迭代器，并在不浪费时间的情况下减少内存使用吗？如果您只需要第 n 个素数或某个范围内的素数计数，则可以将这些作为直接与筛子（位？）一起使用的函数提供。
@GordonBGood，感谢您的评论。是的，你是对的。在当前的实现中有很多可以改进的地方。迭代器是一个，另一个处理更大的素数（例如大于 1e15）。我现在正在研究一个更好地处理后者的实现。在我说别的之前，巨大的功劳归功于 Kim Walisch 和 Tomas Oliviera e Silva。就此而言，也感谢您。您在 SoA 与 SoE 上的帖子非常鼓舞人心，我经常引用它们。当我更新以包含迭代器时，我一定会提到你。
谢谢，我们尽力了。虽然我们有一个沟通渠道，但我会告诉你，我正在研究一种新算法，当用可以编译为高效本机代码的语言实现时，它看起来会等于甚至超过 Kim Walisch 的素数筛例如 Rust、Nim、Haskell 和（当然）C/C++。希望您会在接下来的几周内看到它发布在某个地方。它比primesieve 更容易实现。
@GordonBGood，我迫不及待地想看到它。你将如何释放它？ GitHub？杂志？再次感谢您的所有贡献。
我最近在 Eratosthenes 的 Sieve 上看到了别人的 github 存储库，它似乎非常快（至少在高端桌面 CPU 上）并且因为它是用 C 编写的，尽管 you might like to look at it。我浏览了一下，认为我理解它是如何工作的，但是它非常复杂，并且在 2500 行 C 代码中有点难以阅读。我的开发版本可能要短得多，而且绝对不是 C - 我可能会用 Nim（本质上是 C 前端）编写它，它更简洁，但大多数人都可以轻松阅读。

【解决方案3】：

在 (a) 之前的每个数字 (i) 都会根据检查数字 (i-1) 生成的素数 (n) 列表进行检查

感谢您的建议：

prime = function(a,n){
    n=c(2)
    i=3
    while(i <=a){
      for(j in n[n<=sqrt(i)]){
        r=0
        if (i%%j == 0){
          r=1}
        if(r==1){break}
        
        
      }
      if(r!=1){n = c(n,i)}
      i=i+2
    }
    print(n)
  }

【讨论】：

【解决方案4】：

这种方法应该更快更简单。

allPrime <- function(n) {
  primes <- rep(TRUE, n)
  primes[1] <- FALSE
  for (i in 1:sqrt(n)) {
    if (primes[i]) primes[seq(i^2, n, i)] <- FALSE
  }
  which(primes)
}

n = 1e6 在我的电脑上为 0.12 秒

我在primefactr包中的函数AllPrimesUpTo中实现了这个。

【讨论】：

【解决方案5】：

for (i in 2:1000) {
a = (2:(i-1))
b = as.matrix(i%%a)
c = colSums(b != 0)
if (c == i-2)
 {
 print(i)
 }
 }

【讨论】：

这不是一个有用的答案，如果这段代码必须运行十亿个数字，更是如此。

【解决方案6】：

George Dontas 发布的那个筛子是一个很好的起点。这是一个更快的版本，1e6 素数的运行时间为 0.095 秒，而原始版本的运行时间为 30 秒。

sieve <- function(n)
{
   n <- as.integer(n)
   if(n > 1e8) stop("n too large")
   primes <- rep(TRUE, n)
   primes[1] <- FALSE
   last.prime <- 2L
   fsqr <- floor(sqrt(n))
   while (last.prime <= fsqr)
   {
      primes[seq.int(2L*last.prime, n, last.prime)] <- FALSE
      sel <- which(primes[(last.prime+1):(fsqr+1)])
      if(any(sel)){
        last.prime <- last.prime + min(sel)
      }else last.prime <- fsqr+1
   }
   which(primes)
}

下面是一些在 R 中以尽可能快的速度编码的替代算法。它们比筛子慢，但比提问者的原始帖子快得多。

这是一个使用 mod 但被矢量化的递归函数。它几乎立即返回 1e5，并在 2 秒内返回 1e6。

primes <- function(n){
    primesR <- function(p, i = 1){
        f <- p %% p[i] == 0 & p != p[i]
        if (any(f)){
            p <- primesR(p[!f], i+1)
        }
        p
    }
    primesR(2:n)
}

下一个不是递归的并且再次更快。下面的代码在我的机器上大约在 1.5 秒内启动到 1e6。

primest <- function(n){
    p <- 2:n
    i <- 1
    while (p[i] <= sqrt(n)) {
        p <-  p[p %% p[i] != 0 | p==p[i]]
        i <- i+1
    }
    p
}

顺便说一句，spuRs 包有许多素数查找功能，包括 E 的筛子。还没有检查它们的速度。

虽然我正在写一个很长的答案...如果一个值是素数，您将如何检查 R。

isPrime <- function(x){
    div <- 2:ceiling(sqrt(x))
    !any(x %% div == 0)
}

【讨论】：

您的 isPrime 函数不正确。如果传入 3，则返回 FALSE。
它也不能与 2 一起使用，现在仍然不行。我设计它只是为了快速检查我不知道的数字。我正在考虑用逻辑来修复它，但它很容易修复为 3，并且通过将 floor 更改为 ceiling（完成）只会稍微慢一点。
感谢您的精彩帖子！我想知道您的 isPrime 函数是否会通过将 2:ceiling(sqrt(x)) 替换为 seq_len(ceiling(sqrt(x)))[-1] 来改进
这将是一个有趣的代码子集来测试速度。 seq_len 应该比 : 快，但是否定选择可能会取消它。也许数量很大时会更好。
@John 在筛子模型中偶然我们可以计算范围之间的素数？

【解决方案7】：

上面发布的 isPrime() 函数可以使用 sieve()。只需要检查是否有任何素数

isPrime <- function(x) {
    div <- sieve(ceiling(sqrt(x)))
    (x > 1) & ((x == 2) | !any(x %% div == 0))
}

【讨论】：

sieve 来自哪个包？它返回什么，为什么会这样？
sieve() 函数来自“John”的上述帖子。它在 R 中实现了 Eratosthenes 筛法。您是否在问为什么我们只需要检查素数是否除而忽略非素数？这是基本的数论。
我不是这么问的。我以为您正在使用公开可用包中的功能，抱歉。但这意味着这可能不应该是一个答案，而是对约翰答案的评论。此外，我怀疑找到用于除法的素数非常昂贵，最好不要这样做。作为一个额外的问题，对于x 的哪些值，您在此计算中实际上需要ceiling 而不是floor？ floor 在数学上是正确的。
是的，应该是评论，但我不能评论约翰的回答，因为我的“声誉”Bernstein's primegen。

【解决方案8】：

这是 R 中 Sieve of Eratosthenes 算法的实现。

sieve <- function(n)
{
   n <- as.integer(n)
   if(n > 1e6) stop("n too large")
   primes <- rep(TRUE, n)
   primes[1] <- FALSE
   last.prime <- 2L
   for(i in last.prime:floor(sqrt(n)))
   {
      primes[seq.int(2L*last.prime, n, last.prime)] <- FALSE
      last.prime <- last.prime + min(which(primes[(last.prime+1):n]))
   }
   which(primes)
}

 sieve(1000000)

【讨论】：

这是算法的一个很好的实现，但是因为我们使用的是 R，所以它很慢......请参阅下面的回复。
嗨，我今天读到了一篇关于素数的有趣文章，关于Yitabg Zhang 和his work。事实上，Polish entry on Wikipedia for the Sieve 包含许多有用的实现。我建议你贡献你的。看到 R 中没有可用的实现，我很难过。
这不是我的。我在这里找到了它：rosettacode.org/wiki/Sieve_of_Eratosthenes#R

【解决方案9】：

我推荐primegen，这是 Dan Bernstein 对 Atkin-Bernstein 筛子的实现。它非常快，并且可以很好地扩展到其他问题。您需要将数据传递给程序才能使用它，但我想有办法做到这一点吗？

【讨论】：

看起来像一个 C 实现。通过正确的调整，您可以毫无问题地将 Primegen 合并到 R 中（参见 ?.C）
@Joris：我不知道你能做到这一点！感谢您的提醒。（我想把它称为外部。）
"primegen" 与Kim Walich's primesieve 相比并不是特别快，Kim Walich's primesieve 对于大于十亿的范围确实会显着减慢。它还有一个“C”接口，因此您可以以类似的方式从“R”调用它。
@GordonBGood 同意，不知道我当时在想什么。我现在要投票给 Joseph Wood 的答案。
@Charles，是的，我也刚刚投票给 Joseph Wood 的答案。

【解决方案10】：

你也可以作弊并使用schoolmath包中的primes()函数:D

【讨论】：

不！ Neil Gunther 指出，在学校数学中的 primes() 函数中有一个错误！在这里阅读：perfdynamics.blogspot.com/2010/06/… 也在这里：perfdynamics.blogspot.com/2010/06/… 和这里：xianblog.wordpress.com/2010/06/14/bug-in-schoolmath

【解决方案11】：

我所知道的生成所有素数的最佳方法（不涉及疯狂的数学）是使用Sieve of Eratosthenes。

实现起来非常简单，并且允许您在不使用除法或模数的情况下计算素数。唯一的缺点是它占用大量内存，但可以进行各种优化以提高内存（例如忽略所有偶数）。

【讨论】：