我正在尝试生成1亿以下的质数列表。我正在尝试这个，但是这种结构相当糟糕。有什么建议吗？
a <- 1:1000000000
d <- 0
b <- for (i in a) {for (j in 1:i) {if (i %% j !=0) {d <- c(d,i)}}}
生成一定范围内的素数列表
19
- Thomas
3
6你需要一次性获取所有质数吗？如果需要，似乎可以从这样的地方下载或获取它们：http://primes.utm.edu/lists/small/millions/。请注意，我的翻译只涉及转换原文内容成中文，不包括其他任何解释或信息。 - Chase
2在这里，您可以获取小于1000000000的质数列表 http://www.prime-numbers.org/ - Déjà vu
1如果你简单地将计算向量化，它们会比筛法运行得快得多，因为R的本质如此。（请查看我的回答） - John
11个回答
35
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)
}
- John
5
你的isPrime函数不正确。如果你传入3，它会返回FALSE。 - mchangun
1它也不能用于2，现在仍然不行。我只是为了快速检查我不知道的数字而设计它。我想用逻辑修复它，但对于3来说很容易修复，并且只需要稍微慢一点，通过将“floor”更改为“ceiling”（已完成）。 - John
感谢您的出色文章！我想知道，如果将2：ceiling(sqrt(x))替换为seq_len(ceiling(sqrt(x)))[-1]，您的isPrime函数是否会得到改进。 - De Novo
这将是一个有趣的代码子集，可以用于测试速度。
seq_len应该比:更快，但负选择可能会抵消它。也许当数字非常大时会更好。 - John@John 在筛法模型中，我们可以计算在某个范围内的质数吗？ - Nithish
23
这是 R 语言实现的埃拉托色尼筛法算法。
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)
- George Dontas
3
3这是一个很好的算法实现，但因为我们使用 R 语言，速度很慢...请见下面我的回应。 - John
嗨，今天我读到了一篇有趣的关于Yitabg Zhang和他在质数方面的工作的文章。事实上，维基百科上的波兰筛子条目包含了许多有用的实现。我建议你也做出自己的贡献。我很遗憾地发现R语言中没有可用的实现。 - Konrad
这不是我的。我在这里找到它：http://rosettacode.org/wiki/Sieve_of_Eratosthenes#R - gd047
9
R语言中的质数
原帖要求生成小于十亿的所有质数。目前提供的所有答案都不能胜任此任务，或者需要很长时间才能执行，或者在 R 中当前不可用（请参见 @Charles 的答案）。包 RcppAlgos（我是作者）可以在单线程模式下仅需 1秒 生成所请求的输出。它基于 Kim Walisch 的分段埃拉托色尼筛法（primesieve）。
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个质数）），并且正在进行需要获取第n个质数的项目，则可以选择使用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
从这些比较中，我们可以看到随着n的增大，
RcppAlgos的扩展性更好。 _________________________________________________________
| | 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 |
-------------------------------------------------------------
Primes Over a Range
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（请参见安装R软件包的旧版本）：
## 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解决方案最恰当地证明了这一点。该解决方案简洁、清晰且非常高效。
- Joseph Wood
10
你的代码很好，除了返回找到的大素数数组之外；你能否返回每个分段页面上非合数的迭代器，而不是巨大的结果数组，并减少内存使用而不失时间？如果你只需要第n个素数或在一个范围内的素数计数，这些可以作为直接使用筛法（位？）的函数提供。 - GordonBGood
@GordonBGood，感谢您的评论。是的，您说得对。当前实现中有很多可以改进的地方。迭代器是其中之一，另一个是处理更大的质数（例如大于1e15）。我正在开发一种更好地处理后者的实现。在我说其他任何话之前，要非常感谢Kim Walisch和Tomas Oliviera e Silva。同样，也要感谢您。您在SoA vs SoE上的帖子非常鼓舞人心，我经常参考它们。当我更新以包括迭代器时，我一定会提到您的。 - Joseph Wood
谢谢，我们会尽力而为。虽然我们有沟通渠道，但我想告诉你，我正在研究一种新的算法，看起来它将在编译成高效本地代码的语言（如Rust、Nim、Haskell和当然是C/C++）中实现时，能够达到或甚至超过Kim Walisch的primesieve。希望你能在接下来的几周内在某个地方看到它的发布。它比primesieve要简单一些实现。 - GordonBGood
@GordonBGood，我迫不及待地想看到它。你会如何发布它？Github？期刊？再次感谢您的所有贡献。 - Joseph Wood
1我最近看到了别人在Sieve of Eratosthenes上的github repo，它似乎非常快（至少在高端桌面CPU上），而且由于它是用C编写的，所以我认为你可能会喜欢看看它。我已经浏览过它并且认为我理解了它的工作原理，但它相当复杂，2500行的C代码有点难读。我的开发版本可能会更短，而且肯定不是C - 我可能会用Nim（本质上是C前端）来编写它，这样更简洁，但大多数人都可以轻松阅读它。 - GordonBGood
显示剩余5条评论
6
我所知道的生成所有质数的最佳方法（不涉及疯狂的数学）是使用埃拉托斯特尼筛法。它非常容易实现，可以在不使用除法或取模的情况下计算质数。唯一的缺点是它需要大量内存，但可以进行各种优化以改善内存（例如忽略所有偶数）。
- riwalk
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 函数中。
- F. Privé
3
我推荐使用primegen，它是Dan Bernstein实现的Atkin-Bernstein筛法。它非常快速，并且可以很好地适应其他问题。你需要将数据传递给程序才能使用它，但我想肯定有方法可以做到这一点。
- Charles
5
看起来像是C语言实现。通过适当的调整，你可以在R中轻松地集成primegen（参见?.C）。 - Joris Meys
@Joris：我不知道你可以这样做！谢谢你提醒我。（我一直想象它作为外部调用。） - Charles
与Kim Walich的primesieve相比，“primegen”并不特别快，对于大于十亿的范围，它会显著减慢速度。它还有一个“C”接口，因此您可以以类似的方式从“R”中调用它。 - GordonBGood
1@GordonBGood同意，当时我不确定自己在想什么。现在我要投票支持Joseph Wood的答案。 - Charles
1@Charles，是的，我也刚刚投了Joseph Wood的答案。 - GordonBGood
1
上面发布的isPrime()函数可以使用sieve()。只需要检查小于ceiling(sqrt(x))的质数是否能够整除x，并且没有余数。还需要处理1和2。
isPrime <- function(x) {
div <- sieve(ceiling(sqrt(x)))
(x > 1) & ((x == 2) | !any(x %% div == 0))
}
- John
4
sieve 是哪个包中的？它返回什么，为什么会起作用？ - Matthew Lundbergsieve() 函数来自于上面 "John" 的帖子。它在 R 中实现了埃拉托斯特尼筛法。你是否想知道为什么我们只需要检查质数是否能被整除而忽略非质数？这是基本的数论知识。 - John
1我并不是在问那个。我以为你是在使用一个公开可用的包中的函数，抱歉。但这意味着这可能不应该是一个答案，而应该是对John答案的评论。此外，我怀疑找到用于除法的质数是如此昂贵，以至于最好不要这样做。作为一个额外的问题，对于哪些值的x，您实际上需要使用ceiling而不是floor进行计算？从数学上讲，floor是正确的。 - Matthew Lundberg
1是的，这应该是一条评论，但由于我的“声望”<50，我不能对John的回答进行评论。John发布的sieve()函数不昂贵。如果执行时间是一个问题，那么应该使用Bernstein的primegen。 - John
1
你还可以作弊，使用 schoolmath 包中的 primes() 函数:D
- nico
1
10不行！Neil Gunther 指出了 schoolmath 中 primes() 函数的一个错误！在这里阅读：http://perfdynamics.blogspot.com/2010/06/primes-in-r-part-iii-schoolmath-is.html，另外还可以看看这里：http://perfdynamics.blogspot.com/2010/06/playing-with-primes-in-r-part-ii.html 和这里：http://xianblog.wordpress.com/2010/06/14/bug-in-schoolmath/。 - Henrik
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天。
- Antoni Parellada
4
你的试运行中获取了大约500万个质数，你的上限是多少？我猜测应该是1亿，但这对我来说不太清楚。其次，在你的
pi(100,000,000)语句中似乎缺少一个零（应该是pi(1,000,000,000)）。 - Joseph Wood@JosephWood 是的，就像原帖中所说的那样。感谢指出 pi 转录错误。 - Antoni Parellada
只是好奇，你看到我上面的答案了吗？它显示，不仅可以在一秒钟内生成十亿以下的质数，而且可以在大约5秒钟内获得小于一百亿的质数。它还展示了其他基于R的方法，能够在合理的时间内生成小于十亿的质数。我之所以这样说，是因为你似乎试图表明生成小于十亿的质数的练习是不可行的任务。如果我误解了，请纠正我。 - Joseph Wood
@JosephWood 我一开始没看懂，因为它一开始谈论的是包。但结果证明这个包是你写的！这非常令人印象深刻！ - Antoni Parellada
0
for (i in 2:1000) {
a = (2:(i-1))
b = as.matrix(i%%a)
c = colSums(b != 0)
if (c == i-2)
{
print(i)
}
}
- Hadi Atashi
1
1这不是一个有用的答案，尤其是如果这段代码需要运行数十亿个数字。 - vonbrand
网页内容由stack overflow 提供, 点击上面的可以查看英文原文，
原文链接
原文链接