Project Euler Problem 027

http://projecteuler.net/index.php?section=problems&id=27
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%2027

Euler published the remarkable quadratic formula:

  n^2 + n + 41
  
It turns out that the formula will produce 40 primes for the consecutive values n = 0 to 39.
However, when n = 40, 402 + 40 + 41 = 40(40 + 1) + 41 is divisible by 41,
and certainly when n = 41, 41^2 + 41 + 41 is clearly divisible by 41.
Using computers, the incredible formula  n^2  79n + 1601 was discovered, which produces 80 primes for the consecutive values n = 0 to 79. The product of the coefficients, 79 and 1601, is 126479.
Considering quadratics of the form:

  n^2 + an + b, where |a| < 1000 and |b| < 1000

where |n| is the modulus/absolute value of n
e.g. |11| = 11 and |4| = 4
Find the product of the coefficients, a and b,
for the quadratic expression that produces the maximum number of primes for consecutive values of n,
starting with n = 0.

aとbの条件に当てはまる組み合わせを計算すると400万回弱計算しないといけないけど、あんまり計算したくないので、少し考える。

まず n = 0 のとき、

  f(0) = b

となるので、bは素数でないといけない。1000以下の素数は169個ある。

さらに n が偶数の場合、

  f(n) = 偶数 + 偶数 + b

となるので、b が 2 のときはf(偶数)が偶数になり、ダメっぽい。
したがって、 b は 3 以上 1000 以下の素数 168 個のいずれかになる。

つぎに a について考える。
n = 1 のとき、

  f(1) = 1 + a + b

となり、負数は素数でないので a は -b より大きい数である事が分かる。

さらに n が奇数のとき、

  f(n) = 奇数 + a * 奇数 + 奇数

となるので、a が奇数である事がわかる。(aが偶数だと 奇数 + 偶数 + 奇数 で f(n) = 偶数となる)

そこまで調査する数字を削減すれば十分かな。あとはグルグルと調べる。

import math
import itertools
import datetime


def primelist(limit=1000):
    checked = [ True ] * (limit/2)
    result = [2]
    for i in xrange(len(checked)):
        val = i * 2 + 3
        if checked[i]:
            result.append(val)
            j = val + val + val
            while j < limit:
                checked[(j - 3) / 2] = False
                j += val + val
    return result


cache = {}
def is_prime(n):
    if n < 0: n *= -1
    if n <  2: return False
    if n == 2: return True
    if not n & 0x01: return False
    if n in cache: return cache[n]
    
    x = 3
    limit = math.sqrt(n)
    while x <= limit:
        if n % x == 0:
            cache[n] = False
            return False
        x += 2
    cache[n] = True
    return True



def make_func(a, b):
    def f(n):
        return n * n + a * n + b
    return f



def count_consecutive_length(func, cond):
    for i in itertools.count():
        if not cond(func(i)): break
    return i


def euler027():
    maxcnt = 0
    result = 0
    
    it = iter(primelist(1000))
    it.next()
    
    for b in it:
        for a in xrange(-b, 1000, 2):
            func = make_func(a, b)
            cnt  = count_consecutive_length(func, is_prime)
            if cnt > maxcnt:
                maxcnt = cnt
                result = (a, b)
            
    print result, result[0] * result[1]


begin = datetime.datetime.now()
euler027()    
end   = datetime.datetime.now()
print end - begin

答え: -59231 (a=-61, b=971)
実行時間: 0.569957秒くらい