【アルゴリズム】ミラー・ラビン素性検査

def modPow a, b, m
  v = 1
  p = a % m
  while b > 0 do
    v = (v * p) % m if (b & 1) != 0
    p = (p * p) % m
    b >>= 1
  end
  v
end

def witness a, n
  n1 = n - 1
  s2 = n1 & -n1
  x = modPow a, n1 / s2, n
  return false if x == 1 || x == n1
  while s2 > 1 do
    x = (x * x) % n
    return true if x == 1
    return false if x == n1
    s2 >>= 1
  end
  true
end

def probably_prime? n, k
  # http://en.wikipedia.org/wiki/Miller-Rabin_primality_test
  # n, an integer to be tested for primality
  # k, a parameter that determines the accuracy of the test
  return true if n == 2 || n == 3
  return false if n < 2 || n % 2 == 0
  k.downto(1) do
    return false if witness rand(n - 3) + 2, n
  end
  true
end

# http://www.spoj.pl/problems/PON/
gets.to_i.downto(1) do
  puts probably_prime?(gets.to_i, 1) ? 'YES' : 'NO'
end

01:  #include <stdio.h>
02: #include <stdlib.h>
03: #include <time.h>
04: 
05: typedef unsigned long long U8;
06:  typedef int bool;
07:  
08:  const bool true = 1;
09:  const bool false = 0;
10:  
11:  U8 modMultiply(U8 a, U8 b, U8 m)
12:  {
13:    return a * b % m;
14:  }
15:  
16:  U8 modPow(U8 a, U8 b, U8 m)
17:  {
18:    U8 v = 1;
19:    for (U8 p = a % m; b > 0; b >>= 1, p = modMultiply(p, p, m))
20:      if (b & 1) v = modMultiply(v, p, m);
21:    return v;
22:  }
23:  
24:  bool witness(U8 a, U8 n)
25:  {
26:    U8 n1 = n - 1, s2 = n1 & -n1, x = modPow(a, n1 / s2, n);
27:    if (x == 1 || x == n1) return false;
28:    for (; s2 > 1; s2 >>= 1)
29:    {
30:      x = modMultiply(x, x, n);
31:      if (x == 1) return true;
32:      if (x == n1) return false;
33:    }
34:    return true;
35:  }
36:  
37:  U8 random(U8 high)
38:  {
39:    // http://www.cppreference.com/wiki/c/other/rand
40:   return (U8)(high * (rand() / (double)RAND_MAX));
41:  }
42:  
43:  // http://en.wikipedia.org/wiki/Miller-Rabin_primality_test
44: // n, an integer to be tested for primality
45: // k, a parameter that determines the accuracy of the test
46: bool probablyPrime(U8 n, int k)
47:  {
48:    if (n == 2 || n == 3) return 1;
49:    if (n < 2 || n % 2 == 0) return 0;
50:    while (k-- > 0) if (witness(random(n - 3) + 2, n)) return false;
51:    return true;
52:  }
53:  
54:  // http://www.spoj.pl/problems/PON/
55: int main()
56:  {
57:    srand(time(NULL));
58:    int t;
59:    scanf("%d", &t);
60:    while (t-- > 0)
61:    {
62:      U8 n;
63:      scanf("%lu", &n);
64:      puts(probablyPrime(n, 3) ? "YES" : "NO");
65:    }
66:    return 0;
67:  }