高效素数判断与幂运算-优快云博客

素数一直是个头疼的话题，不过最好还是克服一下。至少知道如何来进行处理。

求a^b

那种老掉牙的for循环就不要写了吧。

这里粘个高效的。

template<typename T>
T _exp(T a, T b)
{
	T odd = 1;
	while (b > 1)
	{
		if (b&1) odd *= a;
		a *= a;
		b >>= 1;
	}
	return a*odd;
}

求a^b%n

其实这两者有共性了。

template<typename T>
T _mod_exp(T a, T b, T n)
{
	T odd = 1;
	while (b)
	{
		if (b&1) odd = (odd * a) %n;
		a = (a*a)%n;
		b >>= 1;
	}
	return odd;
}

费马小定理

直接用维基百科上的东西：
费马小定理是数论中的一个定理：假如a是一个整数，p是一个质数，那么

$a^p \equiv a \pmod{p}$

如果a不是p的倍数，这个定理也可以写成

$a^{p-1} \equiv 1 \pmod{p}$

这个书写方式更加常用。（符号的应用请参见同余。）

费马小定理的逆定理

这个小定理的逆定理是不成立的。但是可以用来筛选素数。

即认为满足

$a^{p-1} \equiv 1 \pmod{p}$

就是素数，当然a在这里可以取很多值。

但是也有一些数偏可以通过这些测试。

二次探测

若 p 为素数，且 0 < x < p，则方程 x * x = 1 ( mod p ) 的解为 x = 1, p - 1 ( mod p )。

事实上 x * x = 1 ( mod p ) 等价于 x * x - 1 = 0 ( mod p )，则 ( x - 1 )( x + 1 ) = 0 ( mod p )。故 p 须整除 x - 1 或 x + 1，由 p 为素数且 0 < x < p，推出 x = 1 ( mod p ) 或 x = p - 1 ( mod p )。

基于费马小定理的随机化测试

基实思想很简单，就取一个a值不行，那么就多取几个a值。直到行为止。ull表示typedef unsigned long long ull;

int _test_base(ull n)
{
	int ret2 = 1, ret3 = 1, ret5 = 1, ret7 = 1;
	int i = 0;
	ull a;
	if (!(n&1) || 0 == n % 3 || 0 == n % 5 || 0 == n % 7) return 0;
	for(i=1;i<8;i++)  //个人以为这里的8还可以取成20~50中的数。
	{ 
		a = 2 + rand()%(n-2);
		if (_mod_exp(a, n - 1, n) != 1) return 0;
	}
	return 1;
}

一般而言，毕竟这个随机化测试虽然概率很大，还是有漏网之鱼。

这里需要用到一个更强的版本

Miller-Rabin素数测试

定义：令 n - 1 = ( 2^s ) * d，其中 s 是非负整数，d 是正奇数。若 a^d =1 ( mod n ) 或 a^((2^r)*d) = -1 ( mod n )，0 <= r <= s - 1，则称 n 通过以 a 为基的 R-M 测试。

定理：若 n 为奇合数，则 n 通过以 a 为基的 R-M 测试的数目最多为 ( n - 1 ) / 4， 1 <= a <= n - 1

不多说了这里给出一个强化版的素数测试方法

#include<stdio.h>
#include<stdlib.h>
#include<string.h>
typedef unsigned long long ull;
 
ull _mod_exp(ull a, ull b, ull n)
{
	ull t = 1, y = a;
	while (b)
	{
		if (b&1) t = (t * y) % n;
		y = (y * y) % n;
		b >>= 1;
	}
	return t;
}
 
//1 for not primer.
int _is_not_primer(ull a, ull k, ull q, ull n)
{
	ull e = 1, i;
	if (1 == _mod_exp(a, q, n)) return 0;
	for (i = 0; i < k; ++i)
	{
		if (n-1 == _mod_exp(a, q*e, n)) return 0;
		e <<= 1;
	}
	return 1;
}
 
//1 for not primer
int _miller_robin_not_primer(ull n)
{
	ull k = 0, q = n - 1;
	if (2 == n || 3 == n || 5 == n || 7 == n) return 0;
	if (!(n&1) || 2 > n) return 1;
	while (!(q&1)) {q >>= 1; ++k;}
 
	if (n < 1373653)
	{
		if (_is_not_primer(2, k, q, n) 
			|| _is_not_primer(3, k, q, n))
			return 1;
	}
	else if (n < (ull)9080191)
	{
		if (_is_not_primer(31, k, q, n) 
			|| _is_not_primer(73, k, q, n))
			return 1;
	}
	else if (n < (ull)4759123141)
	{
		if (_is_not_primer(2, k, q, n) 
			|| _is_not_primer(3, k, q, n)
			|| _is_not_primer(5, k, q, n)
			|| _is_not_primer(11, k, q, n)
			)
			return 1;
	}
	else if (n < (ull)2152302898747)
	{
		if (_is_not_primer(2, k, q, n)
			||_is_not_primer(3, k, q, n)
			||_is_not_primer(5, k, q, n)
			||_is_not_primer(7, k, q, n)
			||_is_not_primer(11, k, q, n)
			)
			return 1;
	}
	else
	{
		if (_is_not_primer(2, k, q, n)
			||_is_not_primer(3, k, q, n)
			||_is_not_primer(5, k, q, n)
			||_is_not_primer(7, k, q, n)
			||_is_not_primer(11, k, q, n)
			||_is_not_primer(31, k, q, n)
			||_is_not_primer(61, k, q, n)
			||_is_not_primer(73, k, q, n)
			)
			return 1;
	}
	return 0;        
}
 
int prime(ull n)
{
	return _miller_robin_not_primer(n) == 0;
}
 
int main(void)
{
	ull n;
	while (scanf("%llu", &n) != EOF && n)
	{
		if (prime(n))
		{
			printf("YES\n");
		}
		else 
		{
			printf("NO\n");
		}
	}
	return 0;
}

值得说明的是，对不同范围的数，采用不同的a值就可以了。节省时间考虑。

通过下面的测试

http://acm.cs.ecnu.edu.cn/problem.php?problemid=1094

有用的链接

http://blog.youkuaiyun.com/hikaliv/article/details/4261948

最后版本的那个思想来自于优快云。我忘了是哪个的了~~声明一下是转的，但代码是我自己照着写的。谢过先。

高手的链接

http://wenku.baidu.com/view/7433bcd9ad51f01dc281f153.html