生成一个特定数量的素数列表

Question

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

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

Answer 1

George Dontas发布的筛子是一个很好的起点.这是一个更快的版本,运行时间为1个6个素数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,2s内返回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)
}

Answer 2

这是R中的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)

Answer 3

我知道生成所有素数(不需要疯狂数学)的最佳方法是使用Eratosthenes的Sieve.

它非常简单,允许您在不使用除法或模数的情况下计算素数.唯一的缺点是内存密集,但可以进行各种优化以改善内存(例如忽略所有偶数).

Answer 4

R中的质数

OP要求生成所有十亿以下的质数。到目前为止提供的所有答案都不能执行此操作，将需要很长的时间执行，或者目前在R中不可用（请参阅@Charles 的答案）。该包RcppAlgos（我是作者）能够仅1 second使用一个线程就可以生成请求的输出。它基于Kim Walisch制作的Eratosthenes筛分筛。

Rcpp算法

library(RcppAlgos)
system.time(primeSieve(10^9))  ## using 1 thread
  user  system elapsed 
 1.218   0.088   1.307

使用多线程

在最新版本（即>= 2.3.0）中，我们可以利用多个线程来加快生成速度。例如，现在我们可以在半秒内生成多达10亿个质数！

system.time(primeSieve(10^9, nThreads = 8))
  user  system elapsed 
 2.239   0.046   0.416

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.0000  1.00000 1.000000 1.00000 1.000000 1.000000   100
  priRcppAlgos 14.0758 14.20469 8.555965 6.99534 7.114415 2.809296   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.000000   100
    priNumbers  19.092499  22.29017  25.79069  22.527344  23.53524 16.439443   100
      priSpuRs 210.980827 204.75970 203.74259 218.533689 218.12819 64.208745   100
     priPrimes  43.572518  40.61982  36.36935  37.234043  37.55404 10.399216   100
 priPrimefactr  35.850982  37.38720  39.47520  35.848167  37.62628 19.540713   100
    priSfsmisc   9.462374  10.54760  10.55800   9.921876  12.25639  3.887074   100
     priMatlab  19.698811  22.70576  25.39655  22.851422  23.63050 15.265014   100
  priJohnSieve  10.149523  10.68950  11.42043  10.437246  12.72949 11.595701   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  28.39102  27.63922  27.96319  27.34067  25.44119  37.72224    20
      priSpuRs 469.06554 469.09160 445.61612 457.34482 419.91417 398.29053    20
     priPrimes 117.11150 111.35547 107.61258 109.10053 102.32481  97.34148    20
 priPrimefactr  46.13612  47.24150  47.65271  47.34762  46.58394  50.10061    20
    priSfsmisc  17.37116  16.99990  17.64440  16.77242  17.10034  25.25716    20
     priMatlab  27.24177  27.17770  28.79239  27.37511  26.70660  36.79823    20
  priJohnSieve  16.83924  17.43330  18.63179  17.83366  17.24865  28.89491    20

灌注至一亿

在接下来的两个标杆，我们只考虑RcppAlgos，numbers，sfsmisc以及sieve通过@约翰功能。

microbenchmark(priRcppAlgos = RcppAlgos::primeSieve(10^8),
               priNumbers = numbers::Primes(10^8),
               priSfsmisc = sfsmisc::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 31.89653 30.93312 30.73546 30.70144 30.20808 28.79867    20
   priSfsmisc 21.13420 20.14822 19.84391 19.77317 19.40612 18.05891    20
 priJohnSieve 21.39554 20.24689 20.34909 20.24419 20.09711 19.16832    20

灌注至十亿

NB我们必须删除的条件if(n > 1e8) stop("n too large")的sieve功能。

## See top section
## system.time(primeSieve(10^9)). 
##  user  system elapsed 
## 1.218   0.088   1.307

## gc()
system.time(numbers::Primes(10^9))
   user  system elapsed 
 32.375  12.129  45.651        ## ~35x slower than RcppAlgos

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

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

从这些比较中，我们看到RcppAlgos随着n变大，缩放比例会好得多。

 _________________________________________________________
|            |   1e6   |   1e7    |   1e8     |    1e9    |
|            |---------|----------|-----------|-----------
| RcppAlgos  |   1.00  |   1.00   |    1.00   |    1.00   |
|   sfsmisc  |   9.92  |  16.77   |   19.77   |   21.20   |
| JohnSieve  |  10.44  |  17.83   |   20.24   |   23.11   |
|   numbers  |  22.53  |  27.34   |   30.70   |   34.93   |
 ---------------------------------------------------------

在范围内灌注

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 
27.319   1.971   5.822

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）。但是，当我们处理较大的数字时，筛分质数将包含许多数字，每个数字段的乘数少于1。这种情况会产生大量开销，因为我们正在执行许多无用的检查，从而污染了缓存。因此，当数量很大时，我们观察到质数的生成要慢得多。观察版本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ásOliveira最初开发的对缓存友好的改进，我们看到了巨大的改进：

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

system.time(primeSieve(1e15, 1e15 + 1e9, nThreads = 8))
 user  system elapsed 
5.037   0.806   0.981   ##  Over 8x faster using multiple threads

带走

1. 有很多很棒的软件包可以生成素数
2. 如果您一般都在寻找速度，则无法找到RcppAlgos::primeSieve，尤其是对于较大的数字。
3. 如果您使用的是小底数，不要比randtoolbox::get.primes。
4. 如果需要一定范围内的填料numbers，则可以使用primes，，＆包RcppAlgos。
5. 良好编程习惯的重要性不能过分强调（例如矢量化，使用正确的数据类型等）。@John提供的纯碱R解决方案最恰当地证明了这一点。它简洁，清晰且非常有效。