primes function

Find all Primes Less Than n

Find all prime numbers aka primes less than $n$.

Uses an obvious sieve method (and some care), working with logical and and integers to be quite fast.

primes(n, pSeq = NULL)

Arguments

• n: a (typically positive integer) number.
• pSeq: optionally a vector of primes (2,3,5,...) as if from a primes() call; must be correct. The goal is a speedup, but currently we have not found one single case, where using a non-NULL pSeq is faster.

Details

As the function only uses max(n), n can also be a vector of numbers.

The famous prime number theorem states that $\pi(n)$, the number of primes below $n$ is asymptotically c("$n /\n$", "$\\log(n)$") in the sense that $lim[n -> Inf] \pi(n) * log(n) / n ~ 1$.

Equivalently, the inverse of $pi()$, the $n$-th prime number $p_n$ is around $n \log n$; recent results (Pierre Dusart, 1999), prove that

Returns

numeric vector of all prime numbers $<= n$.

Author(s)

Bill Venables (<= 2001); Martin Maechler gained another 40% speed, carefully working with logicals and integers.

See Also

factorize. For large $n$, use the list("gmp") package and its isprime and nextprime

functions.

Examples

(p1 <- primes(100))
 system.time(p1k <- primes(1000)) # still lightning fast
 stopifnot(length(p1k) == 168)

 system.time(p.e7 <- primes(1e7)) # still only 0.3 sec (2015 (i7))
 stopifnot(length(p.e7) == 664579)
 ## The famous  pi(n) :=  number of primes <= n:
 pi.n <- approxfun(p.e7, seq_along(p.e7), method = "constant")
 pi.n(c(10, 100, 1000)) # 4 25 168
 plot(pi.n, 2, 1e7, n = 1024, log="xy", axes = FALSE,
      xlab = "n", ylab = quote(pi(n)),
      main = quote("The prime number function " ~ pi(n)))
 eaxis(1); eaxis(2)

## Exploring  p(n) := the n-th prime number  ~=~ n * pnn(n), where
## pnn(n) := log n + log log n
pnn <- function(n) { L <- log(n); L + log(L) }
n <- 6:(N <- length(PR <- primes(1e5)))
m.pn <- cbind(l.pn = ceiling(n*(pnn(n)-1)), pn = PR[n], u.pn = floor(n*pnn(n)))
matplot(n, m.pn, type="l", ylab = quote(p[n]), main = quote(p[n] ~~
        "with lower/upper bounds" ~ n*(log(n) + log(log(n)) -(1~"or"~0))))
## (difference to the lower approximation) / n   --> ~ 0.0426  (?) :
plot(n, PR[n]/n - (pnn(n)-1), type = 'l', cex = 1/8, log="x", xaxt="n")
eaxis(1); abline(h=0, col=adjustcolor(1, 0.5))