Создайте список простых чисел до определенного числа

Question

Я пытаюсь создать список простых чисел ниже 1 миллиарда.Я пытаюсь это сделать, но такая структура довольно дерьмовая.Есть предложения?

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

Это сито, опубликованное Джорджем Донтасом, является хорошей отправной точкой. Вот гораздо более быстрая версия с временем выполнения для 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. Они медленнее, чем сито, но чертовски намного быстрее, чем оригинальный пост спрашивающих.

Вот рекурсивная функция, которая использует мод, но векторизована. Возвращается на 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)
}

Следующий не является рекурсивным и быстрее снова. Приведенный ниже код выполняет на моей машине до 1е6 примерно за 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
}

Кстати, пакет spuRs имеет ряд основных функций поиска, в том числе сито E. Не проверяли, какая у них скорость.

И хотя я пишу очень длинный ответ ... вот как вы можете проверить в R, является ли одно значение простым.

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

Answer 2

Это реализация алгоритма Сито Эратосфена в 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)

Answer 3

Лучший из известных мне способов получения всех простых чисел (без сумасшедшей математики) - это использовать Сито Эратосфена .

Это довольно просто реализовать и позволяет вычислять простые числа без использования деления или модуля.Единственным недостатком является то, что он требует большого объема памяти, но могут быть сделаны различные оптимизации для улучшения памяти (например, игнорируются все четные числа).

Answer 4

Простые числа в R

ОП попросил сгенерировать все простые числа меньше одного миллиарда.Все предоставленные до сих пор ответы либо не способны выполнить это, займут много времени для выполнения, либо в настоящее время недоступны в R (см. answer by @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,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^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

Теперь давайте удалим условие if(n > 1e8) stop("n too large") в функции sieve, чтобы проверить, насколько быстро он может генерировать простые числа до миллиарда, а также функцию primes из пакета sfsmisc, поскольку они были самыми быстрыми после RcppAlgos.

## gc()
system.time(sieve(10^9))
  user  system elapsed 
26.266   3.906  30.201

## gc()
system.time(sfsmisc::primes(10^9))
  user  system elapsed 
22.760   3.302  26.063

Исходя из этого, мы видим, что RcppAlgos масштабируется намного лучше, поскольку n становится больше (то есть 1.406 (см. Выше) примерно на 20-23x быстрее, тогда как былотолько около 10x для 10^6).

Простые числа до 10 миллиардов менее чем за 6 секунд

##  primes less than 10 billion
system.time(tenBillion <- RcppAlgos::primeSieve(10^10))
  user  system elapsed 
28.004   2.045   5.888

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.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

Take Away

1. Есть много отличных пакетов для генерации простых чисел
2. Если вы ищете скорость в целом, естьне соответствует RcppAlgos::primeSieve, особенно для больших чисел.
3. Если вы работаете с небольшими простыми числами, смотрите не дальше, чем randtoolbox::get.primes.
4. Если вам нужны простые числа в диапазоне,пакеты numbers, primes, & RcppAlgos - путь.
5. Важность хороших практик программирования невозможно переоценить (например, векторизация, использование правильных типов данных и т. Д.).Это наиболее удачно продемонстрировано решением с чистым основанием R, предоставленным @John.Это сжато, ясно и очень эффективно.

Answer 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.

Answer 6

Я рекомендую primegen , реализацию Дэна Бернштейна сита Аткина-Бернштейна.Это очень быстро и хорошо масштабируется на другие проблемы.Вам нужно будет передать данные программе, чтобы использовать ее, но я думаю, что есть способы сделать это?

Answer 7

Вы также можете обмануть и использовать функцию primes() в пакете schoolmath: D

Answer 8

Каждое число (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)
  }

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

Answer 10

Функция isPrime (), опубликованная выше, может использовать sieve (). Нужно только проверить, если какой-либо из простые числа <потолок (sqrt (x)) делят x без остатка. Нужно обрабатывать 1 и 2, а также. </p>

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