2
\$\begingroup\$

I implemented the Miller-Rabin prime test in Rust and made a program to generate large primes.

I have also implemented the same program in C and Haskell and the Rust version is the slowest. I am looking for advice on how to improve performance and how to improve my Rust style code.

The code for the program in Rust, Haskell, and C are below is compiled using cargo build --release. The contents of the Rust toml file follow the rust code.

For an input initial number n=10^500 and number of tests k=40, the time taken is 4.5 seconds on my computer (1 sec for C, 2 sec for Haskell) and the answer I get for the next prime is 10^500+961. Each program is run with n and k as the arguments, e.g. cargo run <n> <k>.

Rust code, main.rs, compiled with cargo build --release

use num;
use rand;
use num_bigint::{BigUint, RandBigInt};
use num::FromPrimitive;
use num::{Zero, One};
use std::env;

const TRIAL_DIVISORS : [u32; 167] = [3,   5,   7,  11,  13,  17,  19,  23,  29,  31,  37,  41,
  43,  47,  53,  59,  61,  67,  71,  73,  79,  83,  89,  97, 101,
 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167,
 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239,
 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313,
 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397,
 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467,
 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569,
 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643,
 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733,
 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811, 821, 823,
 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911,
 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997];

fn main() {
    let args: Vec<String> = env::args().collect();
    let n = args[1].parse::<BigUint>().expect("Error reading bignum.");
    let ntests = args[2].parse::<usize>().expect("Error reading ntests.");

    let p = find_prime(n, ntests);

    println!("{}", p);


}

fn find_prime(mut n : BigUint, ntests : usize) -> BigUint {
    // If the input is even, it should be made odd.
    if &n % 2u32 == BigUint::zero() {
        n += 1u32;
    }

    let two : BigUint = BigUint::from_u32(2).unwrap();

    while !mr_isprime(&n, &ntests) {
        n += &two;
    }

    n
}

fn mr_isprime(n : &BigUint, ntests : &usize) -> bool {


    for i in TRIAL_DIVISORS.iter() {
        if n % i  == BigUint::zero() {
            if n==&(BigUint::from_u32(*i).unwrap()) {
                return true
            }
            return false
        }
    }

    let (d,r) = decompose(n);
    let mut rng = rand::thread_rng();
    let two : BigUint = BigUint::from_u32(2).unwrap();

    for _ in 0..*ntests {
        let a: BigUint = rng.gen_biguint_range(&two,&(n-2u16));
        if trial_composite(n, &d, &r, &a) {
            return false;
        }
    }
    true
}


fn trial_composite(n: &BigUint, d: &BigUint,
            r: &usize, a: &BigUint) -> bool {
    let mut x = a.modpow(&d, &n);
    if (x==BigUint::one()) || (x==(n-1u32)) {
        return false;
    }
    let two = BigUint::from_u32(2).unwrap();
    for i in 0..(r-1) {
        let e = d*( &two << i);
        x = a.modpow(&e,n);
        if n - 1u32 == x {
            return false;
        }
    }

    true
}

fn decompose(n : &BigUint) -> (BigUint, usize) {
    // Split number such that
    // n = d*2^r + 1
    let mut d = n - 1u32;
    let mut r : usize = 0;
    while (&d % 2u32).is_zero() {
        r += 1;
        d /= 2u32;
    }
    (d, r)

}

Cargo.toml

[package]
name = "miller_rabin"
version = "0.1.0"
authors = [""]
edition = "2018"

[dependencies]
num = "0.2.0"
num-bigint = { version = "0.2.2", features = ["rand"] }
rand = "0.6.5"

Haskell code, miller_rabin.hs. Compiled with ghc -threaded -O2 miller_rabin.hs -o miller_rabin

module Main where

import System.Random (StdGen, getStdGen, randomRs)
import System.Environment (getArgs)


trial_divisors = [ 3,   5,   7,  11,  13,  17,  19,  23,  29,  31,  37,  41,
  43,  47,  53,  59,  61,  67,  71,  73,  79,  83,  89,  97, 101,
 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167,
 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239,
 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313,
 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397,
 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467,
 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569,
 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643,
 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733,
 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811, 821, 823,
 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911,
 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997]

