素数判定（素数）- HDU 2012
刚学编程的时候，我们大多需要做的一道题，那就是用C语言来判定一个数是否是素数。那时候很自然的会想到，对于数n，直接遍历一下n以下的数x，如果n%x等于0，说明可以被整除，也就不是素数。
这个是一个很正常的思维，因为不会理论数学分析的我们，很难想象得到我们可以用第二种方式解决这个问题，那就是开根号。
计算机专业书《离散数学及其应用》里面有一个委婉的定理，定理如下：
定理2：如果n是一个合数，那么必有一个素因子小于或等于sqrt(n)。
从定理2可知，如果一个整数不能被小于或等于其平方根的素数整除，则它就是素数 。
OK，我们的第二种解法就是遍历小于sqrt(n)的数。
但是，对于搞数学的来说，这个还是太慢了，因为他们研究的数字往往都是天文数字。
于是经过种种努力与机缘巧合，米勒·罗宾两个人研究出了一个测试算法，该算法也因此以他们的名字命名。
米勒·罗宾测试的错误率至多为1/2的s次方，s为迭代次数。数学家可能喜欢把s搞到20-50，程序员弄个4-10就可以了。
目前来说，这个算法是最快的！
这个算法可以看《算法导论》，里面讲得很详细，离散数学里面没有讨论这个算法，可见算法导论在追求性能的理论方面是做到了极致的。
算法思路如下。
另外一个想说的事情是，数学方面的题目由于理论性太强，就比如之前的hash函数，也是数学上面的结论，我们很难真正从本质是去理解为什么，也因此在源码中不做这方面的注释工作，数学方面的结论就当是模板函数，如果你天资聪颖，建议你还是花时间去看看。
好了，来看题。
Problem Description
对于表达式n^2+n+41，当n在（x,y）范围内取整数值时（包括x,y）(-39<=x<y<=50)，判定该表达式的值是否都为素数。
Input
输入数据有多组，每组占一行，由两个整数x，y组成，当x=0,y=0时，表示输入结束，该行不做处理。
Output
对于每个给定范围内的取值，如果表达式的值都为素数，则输出"OK",否则请输出“Sorry”,每组输出占一行。
Sample Input
0 1
0 0
Sample Output
OK
解题说明：对于任意的输入，测试一下即可。
源代码：G++ 0ms
#include<stdlib.h>
#include<stdio.h>
/***************米勒·罗宾***************/
#define MAXN 65
long long n, x[MAXN];
//multi函数网上还有一种实现，速度为O(1)
long long multi(long long a, long long b, long long p) {
long long ans = 0;
while (b) {
if (b & 1LL) ans = (ans + a) % p;
a = (a + a) % p;
b >>= 1;
}
return ans;
}
long long qpow(long long a, long long b, long long p) {
long long ans = 1;
while (b) {
if (b & 1LL) ans = multi(ans, a, p);
a = multi(a, a, p);
b >>= 1;
}
return ans;
}
int miller_rabin(long long n, int iterate) {
if (n == 2) return 1;
if (n % 2 == 0) return 0;
int s = iterate, i, t = 0;
long long u = n - 1;
while (!(u & 1)) {
t++;
u >>= 1;
}
while (s--) {
long long a = rand() % (n - 2) + 2;
x[0] = qpow(a, u, n);
for (i = 1; i <= t; i++) {
x[i] = multi(x[i - 1], x[i - 1], n);
if (x[i] == 1 && x[i - 1] != 1 && x[i - 1] != n - 1) return 0;
}
if (x[t] != 1) return 0;
}
return 1;
}
/***************米勒·罗宾·end***************/
int main()
{
int x, y;
//米勒·罗宾算法迭代次数
const int iterate = 4;
//输入x、y
while (scanf("%d%d", &x, &y) && x | y) {
//是否全部是素数
int is_all_prime = 1;
for (int i = x; i <= y; ++i) {
//进行米勒·罗宾素数测试
if (!miller_rabin(i * i + i + 41, iterate)) {
is_all_prime = 0;
break;
}
}
printf("%s\n", is_all_prime ? "OK" : "Sorry");
}
}