Day 2 - primal

Relevancy: 1.9 stable

When I start learning a new programming language, I like to code at least several solutions to Project Euler problems. These are very math-oriented and may not be the best introduction to general purpose programming, but it's a start. Anyway, it's just fun to solve them! (...and way more fun to solve them in a fast way and not by brute force.)

A lot of Project Euler problems involve prime numbers in some way. These include finding nth prime, efficient factorization or checking whether some curious number is prime or not. You could of course write these mathematical procedures yourself, which is also an educational activity. But I'm lazy. I set out to find some ready-made code and stumbled upon the primal library by Huon Wilson. Incidentally this was the first external dependency I ever used in a Rust program, long before crates.io (back when it was called slow_primes).

So let's see what's in there, shall we?

Prime sieve

The first thing to do is to create a sieve (see Wikipedia on Sieve of Eratosthenes for a detailed explanation of the algorithm). We need to set an upper bound on the sieve. There's a clever way to estimate that bound (see the docs for estimate_nth_prime) but for simplicity I'll hardcode it now to 10000.

Let's actually check some numbers for primality:

day2.rs

extern crate primal;

use primal::Sieve;

day2.rs

let sieve = Sieve::new(10000);
    let suspect = 5273;
    println!("{} is prime: {}", suspect, sieve.is_prime(suspect));
    let not_a_prime = 1024;
    println!("{} is prime: {}", not_a_prime, sieve.is_prime(not_a_prime));
    let n = 1000;

How about finding 1000th prime number?

day2.rs

let n = 1000;
    match sieve.primes_from(0).nth(n - 1) {
        Some(number) => println!("{}th prime is {}", n, number),
        None => println!("I don't know anything about {}th prime.", n),
    }
    println!("{:?}", sieve.factor(2610));

The primes_from() method returns an iterator over all prime numbers generated by this sieve (2, 3, 5, 7...). Iterators in Rust have a lot of useful methods; the nth() method skips over n initial iterations, returning the nth element (or None if we exhausted the iterator). The argument is zero-based, so to find 1000th prime we need to pass 999 to nth().

Factorization

Factorization is a way to decompose a number into its divisors. For example, 2610 = 2 * 3 * 3 * 5 * 29. Here's how we can find it out with primal API:

day2.rs

println!("{:?}", sieve.factor(2610));

When we run this, we'll get:

$ cargo run
Ok([(2, 1), (3, 2), (5, 1), (29, 1)])

What is this? Let's have a look at the result type of factor():

type Factors = Vec<(usize, usize)>;
fn factor(&self, n: usize) -> Result<Factors, (usize, Factors)>

Looks a bit complicated, but remember the Result type. The Ok variant wraps a vector of pairs of numbers. Each pair contains a prime factor and its exponent (how many times it appears in the factorization). In case of an error we'll get a pair (leftover value, partial factorization).

We can use factorization to find the total number of divisors (including compound ones). This is very important in number theory (although for reasons that are outside the scope of this blog).

Consider the following function:

day2.rs

fn num_divisors(n: usize, primes: &Sieve) -> Option<usize> {
    match primes.factor(n) {
        Ok(factors) => Some(factors.into_iter().fold(1, |acc, (_, x)| acc * (x + 1))),
        Err(_) => None,
    }
}

The trick is to multiply all prime factor exponents, incremented before multiplication. See the explanation at Maths Challenge for the curious. So when we call the function on our 2610 example, we'll get Some(24) as a result.

day2.rs

println!("{:?}", num_divisors(2610, &sieve));

Further reading