Willans' formula in Rust

Question

I implemented Willans' formula to calculate the nth prime number in Rust:

use std::f64::consts::PI;
pub trait Factorial {
 fn factorial(&self) -> u128;
}
impl Factorial for u8 {
 fn factorial(&self) -> u128 {
 let mut x = *self as u128;
 let mut product: u128 = 1;
 loop {
 if x <= 1 {
 break;
 }
 product *= x;
 x -= 1;
 }
 product
 }
}
fn willans(n: u8) -> u32 {
 1 + (1..(2_u8.pow(n.into()) + 1))
 .map(|i| {
 (n as f64
 / (1..(i + 1))
 .map(|j| {
 (PI * ((j - 1).factorial() + 1) as f64 / j as f64)
 .cos()
 .powf(2.0)
 .floor()
 })
 .sum::<f64>())
 .powf(1.0 / n as f64)
 .floor()
 })
 .sum::<f64>()
 .floor() as u32
}
fn main() {
 for n in 1..6 {
 println!("{}: {}", n, willans(n));
 }
}

Due to the suboptimal nature of the formula, the implementation is currently limited to the first five primes, since only the factorials occurring in the formula up to this input will fit into a u128. I'd like to have feedback on the implementation of the formula in general and would also appreciate suggestions on how to extend the range of processable inputs.

• Toby Speight: Does Rust have a bignum library of some kind? That would be the obvious way to increase the range.
• Richard Neumann: Thanks for the suggestion. I played around with num but have not yet successfully found an implementation of .cos() for it.

Answer 1

factorial should be a free function

The trait is superfluous, here, it'd be simpler to just implement:

fn factorial(n: u8) -> u128 { ... }

Check the input range to avoid overflow

In the absence of a big number library, you can only compute the factorial of relatively few numbers before overflowing a u128: it would be nice to assert sooner (rather than later) and with a clearer error message.

A quick Python script reveals that factorial(34) will work, but factorial(35) will overflow.

fn factorial(n: u8) -> u128 {
 assert!(n < 35, "The factorial of {n} overflows a u128, n should not exceed 34");
 ...
}

The manual loop in factorial is unnecessary

First of all, the pattern:

loop {
 if <condition> {
 break;
 }
 ...
}

Can be advantageously replaced with a while loop with no loss of generality:

while <condition> {
 ...
}

Secondly, you're really just iterating over a range here, so you may as well use a range in the first place rather than incrementing/decrementing yourself:

fn factorial(n: u8) -> u128 {
 assert!(n < 35, "The factorial of {n} overflows a u128, n should not exceed 34");
 let mut result = 1;
 for i in 2..=n {
 result *= i.into();
 }
 result
}

Although you may also prefer a functional version:

fn factorial(n: u8) -> u128 {
 assert!(n < 35, "The factorial of {n} overflows a u128, n should not exceed 34");
 (2..=n).fold(1u128, |(acc, i)| acc * (i as u128))
}

Use inclusive ranges, unless performance is measured to suffer

Instead of writing x..(y + 1), you can just write x..=y.

In tight loops, x..(y+1) may occasionally perform better. In your case, however, the loop bodies are large and complex enough that it seems unlikely this will be an issue.

Only switch back to x..(y + 1) if a profiler hints at it being a problem.

Split up large expressions

The definition of the willans function is somewhat painful to read, as it ends up being heavily nested.

The core functions (the lambdas in map(|i/j| ...) can be extracted for clarity:

fn willans(n: u8) -> u32 {
 fn willans_inner(i: u8, j: u8) -> f64 {
 (PI * (factorial(j - 1) + 1) as f64 / j as f64)
 .cos()
 .powf(2.0)
 .floor()
 }
 fn willans_outer(n: u8, i: u8) -> f64 {
 let divisor = (1..=i).map(|j| willans_inner(i, j)).sum::<f64>();
 (n as f64 / divisor).pow(1.0 / n as f64).floor()
 }
 1 + (1..(2u8.pow(n.into()) + 1))
 .map(|i| willans_outer(n, i))
 .sum::<f64>()
 .floor() as u32
}

Numerical accuracy

The conversion from u128 to f64 leaves me uncomfortable. A f64 only has 53 bits of mantissa (accounting for the leading 1), meaning that only the top 53 bits of the u128 will be taken into account, and everything else will be thrown away.

I'm not familiar with the formula, so maybe it's accounted for. If that's the case, a prominent comment next to the conversion would be most welcome.

Beyond the 5th prime

In order to calculate the 6th prime, and beyond, you'll need to use a so-called BigNum library.

There are various bindings for GMP which are available, and GMP is old enough itself that it should feature all the pieces of functionality that you need.