久しぶりにProject Euler。
問題の概要
オイラーは次のような驚くべき関数を発見した。
この関数のnに、0から順に整数を入れて行くと、素数が順番に生成されていく。ただし、のとき、結果は41で割り切れてしまう。つまりこの関数は0から39の範囲で素数を生成するのだ。また、 のとき、41で割り切れるというのは自明だろう。
さらに驚くべきことがある。次の式は80個の素数を連続して生成するのだ。
一般化して、次のような関数を考えよう。
ここで、かつとする。最も多くの素数を連続して生成するとの組み合わせは何か。との積として答えよ。
回答例
そのままでは、aとbの組み合わせで200万、それに素数が連続するかどうかの判定が入るので、計算量はかなり多くなる。
しかし、条件から探索範囲について次のことがわかる。
- bは素数である(nが0から始まるため)
- aは-bより大きい(素数であるためには式の結果が正でなければならない)
- nの上限はb-1である(n=bで必ず割り切れるため)
- 式の値の上限はほぼ200万である(。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