main :: IO ()
main = do
    g <- getStdGen
    args <- getArgs 
    let number = read . head $ args :: Integer
        ntests = read $ args !! 1 :: Int
        p = find_prime g ntests number
    putStrLn . show $ p



miller_rabin :: StdGen -> Int -> Integer -> Bool
miller_rabin g k n = if any (\d -> n `mod` d == 0) trial_divisors
                        then if any (\d -> n == d) trial_divisors 
                            then True
                            else False
                        else all (not . trial_composite n d r) $
                             take k (randomRs (2, n-2) g)    
    where 
        (d, r) = decompose (n-1) 0


trial_composite :: Integer -> Integer -> Integer -> Integer -> Bool
trial_composite n d r a = let x = fastPow a d n in
    if (x == 1) || (x==n-1) 
    then False
    else all ((/=) (n-1)) $ map (\i -> fastPow a (d*(2^i)) n) [0..r-1]

decompose :: Integer -> Integer -> (Integer, Integer)
decompose d r 
    | d `mod` 2 == 0 = decompose (d `div` 2) (r+1)
    | otherwise = (d, r) 


fastPow :: Integer -> Integer -> Integer -> Integer
fastPow base 1 m = mod base m
fastPow base pow m | even pow = mod ((fastPow base (div pow 2) m) ^ 2) m
                   | odd  pow = mod ((fastPow base (div (pow-1) 2) m) ^ 2 * base) m




find_prime :: StdGen -> Int -> Integer -> Integer
find_prime g k n 
    | even n    = find_prime_odd g k (n+1)
    | otherwise = find_prime_odd g k n
  where
    find_prime_odd g k n = case miller_rabin g k n of
                            True  -> n
                            False -> find_prime g k (n+2)

C code, miller_rabin.c, compiled with gcc -O2 miller_rabin_find_prime.c -o miller_rabin -lgmp

#include <gmp.h>
#include <time.h>
#include <stdlib.h>
#include <stdio.h>
#include <assert.h>
#include <stdbool.h>
#include <math.h>

#define NDIVISORS 167

const int trial_divisors[NDIVISORS] = {3,   5,   7,  11,  13,  17,  19,  23,  29,  31,  37,  41,
  43,  47,  53,  59,  61,  67,  71,  73,  79,  83,  89,  97, 101,
 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167,
 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239,
 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313,
 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397,
 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467,
 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569,
 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643,
 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733,
 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811, 821, 823,
 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911,
 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997};


bool trial_composite(mpz_t n, mpz_t d, unsigned long int r, 
                    mpz_t a, mpz_t x, mpz_t tmp) {
  mpz_powm(x, a, d, n);
  mpz_sub_ui(tmp,n,1);
  if ((mpz_cmp_ui(x,1)==0) || (mpz_cmp(x,tmp)==0)) {
    return false;
  }
  for (unsigned long int i=0; i<r; i++) {
    mpz_mul_2exp(tmp, d, i);
    mpz_powm(x, a, tmp, n);
    mpz_sub_ui(tmp, n, 1);
    if (mpz_cmp(x,tmp)==0) {
      return false;
    }
  }
  return true;
}

bool mr_test(mpz_t n, int numtests) {


  mpz_t d;
  unsigned long int r=0;
  mpz_init(d);
  mpz_t tmp;
  mpz_init(tmp);
  mpz_sub_ui(d,n,1);

  // Decompose into d*2^r + 1 = n
  while (mpz_divisible_ui_p(d,2)) {
    mpz_fdiv_q_2exp(d,d,1);
    r++;
  }


  // Trial division
  for (int i=0; i<NDIVISORS; i++) {
    if (mpz_divisible_ui_p(n,trial_divisors[i])) {
      if (mpz_cmp_ui(n, trial_divisors[i])==0) {
        return true;
      }
      return false;
    }
  }




  gmp_randstate_t rstate;
  gmp_randinit_default(rstate);

  mpz_t x;

  mpz_init(x);
  mpz_t a;
  mpz_init(a);

  for (int k=0; k<numtests; k++) {
    mpz_sub_ui(tmp, n, 4);
    mpz_urandomm(a, rstate, tmp);
    mpz_add_ui(a,a,2);
    if (trial_composite(n,d,r,a,x,tmp)) {
      return false;
    }

  }

  mpz_clear(d);
  mpz_clear(a);
  mpz_clear(x);
  mpz_clear(tmp);

  return true;
}


