My Rust script that calculates primes seems to freeze at 16777213

Question

I'm new to the Rust language, and I've been working on this, my first Rust project, to get a better grasp of how the lang works. I've got this code that calculates primes, but it stops every time it reaches the value 16,777,213. It doesn't panic, it doesn't crash, it just freezes. I've found another post where someone tried a similar thing in C++ and has the same problem: Prime number generator stops abruptly at 16777213 | StackOverflow, but nobody seems to answer how to fix it.

Here is my code:

fn half(val: i128) -> bool { // For better precision (to prevent false positives in small numbers)
    let mut factors: Vec<i128> = vec![];
    for i in 1i128..(((val as f32)/2.0)+1.0) as i128 {
        let factor: f32 = (val as f32)/(i as f32);
        if factor == factor.round() {
                let factor:i128 = factor as i128;
                if factor != val {
                    factors.push(factor as i128)
                }
            }
        }
    if factors.is_empty() {
        return true
    } else {false}
}
fn root(val:i128) -> bool {
    let mut factors: Vec<i128> = vec![];
    for i in 1i128..((val as f32).sqrt().round()+1.0) as i128 {
        if val > 2 {
            let factor: f32 = (val as f32)/(i as f32);
            if factor == factor.round() {
                let factor: i128 = factor as i128;
                if factor != val {
                    factors.push(factor)
                }
            }
        }
    }
    if factors.is_empty() {
        return true
    } else {false}
}
fn is_prime(val:i128) -> bool {
    if val > 100 {
        root(val)
    }
    else {
        half(val)
    }
}
fn main() {
    //let mut primes: Vec<i128> = vec![];
    let mut i:i128 = 16_700_000;
    loop {
        let prime: bool = is_prime(i);
        if prime {
            println!("{}",i)
            /*
This was the original code, but I removed the vector so that I'd be able to
determine the exact value at which the script breaks at.
            primes.push(i);
            if primes.len() as i32 == 24 {
                println!("{:?}",primes);
                primes.clear()
            }*/
        }
        i += 1
    }
}

Answer 1

You have to replace all f32 by f64. I will try to explain it more deeply, but now found that this fixing problem.

Firstly to simplify we can assume that if val > 100 is always true and exclude half function from considerations.

Now we can rewrite our program to see differences between f32 and f64 versions:

fn root_64(val:i128) -> bool {
    let mut factors: Vec<i128> = vec![];
    for i in 1i128..((val as f64).sqrt().round()+1.0) as i128 {
        if val > 2 {
            let factor: f64 = (val as f64)/(i as f64);
            if factor == factor.round() {
                let factor: i128 = factor as i128;
                if factor != val {
                    factors.push(factor)
                }
            }
        }
    }
    println!("F64 {} {:?}", val, factors);
    if factors.is_empty() {
        return true
    } else {false}
}

fn root_32(val:i128) -> bool {
    let mut factors: Vec<i128> = vec![];
    for i in 1i128..((val as f32).sqrt().round()+1.0) as i128 {
        if val > 2 {
            let factor: f32 = (val as f32)/(i as f32);
            if factor == factor.round() {
                let factor: i128 = factor as i128;
                if factor != val {
                    factors.push(factor)
                }
            }
        }
    }
    println!("F32 {} {:?}", val, factors);
    if factors.is_empty() {
        return true
    } else {false}
}

fn is_prime_64(val:i128) -> bool {
    if val > 100 {
        return root_64(val)
    };
    return false
}
fn is_prime_32(val:i128) -> bool {
    if val > 100 {
        return root_32(val)
    };
    return false
}
fn main() {
    //let mut primes: Vec<i128> = vec![];
    let mut i:i128 = 16_777_210;
    loop {
        let prime_32: bool = is_prime_32(i);
        let prime_64: bool = is_prime_64(i);
        if prime_32 || prime_64 {
            println!("{} 32:{} 64:{}",i, prime_32, prime_64)
        }
        i += 1;
        if i > 16_777_260 {
            break
        }
    }
}

We can see that initially factors lists are equal, but for 16777217 factors lists starts be different:

F32 16777210 [8388605, 3355442, 1677721]
F64 16777210 [8388605, 3355442, 1677721]

...

16777213 32:true 64:true

...

F32 16777217 [16777216, 8388608, 4194304, 2097152, 1048576, 524288, 262144, 131072, 65536, 32768, 16384, 8192, 4096]
F64 16777217 [172961, 65281, 24929]

We can note that

16777217 / (2^24) = 1.0000000596

but

16777216 / (2^24) = 1

So let's focus computing square from number 16777217.

List of factors found by root_32 method indicates that 16777217 was approximated to 16777216 in comparison

factor == factor.round()

where factor is result of division

(val as f32)/(i as f32)

You can even write test for it

#[cfg(test)]
mod tests {
    #[test]
    fn it_works() {
        let result = 16777217 as f32 / 16777216 as f32;
        assert_eq!(result.round(), 1f32);
        assert_eq!(result, 1f32);
    }
}

If you are interested why 2^24 is crucial for this approximation you can read about f32 structure

https://doc.rust-lang.org/std/primitive.f32.html

especially about

pub const MANTISSA_DIGITS: u32 = 24u32
Number of significant digits in base 2.

General advice

Both float division and modulo division are considered as heavy for CPU. You can read why here:

Why is division more expensive than multiplication?

And you should not use them in context of big numbers and computations. Instead I propose to implement Miller Rabin algorithm, that you can check here:

https://github.com/huonw/primal/blob/master/primal-check/src/is_prime.rs

If you really want to push it to limits you probably have to read some Curtis Cooper publications, learn about algorithms for quantum computers and wait when they will be ready :)

https://arxiv.org/abs/1302.6245