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;
}

转载于:https://www.cnblogs.com/WArobot/p/7101805.html