Mersenne Primes Generator

help
1

Below is Crystal code to generate Mersenne Primes.
It's mathematical structure is based on Prime Generator Theory (PGT).
Given an odd exponent n, it will compute the next Mp prime e.g. n = p.

The beauty of this method is you don't need to test|know a list of primes.
It will find the next Mp prime from any starting odd n exponent value.
The math guarantees the primes between n and Mp prime p don't have to be checked.
If anyone is interested, I can provide the papers explaining the math.

# Compile as: crystal build --release --mcpu native MPgenerateLL.cr
# Current Crystal: 1.14.0; LLVM 18.1.6
# Last update: 2024/11/10

require "big"

module Primes
  module Utils
    # Lucas-Lehmer primality test for Mersenne Prims
    def mp_prime?
      p = self.to_big_i
      s, mp = 4.to_big_i, (1.to_big_i << p) - 1
      (p-2).times { s = mp_mod(((s * s) - 2), p, mp) }
      s.divisible_by? mp
    end

    private def mp_mod(n, p, mp)
      while n > mp
        n = (n >> p) + (n & mp)
      end
      n
    end
  end
end

struct Int; include Primes::Utils end

def mpgen(p, steps = 1u64)
  puts "Starting Mp gen exponent = #{p}"
  puts "Starting Mp gen process at step #{steps}"
  r = (1.to_big_i << p) - 1
  modpnk = modpn = 3.to_big_i << p
  (steps-1).times { modpnk = (modpnk << 2) + modpn }
  loop do
    mp = modpnk + r
    print "\rsteps = #{steps}"
    break if (p = mp.bit_length).mp_prime?
    steps += 1
    modpnk = (modpnk << 2) + modpn
  end
  print "\rsteps = #{steps}\n"
  puts "Next Mp prime = #{p}"
  p
end

def tm; t = Time.monotonic; yield; Time.monotonic - t end

def mpgenerate
  p = (ARGV[0].to_u64 underscore: true)
  steps = ARGV.size > 1 ? (ARGV[1].to_u64 underscore: true) : 1u64
  puts tm { mpgen(p, steps) }; puts
end

mpgenerate

I would like to compare it to a well done Rust version.
The Crystal version uses BigInts for arbitray sized integers math.
I don't know what|if the equivaltent is in Rust (native or crates).

If Rust can do the (s * s) operations faster it will be even more optimized.

As you can see, it's not much code.
I'm hoping a good|seasoned Rust programmers can optimally translate it.

Here's some Crystal test output to use as comparison references.
As shown, it shows the time it takes mp_prime? to primality test these Mp primes.
(Times generated on 'noisy' system. Other things going on while testing.)
My Linux laptop: Lenovo Legion Slim 5, AMD Ryzen 7 8845HS, 5.1GHz, 8C|16T

➜  crystal-projects ./MPgenLL 86241
Starting Mp gen exponent = 86241
Starting Mp gen process at step 1
steps = 1
Next Mp prime = 86243                  # M28 - 51,924 digits
00:00:31.036883127

➜  crystal-projects ./MPgenLL 110501
Starting Mp gen exponent = 110501
Starting Mp gen process at step 1
steps = 1
Next Mp prime = 110503                # M29 - 66,530 digits
00:00:57.786394665

➜  crystal-projects ./MPgenLL 132047
Starting Mp gen exponent = 132047
Starting Mp gen process at step 1
steps = 1
Next Mp prime = 132049                # M30 - 79,502 digits
00:01:29.874204215

➜  crystal-projects ./MPgenLL 216089
Starting Mp gen exponent = 216089
Starting Mp gen process at step 1
steps = 1
Next Mp prime = 216091                # M31 - 130,100 digits
00:03:23.980284348

➜  crystal-projects ./MPgenLL 756837
Starting Mp gen exponent = 756837
Starting Mp gen process at step 1
steps = 1
Next Mp prime = 756839                # M32 - 455,663 digits
00:53:11.942810881

➜  crystal-projects ./MPgenLL 859431
Starting Mp gen exponent = 859431
Starting Mp gen process at step 1
steps = 1
Next Mp prime = 859433                # M33 - 517,430 digits
01:11:24.882317036

➜  crystal-projects ./MPgenLL 1257785
Starting Mp gen exponent = 1257785
Starting Mp gen process at step 1
steps = 1
Next Mp prime = 1257787               # M34 - 757,263 digits
02:48:17.749699153

Here's a list of all the known Mp, digits sizes, et al.

2

Rust doesn't have a native (as in language-defined) arbitrary precision integer type, so you will be reaching for a crate. bigdecimal seems to be the most downloaded, but I don't know that for certain. And I don't know if downloads are a good metric for "most optimized", "easiest to use", "most compatible", or any other criteria you might have for choosing a type.

This is implementation-dependent. The operator is implemented with the std::ops::Mul trait.

For what purpose? Microbenchmarks are notoriously misleading. You will need to run benchmarks in the domain and with the data that you intend to target to get a useful comparison. Presumably you can't just use the provided Crystal code, which is why you are asking for Rust. But why compare it to Crystal in the first place?

One observation I can make right away is that the provided Crystal code uses the Lucas-Lehmer primality test, which has worse time complexity than Miller-Rabin. And its Wikipedia page claims that the time complexity can be reduced with the Harvey-Hoeven algorithm. So, even the Crystal code could still be optimized, in theory. (To what degree, I do not know.)

So, why was the Crystal code provided, and why do you want it translated to Rust? An exact translation could end up with comparable performance. Or it could hypothetically perform better by adjusting the algorithm (not limited to the primality test). But this still misses the bigger picture that microbenchmarks measure the code outside of the domain in which it will be used in practice.

5 Likes
3

It looks like Crystal uses gmp under the hood so the Rust equivalent would use the rug crate.

4

Found an implementation on Rosetta Code.

extern crate rug;
extern crate primal;

use rug::Integer;
use rug::ops::Pow;
use std::thread::spawn;

fn is_mersenne (p : usize) {
    let p = p as u32;
    let mut m = Integer::from(1);
    m = m << p;  
    m = Integer::from(&m - 1);
    let mut flag1 = false;
    for k in 1..10_000 {
        let mut flag2 = false;
        let mut div : u32 = 2*k*p + 1;
        if &div >= &m {break; }
        for j in [3,5,7,11,13,17,19,23,29,31,37].iter() {
            if div % j == 0 {
                flag2 = true;
                break;
            }   
        }
        if flag2 == true {continue;}
        if div % 8 != 1 && div % 8 != 7 { continue; }
        if m.is_divisible_u(div) { 
            flag1 = true;
            break;
        }
    }
    if flag1 == true {return ()}
    let mut s = Integer::from(4);
    let two = Integer::from(2);
    for _i in 2..p {
		let mut sqr = s.pow(2);
		s = Integer::from(&Integer::from(&sqr & &m) + &Integer::from(&sqr >> p));
		if &s >= &m {s = s - &m}
		s = Integer::from(&s - &two);
    }
	if s == 0 {println!("Mersenne : {}",p);} 
}

fn main () {
    println!("Mersenne : 2");
    let limit = 11_214;
    let mut thread_handles = vec![];
    for p in primal::Primes::all().take_while(|p| *p < limit) {
        thread_handles.push(spawn(move || is_mersenne(p))); 
    }
    for handle in thread_handles {
        handle.join().unwrap();
    }
}