Generate a list of primes up to a certain number

Question

I'm trying to generate a list of primes below 1 billion. I'm trying this, but this kind of structure is pretty shitty. Any suggestions?

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

Solution 1

This is an implementation of the Sieve of Eratosthenes algorithm in 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)

Solution 2

That sieve posted by George Dontas is a good starting point. Here's a much faster version with running times for 1e6 primes of 0.095s as opposed to 30s for the original version.

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

Here are some alternate algorithms below coded about as fast as possible in R. They are slower than the sieve but a heck of a lot faster than the questioners original post.

Here's a recursive function that uses mod but is vectorized. It returns for 1e5 almost instantaneously and 1e6 in under 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)
}

The next one isn't recursive and faster again. The code below does primes up to 1e6 in about 1.5s on my machine.

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, the spuRs package has a number of prime finding functions including a sieve of E. Haven't checked to see what the speed is like for them.

And while I'm writing a very long answer... here's how you'd check in R if one value is prime.

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

Solution 3

Prime Numbers in R

The OP asked to generate all prime numbers below one billion. All of the answers provided thus far are either not capable of doing this, will take a long a time to execute, or currently not available in R (see the answer by @Charles). The package RcppAlgos (I am the author) is capable of generating the requested output in just over 1 second using only one thread. It is based off of the segmented sieve of Eratosthenes by Kim Walisch.

RcppAlgos

library(RcppAlgos)
system.time(primeSieve(1e9))  ## using 1 thread
  user  system elapsed 
 1.099   0.077   1.176 

Using Multiple Threads

And in recent versions (i.e. >= 2.3.0), we can utilize multiple threads for even faster generation. For example, now we can generate the primes up to 1 billion in under half a second!

system.time(primeSieve(10^9, nThreads = 8))
  user  system elapsed 
 2.046   0.048   0.375

Summary of Available Packages in R for Generating Primes

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

Before we begin, we note that the problems pointed out by @Henrik in schoolmath still exists. Observe:

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

The point is, don't use schoolmath for generating primes at this point (no offense to the author... In fact, I have filed an issue with the maintainer).

Let's look at randtoolbox as it appears to be incredibly efficient. Observe:

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

A closer look reveals that it is essentially a lookup table (found in the file randtoolbox.c from the source code).

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

Where primes.h is a header file that contains an array of "halves of differences between prime numbers". Thus, you will be limited by the number of elements in that array for generating primes (i.e. the first one hundred thousand primes). If you are only working with smaller primes (less than 1,299,709 (i.e. the 100,000th prime)) and you are working on a project that requires the nth prime, randtoolbox is the way to go.

Below, we perform benchmarks on the rest of the packages.

Primes up to One Million

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

Primes up to Ten Million

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

Primes up to One Hundred Million

For the next two benchmarks, we only consider RcppAlgos, numbers, sfsmisc, matlab, and the sieve function by @John.

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

Primes up to One Billion

N.B. We must remove the condition if(n > 1e8) stop("n too large") in the sieve function.

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

From these comparison, we see that RcppAlgos scales much better as n gets larger.

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

The difference is even more dramatic when we utilize multiple threads:

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

And multiplying the table above by the respective median times for the serial results:

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

Primes up to 10 billion in Under 6 Seconds

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

Primes Over a Range of Very Large Numbers:

Prior to version 2.3.0, we were simply using the same algorithm for numbers of every magnitude. This is okay for smaller numbers when most of the sieving primes have at least one multiple in each segment (Generally, the segment size is limited by the size of L1 Cache ~32KiB). However, when we are dealing with larger numbers, the sieving primes will contain many numbers that will have fewer than one multiple per segment. This situation creates a lot of overhead, as we are performing many worthless checks that pollutes the cache. Thus, we observe much slower generation of primes when the numbers are very large. Observe for version 2.2.0 (See Installing older version of R package):

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

And now using the cache friendly improvement originally developed by Tomás Oliveira, we see drastic improvements:

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

Take Away

