【算法】米勒-拉宾素性检验

来源：互联网发布：拳皇mugen出招简化软件编辑：程序博客网时间：2024/04/28 20:40

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 ≤ 2⁶³-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

这道题目要求我们判断大约五百个给定的数是不是素数，其中每个数不超过 2⁶³-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 2^s·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 ← a^d mod n
05: if x = 1 or x = n − 1 then do next LOOP
06: for r = 1 .. s − 1
07: x ← x² 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 + 2^s· d 是素数，则值序列

a ^d mod n，a ^2d mod n，a ^4d mod n，…，a ^{2^sd} mod n

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

该算法用于判断一个大于 3 的奇数 n 是否素数。参数 k 用于决定 n 是素数的概率。
该算法能够肯定地判断 n 是合数，但是只能说 n 可能是素数。
第 01 行，将 n – 1 分解为 2^s·d 的形式，这里 d 是奇数。
第 02 行，将以下步骤(第 03 到 10 行)循环 k 次。
第 03 行，◇在 [2, n - 2] 的范围中独立和随机地选择一个正整数 a 。
第 04 行，◇计算该序列的第一个值：x ← a^d mod n 。
第 05 行，◇如果该序列的第一个数是 1 或者 n - 1，符合上述条件，n 可能是素数，转到第 03 行进行一下次循环。
第 06 行，◇循环执行第 07 到 09 行，顺序遍历该序列剩下的 s – 1 个值。
第 07 行，◇◇计算该序列的下一个值：x ← x² 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)²⁵。这种机会小于 10¹⁵ 分之一。即使我们以这样一个过程验证了十亿个不同的素数，预料出错的概率仍将小于百万分之一。因此如果真出了错，与其说此算法重复地猜测错，倒不如说由于硬件的失灵或宇宙射线的原因，我们的计算机在它的计算中丢了一位。这样的概率性算法使我们对传统的可靠性标准提出一个问号：我们是否真正需要有素性的严格证明。(以上文字引用自 Donald E. Knuth 所著的《计算机程序设计艺术第2卷半数值算法(第3版)》第 359 页“4.5.4 分解素因子”中的“算法P(概率素性检验)”后面的说明)

Ruby 程序

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
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/rand40:    return (U8)(high * (rand() / (double)RAND_MAX));41:  }42:  43:  // http://en.wikipedia.org/wiki/Miller-Rabin_primality_test44:  // n, an integer to be tested for primality45:  // k, a parameter that determines the accuracy of the test46:  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 语言中并没有内置的大整数类型。

0 0