Genera una lista de números primos hasta cierto número

Intento generar una lista de primos por debajo de mil millones. Estoy intentando esto, pero este tipo de estructura es bastante mierda. ¿Alguna sugerencia?

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

Esta es una implementación del algoritmo Sieve of Eratosthenes en 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) 

Ese tamiz publicado por George Dontas es un buen punto de partida. Aquí hay una versión mucho más rápida con tiempos de ejecución para 1e6 primos de 0.095 en comparación con 30s para la versión original.

 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) } 

Aquí hay algunos algoritmos alternativos codificados lo más rápido posible en R. Son más lentos que el tamiz, pero muchísimo más rápidos que la publicación original del cuestionario.

Aquí hay una función recursiva que usa mod pero está vectorizada. Regresa por 1e5 casi instantáneamente y 1e6 en menos de 2s.

 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) } 

El siguiente no es recursivo y más rápido de nuevo. El código de abajo hace primos hasta 1e6 en aproximadamente 1.5s en mi máquina.

 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 } 

Por cierto, el paquete spuRs tiene una serie de funciones principales de búsqueda que incluyen un tamiz de E. No se han comprobado para ver cuál es la velocidad para ellos.

Y mientras escribo una respuesta muy larga ... así es como registrarías R si un valor es primo.

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

La mejor manera que conozco de generar todos los números primos (sin entrar en matemáticas locas) es usar el Tamiz de Eratóstenes .

Es bastante sencillo de implementar y te permite calcular números primos sin usar división o módulo. El único inconveniente es que consume mucha memoria, pero se pueden realizar varias optimizaciones para mejorar la memoria (por ejemplo, ignorando todos los números pares).

Este método debería ser más rápido y simple.

 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) } 

0,12 segundos en mi computadora para n = 1e6

Implementé esto en la función AllPrimesUpTo en el paquete primefactr.

Recomiendo la implementación de Primegen , Dan Bernstein del tamiz Atkin-Bernstein. Es muy rápido y se adaptará bien a otros problemas. Tendrá que pasar los datos al progtwig para usarlo, pero me imagino que hay formas de hacerlo.

También puede hacer trampa y usar la función primes() en el paquete schoolmath : D

El OP pidió generar todos los números primos por debajo de mil millones. Todas las respuestas proporcionadas hasta el momento no son capaces de hacerlo, tomarán mucho tiempo para ejecutarse, o actualmente no están disponibles en R (vea la respuesta de @Charles). El paquete RcppAlgos (soy el autor) es capaz de generar la salida solicitada en poco más de 1 second . Se basa en el tamiz segmentado de Eratóstenes de Kim Walisch .

 library(RcppAlgos) system.time(primeSieve(10^9)) user system elapsed 1.300 0.105 1.406 

Además, a continuación se muestra un resumen de los paquetes (y la función de sieve anterior proporcionada por @John) que pueden generar números primos.

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

Antes de comenzar, notamos que los problemas señalados por @Henrik en la schoolmath todavía existen. Observar:

 ## 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.... Huuuhhh???? sum(schoolmath::primes(start = 1, end = 10^4) %% 7L == 0) [1] 52 

El punto es que no use schoolmath para generar números primos en este punto (sin ofender al autor … De hecho, he archivado un problema con el mantenedor).

Echemos un vistazo a randtoolbox ya que parece ser increíblemente eficiente. Observar:

 ## 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.00000 1.00000 1.000000 1.000000 100 priRcppAlgos 19.37208 18.30095 9.897374 7.639231 7.249031 5.449076 100 

Una mirada más cercana revela que es esencialmente una tabla de búsqueda (que se encuentra en el archivo randtoolbox.c del código fuente).

 #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]; } 

Donde primes.h es un archivo de encabezado que contiene una matriz de "mitades de diferencias entre números primos" . Por lo tanto, estará limitado por la cantidad de elementos en esa matriz para generar primos (es decir, los primeros cien mil primos). Si solo está trabajando con números primos más pequeños (menos de 1,299,709 (es decir, la primo 1,299,709 ) y está trabajando en un proyecto que requiere la nth prima, randtoolbox es el camino a seguir.

Ahora, veamos cómo se comparan el rest de los paquetes en generar números primos de menos de un millón:

 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 21.599399 21.27664 20.17172 20.348373 20.82308 10.992780 100 priSpuRs 223.240897 231.69938 198.16347 215.520998 202.37479 59.608137 100 priPrimes 41.689864 38.43470 33.74977 35.329320 33.21983 17.323242 100 priPrimefactr 39.452661 39.77808 38.64081 37.887232 36.71392 14.549090 100 priSfsmisc 9.667778 11.16356 11.58725 10.862231 11.61987 8.990612 100 priMatlab 21.055741 21.45761 21.46455 21.053058 20.86727 15.029687 100 priJohnSieve 9.316246 10.22913 10.52325 9.825776 10.30900 8.167797 100 

Probemos la velocidad de generación de primos en un rango:

 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.0000 1.0000 1.0000 1.00000 20 priNumbers 137.5877 134.2035 125.1085 133.1225 121.9784 94.93077 20 priPrimes 911.2896 877.6692 806.9694 861.4054 783.9666 568.05759 20 20 

Ahora, eliminemos la condición if(n > 1e8) stop("n too large") en la función sieve para probar qué tan rápido puede generar los números primos por debajo de mil millones, así como la función de primes del paquete sfsmisc , ya que fueron más rápido después de RcppAlgos .

 system.time(sieve(10^9)) user system elapsed 26.703 4.582 31.328 system.time(sfsmisc::primes(10^9)) user system elapsed 24.772 4.521 29.436 

A partir de esto, vemos que RcppAlgos escala mucho mejor a medida que n se hace más grande (es decir, 1.406 (ver arriba) es aproximadamente 20-22x más rápido mientras que solo fue alrededor de 10x para 10^6 ).

Y solo para una buena medida:

 ## primes less than 10 billion system.time(tenBillion <- RcppAlgos::primeSieve(10^10)) user system elapsed 16.510 1.784 18.296 length(tenBillion) [1] 455052511 ## Warning!!!... Large object created tenBillionSize <- object.size(tenBillion) print(tenBillionSize, units = "Gb") 3.4 Gb 

Para llevar

  1. Hay muchos paquetes excelentes disponibles para generar números primos
  2. Si está buscando velocidad en general, no hay coincidencia con RcppAlgos::primeSieve .
  3. Si está trabajando con primos pequeños, no busque más que randtoolbox::get.primes .
  4. Si necesita números primos en un rango, los numbers paquetes, numbers primes y RcppAlgos son el camino a seguir.
  5. La importancia de las buenas prácticas de progtwigción no se puede exagerar (por ejemplo, vectorización, uso de tipos de datos correctos, etc.). Esto se demuestra muy bien con la solución de base pura R proporcionada por @John. Es conciso, claro y muy eficiente.

La función isPrime () publicada anteriormente podría usar tamve (). Solo se necesita verificar si alguno de los números primos

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

Cada número (i) antes (a) se compara con la lista de números primos (n) generados al verificar el número (i-1)

Gracias por las sugerencias:

 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) }