int main(int argc, char *argv[]){
  srand(time(NULL));


  mpz_t n;
  int flag;
  mpz_init(n);
  mpz_set_ui(n,0);
  flag = mpz_set_str(n, argv[1], 10);
  assert(flag==0);

  int k = atoi(argv[2]);

  if (mpz_divisible_ui_p(n,2)) {
    mpz_add_ui(n,n,1);
  }

  bool p;
  while (!mr_test(n,k)) {
    mpz_add_ui(n,n,2);
  }

  mpz_out_str(stdout, 10, n);
  printf("\n");

  mpz_clear(n);
    return 0;
}
\$\endgroup\$
  • \$\begingroup\$ I can add the C and Haskell code, I can do that in a few hours (wanted to focus on the rust code for this question) \$\endgroup\$ – user668074 Jun 25 at 12:09
  • \$\begingroup\$ Note that the other languages can imply comparative-review. Make sure to mark the variant you want to get reviewed. \$\endgroup\$ – Zeta Jun 25 at 13:51
  • \$\begingroup\$ @Stargateur Haskell's Integer is an arbitrary large integer, there's no external library required (although GHC usually uses GMP in its implementation). \$\endgroup\$ – Zeta Jun 25 at 13:53
  • \$\begingroup\$ @Stargateur, I used GMP bigint in C and Haskell’s builtin Integer type which has arbitrary size \$\endgroup\$ – user668074 Jun 25 at 14:26
  • 1
    \$\begingroup\$ it's hard to follow without explanation of your code, I don't know how this is suppose to work, but I see some things that could be a problem, you often build a BigUint that you already build before. Also each time your function is call you generate a new rng. \$\endgroup\$ – Stargateur Jun 26 at 4:54
4
\$\begingroup\$

I am looking for advice on how to improve performance and how to improve my Rust style code.

Not so applicable, but some minor C comments

Over specifying array size

trial_divisors[] array size is specified with a constant and initialized with maybe the correct number of initializers. Avoid that maybe.

Instead initialize, then form the size.

const int trial_divisors[] = {
  3,   5,   7,  11,  13,  17,  19,  23,  29,  31,  37,  41,
  43,  47,  53,  59,  61,  67,  71,  73,  79,  83,  89,  97, 101,
  // ...
  919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997
};

#define NDIVISORS (sizeof trial_divisors / sizeof trial_divisors[0])

Add comments

The goal, algorithm, use restrictions, etc of the functions are not obvious - some more light commentary is warranted.

Vertical white-spaces

Blank lines in functions appeared excessive.

\$\endgroup\$
  • 1
    \$\begingroup\$ Thank you for the advise. I was not aware of that method of finding the size of an array. \$\endgroup\$ – user668074 Jun 26 at 4:47
4
\$\begingroup\$

I'm not a Rust user, so I can't say much about style, but I can see an optimisation in trial_composite:

        x = a.modpow(&e,n);

You already have \$x^{e/2}\$, so all you need to do is square it. I.e. this line should be (suitably corrected to compile)

        x = x.modpow(2u32,n);

or

        x = x * x % n;

There may also be a further, minor, optimisation in pulling out a local constant for n - 1u32 so as to avoid having to do the subtraction each time round the loop.


There is one point of style which I think is pretty universal:

            if n==&(BigUint::from_u32(*i).unwrap()) {
                return true
            }
            return false

is an overly complicated way of writing

            return n==&(BigUint::from_u32(*i).unwrap())
\$\endgroup\$
  • 1
    \$\begingroup\$ So changing the way of calculating the exponential to x = x * x % n; had a significant improvement on the performance. From 4.5 seconds to 1.5 seconds. Can't remember why I used this method, but I first used it in the C version (C GMP has a specific function for a*d^e). \$\endgroup\$ – user668074 Jun 26 at 9:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.