Project Euler 27 Quadratic primes（ 米勒测试 + 推导性质 ）

---

题意：

欧拉发现了这个著名的二次多项式：

f(n) = n² + n + 41

对于连续的整数n从0到39，这个二次多项式生成了40个素数。然而，当n = 40时40² + 40 + 41 = 40(40 + 1) + 41能够被41整除，同时显然当n = 41时，41² + 41 + 41也能被41整除。

随后，另一个神奇的多项式n² − 79n + 1601被发现了，对于连续的整数n从0到79，它生成了80个素数。这个多项式的系数-79和1601的乘积为-126479。

考虑以下形式的二次多项式：

• n² + an + b, 满足|a| < 1000且|b| < 1000

• 其中|n|指n的模或绝对值
  例如|11| = 11以及|−4| = 4

这其中存在某个二次多项式能够对从0开始尽可能多的连续整数n都生成素数，求其系数a和b的乘积。

思路：根据上面的式子 f(n) 可以进行推导几个性质：

1. 当 n = 0 时 f(0) = b ，因为f(n)是素数所以 b 一定是素数
2. 当 n = 1 时 f(1) = a + b + 1，因为f(n) > 0 所以 a > -b - 1
3. 当 gcd(a , b) = d 且 d != 1时，一定会枚举到一个数 n' = minP(d) (代表d的最小素因数) 使得 f(n') 为合数，即当 a , b 不互质时，最多生成 minP(d) 个连续的素数。

根据这三条性质就可以优化代码，降低程序运行时间。

这个题中米勒测试没有增加二次探测，所以判断素数还是有几率出错的，虽然这个几率小的惊人，但万一非洲人呢 : )

---

/*************************************************************************
    > File Name: euler027.c
    > Author:    WArobot 
    > Blog:      http://www.cnblogs.com/WArobot/ 
    > Created Time: 2017年06月30日 星期五 20时42分40秒
 ************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <inttypes.h>

#define MAX_ROUND 30
#define MAX_N 1000

int32_t prime[MAX_N + 10] = {0};
int32_t primeList[MAX_N + 10] = {0};

void Init() {
	prime[1] = 1;
	for (int32_t i = 2 ; i <= MAX_N ; i++) {
		if (!prime[i]) { 
			prime[i] = i;
			primeList[ ++primeList[0] ] = i; 
		}
		for (int32_t j = 1 ; j <= primeList[0] ; j++) {
			if (i * primeList[j] > MAX_N)	break;
			prime[i * primeList[j]] = primeList[j];
			if (i % primeList[j] == 0)		break;
		}
	}
}

int32_t quick_power(int32_t a , int32_t b , int32_t mod) {
	int32_t ret = 1;
	while(b) {
		if (b & 1) ret = ret * a % mod;
		a = a * a % mod;
		b >>= 1;
	}
	return ret % mod;
}
bool R_M_TEST(int32_t x) {
	if (x <= 1)		return false;
	srand(time(NULL));
	for (int32_t i = 0 ; i < MAX_ROUND ; i++) {
		int32_t a = (rand()) % (x - 1) + 1;
		if (quick_power(a , x - 1 , x) != 1)	return false;
	}
	return true;
}

int32_t HowManyPrime(int32_t a , int32_t b) {
	int32_t i = 0;
	while (R_M_TEST(i * i + a * i + b)) i++;
	return i;
}

int32_t gcd(int32_t a , int32_t b) {
	return b == 0 ? a : gcd(b , a % b);
}

int32_t main() {
	int32_t a , b , d , maxQuadraticPrime = 0 , ans , manyPrime;
	Init();
	for (int32_t i = 1 ; i <= primeList[0] ; i++) {
		b = primeList[i];
		for (a = -b - 1 ; a < MAX_N ; a++) {
			if (prime[a + b + 1] != a + b + 1)	continue;
			d = gcd(a , b);
			manyPrime = HowManyPrime(a , b);
			if (d != 1 && maxQuadraticPrime >= prime[d])	continue;
			if (maxQuadraticPrime < manyPrime) {
				maxQuadraticPrime = manyPrime;
				ans = a * b;
			}
		}
	}
	printf("ans = %d\n",ans);
	return 0;
}