Task — 10001st prime number

Write a script to generate the 10001st prime number.

Analysis

The first question to ask is whether the 10001st prime can be represented as an integer in Perl. The largest number which Perl treats as an integer is given by ~1 + 1 and on my computer that is 18446744073709551615 or about 2 * (10**19). Primes get a bit thin on the ground as they get larger, but I think we can safely postulate that we will find 10**5 of them before we hit that limit.

My 'go to' algorithm for primes is the sieve of Eratosthenes. Born in 276 BC, Eratosthenes was the first person to estimate the circumference of the earth - remarkably accurately - while managing the largest library in the world, in Alexandria.

His sieve works like this. Write down all the numbers from 2 to however large you want - say n. Take the leftmost number - 2 - and rub out all its multiples - 2, 4, 6, 8 and so on. Now do it again: the leftmost number remaining is 3, because you rubbed out 2, so now you rub out 3, (not 6 because you already rubbed it out), 9, 15 and so on. You repeat that until the leftmost number is more than the square root of n, and what remains is a list of all the prime numbers less than or equal to n.

So that gets us all the primes less than n, but what we want is the nth prime number. Eratosthenes unfortunately neglected to address that problem, but we can improvise.

After a bit of experimenting I came up with this:

• Make a sieve out of 1-1000
• Count the primes
• If you haven't reached 10001 of them then ...
• Make a sieve of 1-2000

and so on This works fine and gets the answer in 3 seconds on my computer. However if we want, say, the 1000001th prime it takes about 27 seconds, So I improved on it by extending the existing sieve by 1000 rather than making a new one. This makes the code somewhat harder to follow, but it reduced the time for the 10001th prime to under a second and the million and first one to under 7 seconds.

Try it 

#!/usr/bin/perl

# Peter Campbell Smith - 2022-01-04
# PWC 146 task 1

use v5.28;
use warnings;
use strict;

my ($seeking, $prime_index, $from, $to, $test, $factor, $multiple, @not_a_prime,
    @prime, $time, $start);

# initialise
$seeking = 10001;
$time = time;
$prime_index = 0;

# find primes in ranges of 1000, and count them until we get to (or past) $seeking
for ($from = 1; $prime_index <= $seeking; $from += 1000) {
    $to = $from + 999;

    # loop over all the possible factors, ie primes < sqrt($to)
    for $factor (2 .. int(sqrt($to))) {
        next if defined $not_a_prime[$factor];   # already known as not a prime
        
        # mark all the multiples of $factor as non-primes
        $start = int($from / $factor);           # multiples less than $start have already been done
        $start = 2 if $start < 2;
        for ($multiple = $start; $factor * $multiple <= $to; $multiple ++) {
            $not_a_prime[$factor * $multiple] = 1;
        }
    }
    
    # now enumerate the primes in this range
    for $test ($from .. $to) {
        next if $test == 1;                       # 1 is not regarded as a prime
        if (not defined $not_a_prime[$test]) {    # $test is a new prime
            $prime[++ $prime_index] = $test;
        }
    }
}

say qq[Prime no $seeking is $prime[$seeking]\nFound in ] . (time - $time) . qq[ seconds];
say qq[This assumes that that first prime is 2, as per Wikipedia. If you reckon that
the first prime is 1 then prime number $seeking is $prime[$seeking - 1]];

Output

Prime no 10001 is 104743
Found in 0 seconds
This assumes that that first prime is 2, as per Wikipedia. If you reckon that
the first prime is 1 then prime number 10001 is 104729