1. There are many great packages available for generating primes
2. If you are looking for speed in general, there is no match to RcppAlgos::primeSieve, especially for larger numbers.
3. If you are working with small primes, look no further than randtoolbox::get.primes.
4. If you need primes in a range, the packages numbers, primes, & RcppAlgos are the way to go.
5. The importance of good programming practices cannot be overemphasized (e.g. vectorization, using correct data types, etc.). This is most aptly demonstrated by the pure base R solution provided by @John. It is concise, clear, and very efficient.

Solution 4

Best way that I know of to generate all primes (without getting into crazy math) is to use the Sieve of Eratosthenes.

It is pretty straightforward to implement and allows you calculate primes without using division or modulus. The only downside is that it is memory intensive, but various optimizations can be made to improve memory (ignoring all even numbers for instance).

Solution 5

This method should be Faster and simpler.

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 second on my computer for n = 1e6

I implemented this in function AllPrimesUpTo in package primefactr.

View more solutions

Comments

[Henrik]

NO! There is a bug in the primes() function in schoolmath pointed at by Neil Gunther! Read here: perfdynamics.blogspot.com/2010/06/… see also here: perfdynamics.blogspot.com/2010/06/… and here: xianblog.wordpress.com/2010/06/14/bug-in-schoolmath

[John]

This is a good implementation of the algorithm but because we're using R it's SLOW... see my response below.

[Joris Meys]

Looks like a C-implementation. With the right tweaking you can incorporate primegen in R without problem (see ?.C)

[Charles]

@Joris: I didn't know you could do that! Thanks for the heads-up. (I was imagining calling it as an external.)

[mchangun]

Your isPrime function is not correct. If you pass 3 in, it returns FALSE.

[John]

It didn't work with 2 either, and still doesn't. I designed it just for numbers I don't know as a quick check. I was thinking of fixing it with logicals but it's easily fixed for 3, and only made a tiny bit slower, by changing floor to ceiling (done).

[Matthew Lundberg]

What package is sieve from? What does it return, and why does this work?

[John]

The sieve() function is from the post above by "John." It implements the Sieve of Eratosthenes in R. Are you asking why we only have to check if primes divide and ignore non-primes? That's basic number theory.

[Matthew Lundberg]

I was not asking that. I thought you were using a function in a publicly-available package, sorry. But that means that this probably shouldn't be an answer, but a comment on John's answer. In addition, I suspect that finding the primes to use for division is so expensive, that you are better off not doing so. As a bonus question, for what values of x do you actually need ceiling rather than floor in this computation? floor is correct mathematically.

[John]

Yes, it should be a comment, but I am not allowed to comment on John's answer, since my "reputation" is < 50. The sieve() function John posted is not expensive. If execution time is a concern, than one should use Bernstein's primegen.

[Konrad]

Hi, I was reading about that today as I came across an interesting article about Yitabg Zhang and his work on the prime numbers. As a matter fact, Polish entry on Wikipedia for the Sieve contains numerous useful implementations. I suggest that you contribute yours. I was sadden to see that there were no implementation in R available.

[George Dontas]

It's not mine. I have found it here: rosettacode.org/wiki/Sieve_of_Eratosthenes#R

[vonbrand]

That isn't an useful answer, even more so if this code has to run for billions of numbers.

[De Novo]

Thanks for an excellent post! I wonder if your isPrime function would be improved by replacing 2:ceiling(sqrt(x)) with seq_len(ceiling(sqrt(x)))[-1]

[John]

That would be an interesting subset of the code to test for speed. seq_len should be faster than : however the negative selection might cancel that out. Maybe it will be better when the number is very large.

[GordonBGood]

"primegen" isn't particularly fast as compared to Kim Walich's primesieve, which does slow down dramatically for ranges larger than a billion. It also has a "C" interface so you could call it from "R" in a similar way.

[Charles]

@GordonBGood Agreed, not sure what I was thinking at the time. I’m going to vote for Joseph Wood’s answer now.

[GordonBGood]

@Charles, yes, I just voted for Joseph Wood's answer, too.

