现在的位置: 首页 > 综合 > 正文

MillerRobin(概率测素数)

2013年01月21日 ⁄ 综合 ⁄ 共 1667字 ⁄ 字号 评论关闭

学习知识,是一个很享受的过程。

尤其是那些自己以前听没听过,但是却又是非常神奇的一种思想的时候。

今天上午学习了,快速概率测素数的算法------MillerRobin(),适用于测试单个素数,出错概率比计算机本身出错的概率还要低

为(1/4)^(s),一般s取50就可以认为是准确测试出了。

算法是基于费马小定理(format),二次探测定理(x*x % p == 1 ,若P为素数,则x的解只能是x = 1或者x = p - 1)加上迭代乘法判断的Miller算法

共同构成的。我觉得我要讲起来证明必然会贻笑大方。看下面别人给出的证明及思想吧。


费尔马小定理:如果p是一个素数,且0<a<p,则a^(p-1)%p=1.

            利用费尔马小定理,对于给定的整数n,可以设计素数判定算法,通过 计算d=a^(n-1)%n来判断n的素性,当d!=1时,n肯定不是素数,当d=1时,n   很可能是素数.

 

二次探测定理:如果p是一个素数,且0<x<p,则方程x^2%p=1的解为:x=1或    x=p-1.

            利用二次探测定理,可以再利用费尔马小定理计算a^(n-1)%n的过程 中增加对整数n的二次探测,一旦发现违背二次探测条件,即得出n不是素数的结论.

    

    如果n是素数,则(n-1)必是偶数,因此可令(n-1)=m*(2^q),其中m是正奇数( 若n是偶数,则上面的m*(2^q)一定可以分解成一个正奇数乘以2的k次方的形式 ),q是非负整数,考察下面的测试:

    序列:

         a^m%n; a^(2m)%n; a^(4m)%n; …… ;a^(m*2^q)%n

    把上述测试序列叫做Miller测试,关于Miller测试,有下面的定理:

定理:若n是素数,a是小于n的正整数,则n对以a为基的Miller测试,结果为真.

Miller测试进行k次,将合数当成素数处理的错误概率最多不会超过4^(-k).

我只给出我的代码,我觉得我的还是比较简洁的,呵呵,如果有不足,请指出,共同进步,谢谢 0.-!


#include <stdio.h>
#include <string.h>
#include <time.h>
#include <iostream>
#include <string>

using namespace std;

int N;

int witness(int a, int n)//随机生成的a,来检测n的素性 
{
	int ans = 1;
	int t = n - 1;//这里需要注意,你如果没有改变乘方的次数的话,最后的判断就是(ans == a) ? 0 : 1;
				 // 并且还要另外开辟空间来存储开始的a,比较麻烦,所以就这样了; 
	int x;
	while (t)
	{
		if (t & 1)
		{
			ans = (long long int)ans * a % n;
		}
		x = a;//从这里开始就是迭代乘法,验证二次验证定理 
		a = (long long int)a * a % n;//这里就相当于 x*x % m = 1 
		if (a == 1 && x != 1 && x != (n - 1))
		{
			return 1; // 这里需要注意,返回一的话就说明,追踪过程中,出现了不是素数的依据. 
		}
		t >>= 1;
	}
	return (ans == 1) ? 0 : 1;
}
		

int MillerRobin(int n, int s) // 一般s取50就可以避免所有的偶然性了. 
{
	if (n == 2)
	{
		return 1;
	}
	if (n < 2 || !(n & 1))
	{
		return 0;
	}
	int a;
	for (int i = 0; i < s; i++)
	{
		a = (long long int )rand() * (n - 2) / RAND_MAX + 1; //这样生成的随机数就是真正的随机数了 
		if (witness(a, n))
		{
			return 0;
		} 
	}
	return 1;
}

int main()
{
	while (scanf("%d", &N) != EOF)
	{
		if (N == 0)
		{
			break;
		}
		if (MillerRobin(N, 50))
		{
			printf("%d is a prime!\n", N);
		}
		else
		{
			printf("%d is not a prime!\n", N);
		}
	}
//	system("pause");
	return 0;
}

抱歉!评论已关闭.