prime.MillerRabin

prime.MillerRabin.rs

struct Rand(u64);
impl Rand {
    fn next(&mut self) -> u64 {
        let a = 1000000007;
        let b = 127;
        let m = 1 << 20;
        self.0 = (self.0 * a + b) % m;
        self.0
    }
}

fn powmod(a: u64, k: u64, m: u64) -> u64 {
    if k == 0 {
        1
    } else {
        let t = powmod(a, k / 2, m);
        let t2 = (t * t) % m;
        (t2 * (if a % 2 == 0 { 1 } else { a })) % m
    }
}

// is_prime?
fn mrtest(n: u64) -> bool {
    if n < 2 {
        false
    } else if n < 4 {
        true
    } else if n % 2 == 0 {
        false
    } else {
        let mut rand = Rand(0);
        let mut d = n - 1;
        while d % 2 == 0 { d /= 2 }
        for _ in 0..100 {
            let a = rand.next() % (n - 1) + 1;
            let mut y = powmod(a, d, n);
            let mut t = d;
            while t != n - 1 && y != 1 && y != n - 1 {
                y = (y * y) % n;
                t <<= 1;
            }
            if y != n - 1 && t % 2 == 0 {
                return false;
            }
        }
        true
    }
}

fn main() {
    const M: u64 = 127;
    let bl = mrtest(M);
    println!("{}", if bl { "YES" } else { "NO" });
}

prime.MillerRabin.cc

namespace Prime {
  template<typename Int>
  Int powmod(Int a, Int k, Int m) {
    if (k == 0) return 1;
    Int t = powmod(a, k/2, m);
    t = (t * t) % m;
    if (k&1) t = (t * a) % m;
    return t;
  }

  // Miller-Rabin test
  template<typename Int=int>
  bool test(Int n) {
    if (n<2) return false;
    if (n<4) return true;
    if (!(n&1)) return false;

    random_device dev;
    mt19937 rand(dev());
    Int d = n-1;
    while (!(d&1)) d/=2;
    rep (_, 100) {
      Int a = rand()%(n-1) + 1;
      Int y = powmod<Int>(a,d,n);
      Int t = d;
      while (t!=n-1 and y!=1 and y!=n-1) {
        y = (y * y) % n;
        t <<= 1;
      }
      if (y!=n-1 and !(t&1)) return false;
    }
    return true;
  }
}
bool isprime(int n) { // naiive
  if (n<2) return false;
  for (int i = 2; i*i <= n; ++i)
    if (n % i == 0) return false;
  return true;
}

int main() {
  const int M = 10000;
  for (int i = 0; i < M; ++i) {
    if(isprime(i) != Prime::test(i)) {
      cout << i << ' ' << isprime(i) << ' ' << Prime::test(i) << endl;
      return 1;
    }
  }
  return 0;
}