HDu 2138 How many prime numbers 高效Miller素数测试

题目就是给出一组数，让我们测试其中有多少个是素数。

求素数有测试sqrt(n)个数的方法，有筛子方法，不过对于本题这样的题目来说就都不是高效的。

本题使用Miller Rabin素数测试法，效率奇高，对于不是极其大的整数测试都几乎是常数时间。令人神往的算法啊。

网上有个程序，好像是什么吉林的模板程序，不过我一直没看懂他是什么思路写的，是个AC的程序，不过却是错误的，呵呵，因为程序一直把9当做素数。

于是上网查找了其中原理，自己写了个程序，效率和他的差不多一样，通过时间基本无差别，不过我的思路是按照wiki上的写的。

wiki介绍例子：

Example[edit]

Suppose we wish to determine if n = 221 is prime. We write n ? 1 = 220 as 2²·55, so that we have s = 2 and d = 55. We randomly select
a number a such that a < n, say a = 174. We proceed to compute:

• a^20·d mod n = 174⁵⁵ mod 221 = 47 ≠ 1, n ? 1
• a^21·d mod n = 174¹¹⁰ mod 221 = 220 = n ? 1.

Since 220 ≡ ?1 mod n, either 221 is prime, or 174 is a strong liar for 221. We try another random a, this time choosing a = 137:

• a^20·d mod n = 137⁵⁵ mod 221 = 188 ≠ 1, n ? 1
• a^21·d mod n = 137¹¹⁰ mod 221 = 205 ≠ n ? 1.

Hence 137 is a witness for the compositeness of 221, and 174 was in fact a strong liar. Note that this tells us nothing about the factors of 221 (which are 13 and 17). However, the example with 341 in the next section shows how these calculations can sometimes
produce a factor of n.

出处：http://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test

参照上面例子就可以写程序了，挺难写的一个程序，运用好数论知识和位运算才能写出高效的程序：

代码乃属原创，转载请注明作者：

作者 靖心 http://blog.csdn.net/kenden23/article/details/30478591

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <cmath>
__int64 s, d;

void preCalsd(__int64 N)
{
	s = 0;
	N--;
	while ((N & 1) == 0)
	{
		s++;
		N >>= 1LL;
	}
	d = N;
}

__int64 fastPowMod(__int64 base, __int64 num, __int64 mod)
{
	__int64 ans = 1LL;
	while (num)
	{
		if (num & 1) ans = ans * base % mod;
		base = base * base % mod;
		num >>= 1LL;
	}
	return ans % mod;
}

bool witness(__int64 a, __int64 N)
{
	__int64 num = d, m = -1;
	for (int i = 0; i < s; i++)
	{
		num = d * (1LL << i);
		m = fastPowMod(a, num, N);
		if (m == 1 || m == N-1) return false;
	}
	return true;
}

bool miller_rabin(__int64 N)
{
	if(N == 2)	return true;
	if(N == 1 || ((N&1) == 0))	return false;

	preCalsd(N);
	for(int i = 0; i < 3; i++) //i < 50
	{
		__int64 a = 2 + rand() % (N-2);
		if(witness(a, N)) return false;
	}
	return true;
}

int main()
{
	srand(unsigned(time(NULL)));
	int n,cnt;
	__int64 a;
	while(scanf("%d", &n)!=EOF)
	{
		cnt = 0;
		while(n--)
		{
			scanf("%I64d", &a);
			if(miller_rabin(a))
			{
				cnt++;
				//printf("%d ", a);
			}
		}
		printf("%d\n",cnt);
	}
	return 0;
}

HDu 2138 How many prime numbers 高效Miller素数测试,布布扣,bubuko.com