Я пытаюсь создать список простых чисел ниже 1 миллиарда. Я пытаюсь это сделать, но такая структура довольно дерьмовая. Любые предложения?
a <- 1:1000000000
d <- 0
b <- for (i in a) {for (j in 1:i) {if (i %% j !=0) {d <- c(d,i)}}}
Создайте список простых чисел до определенного числа
Ответ 1
Это реализация алгоритма Sieve of Eratoshenes в 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)
Ответ 2
Это сито, отправленное Джорджем Дантасом, является хорошей отправной точкой. Здесь гораздо более быстрая версия с временем выполнения для 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 почти мгновенно и 1e6 менее чем за 2 секунды.
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)
}
Следующий не рекурсивный и быстрый. В приведенном ниже коде на моей машине штрихи составляют до 1e6 примерно 1,5 с.
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
}
BTW, пакет spuRs имеет ряд основных функций поиска, включая сито E. Не проверял, чтобы посмотреть, какая скорость для них похожа.
И хотя я пишу очень длинный ответ... вот как вы проверите R, если одно значение является простым.
isPrime <- function(x){
div <- 2:ceiling(sqrt(x))
!any(x %% div == 0)
}
Ответ 3
Лучший способ, которым я знаю, генерировать все простые (без сумасшедшей математики), использовать Сито Эратосфена.
Это довольно просто реализовать и позволяет вычислять простые числа без использования деления или модуля. Единственным недостатком является интенсивность памяти, но для улучшения памяти могут быть сделаны различные оптимизации (без учета всех четных чисел).
Ответ 4
Простые числа в R
ОП попросил сгенерировать все простые числа ниже одного миллиарда. Все ответы, предоставленные до сих пор, либо не способны сделать это, могут занять много времени, либо в настоящее время недоступны в R (см. Ответ @Charles). Пакет RcppAlgos (я являюсь автором) способен генерировать запрошенный вывод всего за 1 second используя только один поток. Он основан на сегментированном сите Эратосфена Ким Валиш.
RcppAlgos
library(RcppAlgos)
system.time(primeSieve(10^9)) ## using 1 thread
user system elapsed
1.218 0.088 1.307
Использование нескольких потоков
А в последних версиях (т.е. >= 2.3.0) мы можем использовать несколько потоков для еще более быстрой генерации. Например, теперь мы можем генерировать простые числа до 1 миллиарда менее чем за полсекунды!
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 1,299,709 (т. 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 от @John.
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
Простые до 10 миллиардов менее чем за 6 секунд
## 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). Однако, когда мы имеем дело с большими числами, простые числа просеивания будут содержать много чисел, которые будут иметь менее одного кратного на сегмент. Эта ситуация создает много накладных расходов, так как мы выполняем много бесполезных проверок, которые загрязняют кеш. Таким образом, мы наблюдаем гораздо более медленную генерацию простых чисел, когда числа очень велики. Обратите внимание на версию 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
И теперь, используя улучшение, дружественное к кешу, первоначально разработанное Томасом Оливейрой, мы видим радикальные улучшения:
## 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
Увезти
- Есть много отличных пакетов, доступных для генерации простых чисел
- Если вы ищете скорость в целом,
RcppAlgos::primeSieve, особенно для больших чисел. - Если вы работаете с небольшими простыми числами, смотрите не дальше, чем
randtoolbox::get.primes. - Если вам нужны простые числа из диапазона,
numbersпакетов,primesиRcppAlgos- это то, что вам нужно. - Важность хороших практик программирования невозможно переоценить (например, векторизация, использование правильных типов данных и т.д.). Это наиболее удачно продемонстрировано решением с чистым основанием R, предоставленным @John. Это сжато, ясно, и очень эффективно.
Ответ 5
Этот метод должен быть более быстрым и простым.
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 секунды на моем компьютере для n = 1e6
Я реализовал это в функции AllPrimesUpTo в пакете primefactr.
Ответ 6
Я рекомендую primegen, внедрение Дан Бернштейна в сите Аткина-Бернштейна. Он очень быстрый и хорошо масштабируется для других проблем. Вам нужно будет передать данные программе, чтобы использовать ее, но я полагаю, что есть способы сделать это?
Ответ 7
Вы также можете обмануть и использовать функцию primes() в пакете schoolmath: D
Ответ 8
Функция isPrime(), опубликованная выше, может использовать сито(). Нужно только проверить, не простые числа < потолок (sqrt (x)) делит x без остатка. Необходимо также обрабатывать 1 и 2.
isPrime <- function(x) {
div <- sieve(ceiling(sqrt(x)))
(x > 1) & ((x == 2) | !any(x %% div == 0))
}
Ответ 9
for (i in 2:1000) {
a = (2:(i-1))
b = as.matrix(i%%a)
c = colSums(b != 0)
if (c == i-2)
{
print(i)
}
}
Ответ 10
Каждое число (i) перед (a) проверяется на список простых чисел (n), сгенерированных путем проверки числа (i-1)
Спасибо за предложения:
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)
}