银河

SKYIV STUDIO
  博客园 :: 首页 :: 博问 :: 闪存 :: :: :: 订阅 订阅 :: 管理 ::
【算法】米勒-拉宾素性检验

问题描述

问题来源于 Sphere Online Judge (SPOJ) 网站的的“Prime or Not”题目：

Prime or Not: Given the number, you are to answer the question: "Is it prime?"

Input:
t – the number of test cases, then t test cases follows. [t ≤ 500]
Each line contains one integer: N [2 ≤ N ≤ 263-1]

Output: For each test case output string "YES" if given number is prime and "NO" otherwise.

Time limit: 21s

Source limit: 5,000 Bytes

这道题目要求我们判断大约五百个给定的数是不是素数，其中每个数不超过 263-1 。时间限制是 21 秒。如果通过试除法进行因子分解，肯定会超时。如果使用筛法，肯定会造成内存超限。

算法描述

首先，让我们来看看跟素数有关的费马定理：

如果 p 是素数，则对于所有整数 a，a p ≡ a (mod p)。

根据上述的费马定理，每当 p 是素数以及 a 不是 p 的倍数时，我们有 a p - 1 ≡ 1 (mod p) 。而且有有效的方法计算 a n - 1 mod n ，且只需要 O(log n) 个模 n 乘法运算。因此我们可以确定，当这个关系不成立时，n 不是素数。对于一个给定的数的非素性来说，费马定理是一个强有力的检验。当 n 不是素数时，总是有可能来求 a < n 的一个值，使得 a n - 1 ≠ 1 (mod n) 。事实上，经验证明，这样一个值几乎总能非常快地求出。有某些稀少的 n 值，它们经常使得 a n - 1 ≡ 1 (mod n) ，但在此情况下，n 有小于 n 1/3 的因子。

我没有找到确定一个很大的数是否素数的有效的算法。但是 Miller-Rabin primatlity test 算法能够以很高的概率来检验一个很大的数是否素数。该算法描述如下：

Input: n > 3, an odd integer to be tested for primality;
Input: k, a parameter that determines the accuracy of the test
Output: composite if n is composite, otherwise probably prime
01: write n − 1 as 2s·d with d odd by factoring powers of 2 from n − 1
02: LOOP: repeat k times:
03:   pick a randomly in the range [2, n − 2]
04:   x ← ad mod n
05:   if x = 1 or x = n − 1 then do next LOOP
06:   for r = 1 .. s − 1
07:     x ← x2 mod n
08:     if x = 1 then return composite
09:     if x = n − 1 then do next LOOP
10:   return composite
11: return probably prime

构成该算法的思想是，如果 a d ≠ 1 (mod n) 以及 n = 1 + 2s · d 是素数，则值序列

a d mod n，a 2d mod n，a 4d mod n，…，a 2s d mod n

将以 1 结束，而且在头一个 1 的前边的值将是 n – 1 (当 p 是素数时，对于 y 2 ≡ 1 (mod p) ，仅有的解是 y ≡ ±1 (mod p)，因为 (y + 1)(y - 1)必须是 p 的倍数)。注意，如果在该序列中出现了 n – 1，则该序列中的下一个值一定是 1。因为：(n – 1)2  ≡  n2 – 2n + 1  ≡  1  (mod n)。在该算法中：

  • 该算法用于判断一个大于 3 的奇数 n 是否素数。参数 k 用于决定 n 是素数的概率。
  • 该算法能够肯定地判断 n 是合数，但是只能说 n 可能是素数。
  • 第 01 行，将 n – 1 分解为 2s·d  的形式，这里 d 是奇数。
  • 第 02 行，将以下步骤(第 03 到 10 行)循环 k 次。
  • 第 03 行，◇在 [2, n - 2] 的范围中独立和随机地选择一个正整数 a 。
  • 第 04 行，◇计算该序列的第一个值：x ← ad mod n 。
  • 第 05 行，◇如果该序列的第一个数是 1 或者 n - 1，符合上述条件，n 可能是素数，转到第 03 行进行一下次循环。
  • 第 06 行，◇循环执行第 07 到 09 行，顺序遍历该序列剩下的 s – 1 个值。
  • 第 07 行，◇◇计算该序列的下一个值：x ← x2 mod n 。
  • 第 08 行，◇◇如果这个值是 1 ，但是前边的值不是 n - 1，不符合上述条件，因此 n 肯定是合数，算法结束。
  • 第 09 行，◇◇如果这个值是 n - 1，因此下一个值一定是 1，符合上述条件，n 可能是素数，转到第 03 行进行下一次循环。
  • 第 10 行，◇发现该序列不是以 1 结束，不符合上述条件，因此 n 肯定是合数，算法结束。
  • 第 11 行，已经对 k 个独立和随机地选择的 a 值进行了检验，因此判断 n 非常有可能是素数，算法结束。

在一次检验中，该算法出错的可能顶多是四分之一。如果我们独立地和随机地选择 a 进行重复检验，一旦此算法报告 n 是合数，我们就可以确信 n 肯定不是素数。但如果此算法重复检验 25 次报告都报告说 n 可能是素数，则我们可以说 n “几乎肯定是素数”。因为这样一个 25 次的检验过程给出关于它的输入的错误信息的概率小于 (1/4)25。这种机会小于 1015 分之一。即使我们以这样一个过程验证了十亿个不同的素数，预料出错的概率仍将小于百万分之一。因此如果真出了错，与其说此算法重复地猜测错，倒不如说由于硬件的失灵或宇宙射线的原因，我们的计算机在它的计算中丢了一位。这样的概率性算法使我们对传统的可靠性标准提出一个问号：我们是否真正需要有素性的严格证明。(以上文字引用自 Donald E. Knuth 所著的《计算机程序设计艺术 第2卷 半数值算法(第3版)》第 359 页“4.5.4 分解素因子”中的“算法P(概率素性检验)”后面的说明) 

Ruby 程序

根据上述算法，编写出相应的 Ruby 程序如下：

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

在该网站提交，结果是“accepted”，运行时间 0.23 秒，内存占用 4.7 MB ，目前在 Ruby 语言中排名第三位。

C 程序

相应的 C 语言程序如下所示：

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:  }

在该网站提交，运行结果是“wrong answer”。这是由于源程序中第 11 到 14 行的 modMultiply 函数运算溢出造成的，虽然该函数的参数和返回值都能够表示题目要求的数值范围(不超过 64 bits 整数的范围)，但是其中的乘法的中间的结果需要 128 bits 的整数才能表达，造成了溢出。而 C 语言中并没有内置的大整数类型。

参考资料

  1. Wikipedia: Fermat’s little theorem
  2. Wikipedia: Miller-Rabin primality test
  3. Wikipedia: Sieve of Eratosthenes
  4. The Art of Computer Programming

算法和数据结构目录

posted on 2010-07-10 21:39  银河  阅读(12066)  评论(16)  编辑  收藏  举报