Faster Sieve of Eratosthenes

Question

This is an implementation of the Sieve of Eratosthenes :

• It takes advantages of the fact that all primes from 5 and above can be written as 6X-1 or 6X+1,

• For better space complexity, it uses a pretty accurate upperbound. Better estimations of the upper bound can be found here. I've observed a very slight increase in performance with this.

---

func eratosthenesSieve(to n: Int) -> [Int] {

    guard 2 <= n else { return [] }

    var composites = Array(repeating: false, count: n + 1)
    var primes: [Int] = []
    let d = Double(n)
    let upperBound = Int((d / log(d)) * (1.0 + 1.2762/log(d)))
    primes.reserveCapacity(upperBound)
    let squareRootN = Int(d.squareRoot())

    //2 and 3
    var p = 2
    let twoOrThree = min(n, 3)
    while p <= twoOrThree {
        primes.append(p)
        var q = p * p
        let step = p * (p - 1)
        while q <= n {
            composites[q] = true
            q += step
        }
        p += 1
    }

    //5 and above
    p += 1
    while p <= squareRootN {
        for i in 0..<2 {
            let nbr = p + 2 * i
            if !composites[nbr] {
                primes.append(nbr)
                var q = nbr * nbr
                var coef = 2 * (i + 1)
                while q <= n  {
                    composites[q] = true
                    q += coef * nbr
                    coef = 6 - coef
                }
            }
        }
        p += 6
    }

    while p <= n {
        for i in 0..<2 {
            let nbr = p + 2 * i
            if nbr <= n && !composites[nbr] {
                primes.append(nbr)
            }
        }
        p += 6
    }

    return primes
}

---

It was inspired by this code by Mr Martin.

Using the same benchmarking code in that answer, adding a fourth fractional digit in the timing results, plus some formatting, here are the results :

        ---------------------------------------------------------------
        |                |       Nbr   |          Time (sec)          | 
        |      Up to     |        of   |------------------------------|
        |                |    Primes   |   Martin's   |      This     |
        |----------------|-------------|------------------------------|
        |        100_000 |        9592 |       0.0008 |        0.0004 |
        |----------------|-------------|--------------|---------------|
        |      1_000_000 |      78_498 |       0.0056 |        0.0026 |
        |----------------|-------------|--------------|---------------|
        |     10_000_000 |     664_579 |       0.1233 |        0.0426 |
        |----------------|-------------|--------------|---------------|
        |    100_000_000 |   5_761_455 |       1.0976 |        0.5089 |
        |----------------|-------------|--------------|---------------|
        |  1_000_000_000 |  50_847_534 |      12.1328 |        5.9759 |
        |----------------|-------------|--------------|---------------|
        | 10_000_000_000 | 455_052_511 |     165.5658 |       84.5477 |
        |----------------|-------------|--------------|---------------|

Using Attabench, here is a visual representation of the performance of both codes while n is less than 2^16:

---

One thing I observe is some elements in the composites array are marked with true multiple times. This is expected (but unwanted) behavior since 6X-1 or 6X+1 aren't all primes.

What I'm looking for is making this Sieve of Eratosthenes quicker. I'm well aware of faster methods of finding primes.

Naming, code clarity, conciseness, consistency, etc, are welcome but are not the main point here.

Answer 1

It takes advantages of the fact that all primes from 5 and above can be written as 6X-1 or 6X+1,

I don't think it does, really. It structures the code around that fact, but to take advantage of it, at a minimum you should replace

    while p <= twoOrThree {
        primes.append(p)
        var q = p * p
        let step = p * (p - 1)
        while q <= n {
            composites[q] = true
            q += step
        }
        p += 1
    }

with

    while p <= twoOrThree {
        primes.append(p)
        p += 1
    }

which in my testing gives a significant speedup.

---

To maximise the advantage, you could reduce composites to only store flags for \$6X \pm 1\$. Proof of concept code (could be tidier):

    var pidx = 1
    p = 5
    while p <= squareRootN {
        if !composites[pidx] {
            primes.append(p)

            var qidx = 3 * pidx * (pidx + 2) + 1 + (pidx & 1)
            let delta = p << 1
            let off = (4 - 2 * (pidx & 1)) * pidx + 1
            while qidx < composites.count {
                composites[qidx - off] = true
                composites[qidx] = true
                qidx += delta
            }
            if qidx - off < composites.count {
                composites[qidx - off] = true
            }
        }

        pidx += 1
        p += 2 + 2 * (pidx & 1)
    }

    while p <= n {
        if !composites[pidx] { primes.append(p) }
        pidx += 1
        p += 2 + 2 * (pidx & 1)
    }

This gives a moderate speedup in my testing.