Project Euler 27 - 二次式素数 - もうカツ丼はいいよな

問題の概要

オイラーは次のような驚くべき関数を発見した。

$n^2 + n + 41$

この関数のnに、0から順に整数を入れて行くと、素数が順番に生成されていく。ただし、 $n=40$ のとき、結果は41で割り切れてしまう。つまりこの関数は0から39の範囲で素数を生成するのだ。また、 $n=41$ のとき、41で割り切れるというのは自明だろう。

さらに驚くべきことがある。次の式は80個の素数を連続して生成するのだ。

$n^2 - 79n + 1601$

一般化して、次のような関数を考えよう。

$n^2 + an + b$

ここで、 $|a| \lt 1000$ かつ $|b| \leq 1000$ とする。最も多くの素数を連続して生成する $a$ と $b$ の組み合わせは何か。 $a$ と $b$ の積として答えよ。

回答例

そのままでは、aとbの組み合わせで200万、それに素数が連続するかどうかの判定が入るので、計算量はかなり多くなる。

しかし、条件から探索範囲について次のことがわかる。

bは素数である（nが0から始まるため）
aは-bより大きい（素数であるためには式の結果が正でなければならない）
nの上限はb-1である（n=bで必ず割り切れるため）
式の値の上限はほぼ200万である（ $1000^ 2 + 1000 \times 1000 + 1000$ 。bは素数なのでもう少し小さい。)

まず、4より、200万以下の素数のリストがあれば充分ということが分かる。この程度ならエラトステネスの篩で高速に生成できる。

limit <- 2e6
sieve <- logical(limit) 
sieve[1] <- TRUE # prime = FALSE
sieve[seq(4, limit, by = 2)] <- TRUE
for (i in seq(3, sqrt(limit), by = 2)) {
  if(!sieve[i]) sieve[seq(i*2, limit, by = i)] <- TRUE
}
primes <- (1:limit)[!sieve]

あとは、無駄な範囲を避けながらループを回せば良い。

chk_plen <- function(a, b) {
  cnt <- 0
  for (n in 1:(b-1)) {
    if(n^2 + a*n + b < 1) break
    if(sieve[n^2 + a*n + b]) break
    cnt <- cnt + 1
  }
  return(cnt)
}

max_plen <- 0
max_a <- 0
max_b <- 0
for(b in primes[primes < 1000]) {
  for(a in (-b+1):999) {
    plen <- chk_plen(a, b)
    if (plen > max_plen){
      max_plen <- plen
      max_a <- a
      max_b <- b
    }
  }
}
max_a * max_b