primes: Find all Primes Less Than n

In mmaechler/sfsmisc: Utilities from 'Seminar fuer Statistik' ETH Zurich

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.

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.

Value

numeric vector of all prime numbers \le 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 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))

mmaechler/sfsmisc documentation built on Feb. 28, 2024, 4:18 a.m.