[GordonBGood]

your code is good other than returning the huge array of found primes; couldn't you return an iterator over the non-composites per segment page rather than the huge results array and reduce the memory use without losing time? And if all you need are the nth prime or the count of primes over a range, these can be provided as functions that work with the sieve (bits?) directly.

[Joseph Wood]

@GordonBGood, thanks for your comment. Yes, you are correct. There is a lot that can be improved in the current implementation. The iterator is one, the other is dealing with larger primes (e.g. greater than 1e15). I’m working on an implementation right now that better deals with the latter. Before I say anything else, huge credit goes to Kim Walisch and Tomas Oliviera e Silva. For that matter, credit to you as well. Your posts on SoA vs SoE were very inspiring and I reference them frequently. When I update to include iterators I will be sure to mention you.

[GordonBGood]

Thank you, we do our best. While we have a channel of communication, I'll let you know that I'm working on a new algorithm that looks like it will equal or even beat Kim Walisch's primesieve by a little bit when inplemented in languages that can compile to efficient native code such as Rust, Nim, Haskell, and (of course) C/C++. Hopefully you'll see it posted somewhere within the next few weeks. It is somewhat simpler to implement than primesieve.

[Joseph Wood]

@GordonBGood, I cannot wait to see it. How will you release it? Github? Journal? Again, thanks for all your contributions.

[GordonBGood]

it may be big enough to justify a GitHub repo at over 500 lines of Haskell code (I hope you read Haskell). I'm writing in Haskell to show that Haskell can be as fast as C/C++, but it's the algorithm that gives the speed; it can be written in JVM languages or CLI languages like F# and C# although then we can't use some optimizations such as extreme loop unrolling and the code will be slower with JIT compilation and enforced array bounds checks. Runs almost as fast as the first class in Julia though, which is pretty amazing.

[GordonBGood]

I recently saw someone else's github repo on the Sieve of Eratosthenes, which seems to be very fast (at least on high end desktop CPU's) and since it is written in C I though you might like to look at it. I've skimmed through it and think I understand how it works, but it is quite complex and a bit hard to read at 2500 lines of C code. My developing version is likely to be much shorter and definitely not C - I'm may write it in Nim (essentially a C front end), which is more concise, yet most can read it easily.

[Joseph Wood]

For your trial run to get the ~5 million primes, what was your upper bound? I'm guessing it was one hundred million but that isn't clear to me. Secondly, it looks like you are missing a zero in your pi(100,000,000) statement (It should be pi(1,000,000,000)).

[Antoni Parellada]

@JosephWood Yes, as in the OP. Thanks for pointing out the pi transcription error.

[Joseph Wood]

Just curious, did you see my answer above? It shows that, not only can you generate primes under a billion in under a second, but you can get the primes under ten billion in around 5 seconds. It also shows other base R methods that are able to generate primes under a billion in a reasonable amount of time. I say all this because it appears you are trying to show that this exercise of generating primes under a billion is not a feasible tasks. Please correct me if I misunderstand.

[Antoni Parellada]

@JosephWood I hadn't read it because it started off talking about packages. But it turns out the package is yours! That is very impressive!

[Antoni Parellada]

Why can't I install your package on Google Colab?

[Joseph Wood]

@AntoniParellada I've never used Google Colab until just now and looks like there is an issue gmp which RcppAlgos relies on. Are you seeing Warning message in install.packages("gmp"): “installation of package ‘gmp’ had non-zero exit status”? Just try installing gmp and you will see the error.

[Antoni Parellada]

Yes, that is the error. I try installing gmp separately without success. Google Colab is awesome if your IT guys think that R packages can be a threat to the organization, because you work on a VM on the Google servers.

[Joseph Wood]

I honestly have no idea how to remedy this. My guess is that the C library underneath gmp isn't installed on that machine. That's just a guess as I have no experience with Google Colab. This could be worthy of a proper SO question.

[Nithish]

@John in sieve model by chance we can do calculate prime between ranges ?