Miller_Rabbin算法判断大素数

普通的素数测试我们有O(√ n)的试除算法。事实上，我们有O(s*log³n)的算法。

下面就介绍一下Miller_Rabbin算法思想：

定理一：假如p是质数，且(a,p)=1，那么a^(p-1)≡1(mod p)。即假如p是质数，且a,p互质，那么a的(p-1)次方除以p的余数恒等于1。（费马小定理）

定理二：如果p是一个素数，那么对于x(0<x<p)，若x^2 mod p 等于1，则x=1或p-1。

它利用了费马小定理，即：如果p是质数，且a，p互质，那么a^(p-1) mod p恒等于1。也就是对于所有小于p的正整数a来说都应该复合a^(p-1) mod p恒等于1。那么根据逆否命题，对于一个p，我们只要举出一个a（a<p）不符合这个恒等式（就是a^(p-1) mod p恒等于1式子），则可判定p不是素数。Miller-rabin算法就是多次用不同的a来尝试p是否为素数。

但是每次尝试过程中还做了一个优化操作，以提高用少量的a检测出p不是素数的概率。这个优化叫做二次探测。它是根据一个定理：如果p是一个素数，那么对于x(0<x<p)，若x^2 mod p 等于1，则x=1或p-1。逆否命题：如果对于x(0<x<p)，若x^2 mod p 不等于1，则p不是素数。根据这个定理，我们要计算a^(p-1) mod p是否等于1时，可以这样计算，设p-1=(2^t) * k。我们从a^k开始，不断将其平方直到得到a^(p-1)，一旦发现某次平方后mod p等于1了，那么说明符合了二次探测定理的逆否命题使用条件，立即检查x是否等于1或p-1，如果不是则可直接判定p为合数。

一些问题：

1、对于p-1=(2^t) * k这个t的找法，可以记录一下x（x表示(p-1)）的二进制形式下末尾有多少0，比如二进制(1000)B代表的是十进制(8)D,8的二进制数形式下末尾0有3个，那么也就是8=(2^3)*1;再例如二进制(11000)B代表的是十进制(24)D,24的二进制数形式下末尾0有3个，那么也就是24=(2^3)*3

2、对于随机数a的生成，因为我们这个算法判断的是p这个数是不是素数，那么素数和任意一个数都互质，题目又有要求a<p所以就随机生成一个[1,p-1]范围内的数就行

3、这个算法的复杂度是O(s*log³n)，这个s是我们提前定好的，因为上面说了我们用的是费马小定理的逆否命题，但是该定理的逆命题是不一定成立的，但是令人可喜的是大多数情况是成立的。所以我们可以多找几个数a来验证一下p是否是素数，找几个数，这个多少数就是s

可以证明，使用以上两个定理以后，检验s次出错的概率至多为2^(-s)，所以这个算法是很可靠的。

 1 #include<stdio.h>

 2 #include<string.h>

 3 #include<iostream>

 4 #include<algorithm>

 5 #include<queue>

 6 #include<map>

 7 #include<vector>

 8 #include<math.h>

 9 #define mem(a,x) memset(a,x,sizeof(a))

10 using namespace std;

11 typedef long long LL;

12 const int maxn=50005;

13 const int mod=26;

14 const int INF=0x3f3f3f3f;

15 const int Times = 10;

16 const int N = 5500;

17 LL ct, cnt;

18 LL fac[N], num[N];

19 LL gcd(LL a, LL b)  //求两数最大公因子

20 {

21     return b? gcd(b, a % b) : a;

22 }

23 LL multi(LL a, LL b, LL m)  //快速乘

24 {

25     LL ans = 0;

26     a %= m;

27     while(b)

28     {

29         if(b & 1)

30         {

31             ans = (ans + a) % m;

32             b--;

33         }

34         b >>= 1;

35         a = (a + a) % m;

36     }

37     return ans;

38 }

39 LL pow(LL a, LL b, LL m)  //快速幂

40 {

41     LL ans = 1;

42     a %= m;

43     while(b)

44     {

45         if(b & 1)

46         {

47             ans = multi(ans, a, m);

48             b--;

49         }

50         b >>= 1;

51         a = multi(a, a, m);

52     }

53     return ans;

54 }

55 bool Miller_Rabin(LL n)  //判断是不是素数

56 {

57     if(n == 2) return true;

58     if(n < 2 || !(n & 1)) return false;

59     LL m = n - 1;

60     int k = 0;

61     while((m & 1) == 0)

62     {

63         k++;  //这个k就是我们讲的时候的t

64         m >>= 1;  //这个m就是k

65     }

66     for(int i=0; i<Times; i++) //Times就是我们事先定义的s（要找a的个数）

67     {

68         LL a = rand() % (n - 1) + 1;  //找一个[1,n-1]内的任意数

69         LL x = pow(a, m, n);

70         LL y = 0;

71         for(int j=0; j<k; j++)

72         {

73             y = multi(x, x, n);

74             if(y == 1 && x != 1 && x != n - 1) return false;

75             x = y;

76         }

77         if(y != 1) return false;

78     }

79     return true;

80 }

81 int main()

82 {

83     LL n;

84     while(cin>>n)

85     {

86         if(Miller_Rabin(n))

87             printf("是素数\n");

88         else printf("不是素数\n");

89     }

90     return 0;

91 }

巴特西

Miller_Rabbin算法判断大素数

最新文章

热门文章