sfsmisc (version 1.1-5)

primes: Find all Primes Less Than n

Description

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.

Usage

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.

Value

numeric vector of all prime numbers \(\le n\).

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 \(n / \log(n)\) in the sense that \(\lim_{n \to \infty}{\pi(n) \cdot \log(n) / n \sim 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 $$\log n + \log\log n - 1 < \frac{p_n}{n} < \log n + \log \log n \quad\mathrm{for } n \ge 6.$$

See Also

factorize. For large \(n\), use the gmp package and its isprime and nextprime functions.

Examples

# NOT RUN {
 (p1 <- primes(100))
 system.time(p1k <- primes(1000)) # still lightning fast
 stopifnot(length(p1k) == 168)
# }
# NOT RUN {
 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)
# }
# NOT RUN {
## Exploring  p(n) := the n-th prime number :
## 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))))
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))
# }