Project Euler 027を解いてみる。「二次式素数」

Project Euler 027

027

オイラーは以下の二次式を考案している:

n^2 + n + 41.

この式は, n を0から39までの連続する整数としたときに40個の素数を生成する.
しかし, n = 40 のとき 40² + 40 + 41 = 40(40 + 1) + 41 となり41で割り切れる.
また, n = 41 のときは 41² + 41 + 41 であり明らかに41で割り切れる.
計算機を用いて, 二次式

n^2 - 79n + 1601

という式が発見できた.

これは n = 0 から 79 の連続する整数で80個の素数を生成する. 係数の積は, -79 × 1601 で -126479である.
さて, |a| < 1000, |b| ≤ 1000 として以下の二次式を考える (ここで |a| は絶対値): 例えば |11| = 11 |-4| = 4である.

n^2 + an + b

n = 0 から始めて連続する整数で素数を生成したときに最長の長さとなる上の二次式の, 係数 a, b の積を答えよ.

->次の問題

考え方

今回も総当りで解いていきました。
a,bそれぞれの探索範囲ですが、
n=0の時、n²+an+b=bなので、bの探索範囲は1000未満の素数となります。
次に、n=1の時、1+a+bが素数、つまり最低限1+a+b≧2となるので、aの探索範囲は1-b以上1000未満とします。
毎回素数判定を行うと時間がかかるので、必要な素数の集合は先に作ってしまい、素数判定はその中に含まれるか否かで行います。
n=bの時、与式はbで割り切れるのでnはb未満、そしてbは1000未満の素数ということで997以下となります。
よって、生成されうる素数は最大でも(997² * 2 + 997)未満なのでここまでの素数のリストを使用します。
あとはforでループをぶん回します。

###　追記
b=2とするとn=2のときに必ず2で割れてしまうので、bは奇数
bが奇数かつnが奇数の時、aが偶数だと奇数+偶数+奇数=偶数で2で割れてしまうのでaは奇数となります。
aの探索範囲は奇数のみでいいのでもっと減らせますね。

コード

euler027.py

import numpy as np
from math import ceil


def main():
    limit = 1000
    # n = 0のとき素数となるのでbは素数
    b_arr = get_prime_arr(limit)
    # 生成されうる素数は b **2 + b **2 + b より小さいはず、先に素数の集合を作っておく
    prime_set = set(get_prime_arr(b_arr[-1] ** 2 * 2 + b_arr[-1]))
    # 答え用の変数
    max_a, max_b, max_n = 0, 0, 0
    for b in b_arr:
        # n = 1のとき、1+a+b>=2なので探索範囲はa>= -b+1となる奇数
        for a in range(-b, limit, 2):
            # n = bのとき、bで割り切れるのでnの範囲は0以上b未満
            for n in range(b):
                if n ** 2 + n * a + b not in prime_set:
                    break
            if n > max_n:
                max_n = n
                max_a = a
                max_b = b
    print(f'a={max_a},b={max_b}のとき、最長で{max_n}の素数を生成')
    print('Answer:', max_a * max_b)
    
    
def get_prime_arr(n: int) -> np.ndarray:
    """自然数n未満の素数を配列で返す関数
    Args:
        n: 返す素数の最大値を規定、最大値はn-1
    Returns:
        numpy.ndarray: n未満の素数の配列

    """
    boolean_arr = np.ones((n,), dtype=np.bool)  # n個の要素をもつのboolean配列(0~n-1に対応)を生成、全てTrueとする
    boolean_arr[0:2] = False  # [0,1]をFalseとする
    boolean_arr[4::2] = False  # 2以外の偶数をFalseとする
    boolean_arr[6::3] = False  # 3以外の3の倍数をFalseとする
    for i in range(5, ceil(n ** 0.5), 6):  # iは初項5、公差6の等差数列とする
        if boolean_arr[i]:  # 6k-1に対応。iに一致する要素がTrue=素数のとき
            boolean_arr[i * 2::i] = False  # iの倍数の要素(非素数)をFalseとする
        if boolean_arr[i + 2]:  # 6k+1に対応。
            boolean_arr[(i + 2) * 2::(i + 2)] = False
    return np.arange(n)[boolean_arr]  # boolean配列を整数の配列に変換して返す



if __name__ == '__main__':
    from time import time as t
    st = t()
    main()
    print(t() - st, 'sec')

paizaにて実行
結果 a=-61,b=971のとき、最長で71の素数を生成 answer= -59231 0.2921280860900879
はじめはprime_setのところをnumpy.ndarrayをそのまま使っており、45秒ほどかかってしまっていました。単に含まれるか否かの判定であればsetが格段に早いですね。

ちなみにlimit=50000までやってみましたが、
a=-79,b=1601のとき、最長で80の素数を生成がMaxでした。あまり大きい値では条件に合うa,bは見つかりませんでした(さらに大きくしたらどうなるかはわかりません)。