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

Posted 2020-09-14 jzdwajue

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

求素数有測试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;
}