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

标签：

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

算法描述

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

根据上述的费马定理，每当 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 算法能够以很高的概率来检验一个很大的数是否素数。该算法描述如下：

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

将以 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 程序如下：

在该网站提交，结果是“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 语言中并没有内置的大整数类型。

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

标签：

原文地址：http://blog.csdn.net/u013948191/article/details/47270979