12.2 Solution

12.2.1 Trial Division

We begin by implementing trial division algorithm. The above is a ‘brute force’ approach, namely we divide our candidate number, n, by a series of integers (remember, a prime number is an integer greater than 1 that is divisible only by 1 and itself). The method has a small improvement, i.e. the divisors are in the range from 2 up to \(\sqrt{n}\).

function isPrime(n::Int)::Bool
    if n < 4
        return n < 2 ? false : true
    end
    upTo::Int = sqrt(n) |> ceil |> Int
    for i in 2:upTo
        if n % i == 0
            return false
        end
    end
    return true
end

isPrime (generic function with 1 method)

The function is pretty straightforward. The first two prime numbers are 2 and 3. Hence, the if n < 4 block that uses a ternary expression. Next, per algorithm we establish the maximum value of the divisor (upto). To that end we calculate sqrt(n), round it to up to the next whole number (ceil) and convert it to integer (Int). Then, we test if the divisors in the range 2:upTo divide our number n evenly (n % i == 0), if so we return false right away, otherwise the number is prime (return true).

With our isPrime ready, generating a sequence of primes is a breeze.

function getPrimesV1(upTo::Int)::Vec{Int}
    @assert upTo > 1 "upTo must be > 1"
    return filter(isPrime, 2:upTo)
end

getPrimesV1 (generic function with 1 method)

All we had to do was to filter out primes from the sequence of numbers within the desired range 2:upTo. Time for a simple test.

# prime numbers up to 100 from:
# https://en.wikipedia.org/wiki/List_of_prime_numbers
primesFromWiki = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31,
    37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]

getPrimesV1(100) == primesFromWiki

true

Works like a charm.

12.2.2 Sieve of Eratosthenes

Next, we will implement the sieve of Eratosthenes algorithm that finds primes by

… iteratively marking as composite (i.e., not prime) the multiples of each prime, starting with the first prime number, 2.

(see the Wikipedia’s page from the link above).

Let’s get down to it.

function getPrimesV2(upTo::Int)
    @assert upTo > 1 "upTo must be > 1"
    nums::Vec{Int} = 1:upTo |> collect
    isPrimeTests::Vec{Bool} = ones(Bool, upTo)
    isPrimeTests[1] = false # first prime is: 2
    for num in nums
        if isPrimeTests[num]
            # numMultiples - local variable visible only in for loop
            numMultiples = (num*2):num:upTo # multiples of num
            isPrimeTests[numMultiples] .= false
        end
    end
    return nums[isPrimeTests]
end

getPrimesV2 (generic function with 1 method)

This time we begin by defining two vectors nums which contains all the examined integers from a given range and isPrimeTests that contains indication of whether a number in nums is prime. Initially, all the numbers are considered to be primes (ones(Bool, upTo)) except for 1 (isPrimeTests[1] = false). Then, we examine every number (num) out of the candidates (nums). If a number is prime (if isPrimeTests[num]) we set all its multiples (numMultiples) as not primes (false) using broadcasting (.) with the assignment (=) operator. Once finished we return only the candidates that passed the test (return nums[isPrimeTests]) using boolean indexing of a vector.

Time for a simple test.

getPrimesV2(100) == primesFromWiki

true

The test passed. Before we close this chapter let’s see if the two functions give identical results.

getPrimesV1(1000) == getPrimesV2(1000)

true

Yes, they do.