Project Euler - Problem 27

問題

しばらく止まってましたが今日から再開。

• 原文

  Considering quadratics of the form:

  n² + an + b, where |a| < 1000 and |b| < 1000

  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| < 1000, |b| < 1000 として以下の二次式を考える (ここで|a|は絶対値):

  n² + an + b

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

解答

最大探索範囲は-999 <= a <= 999、-999 <= b <= 999なので、およそ4,000,000通りの係数の組合せを試すことになります。組合せ毎に数列を生成して、それが素数か判定するわけですからたまりません。簡単な検討を加えて範囲を絞りましょう。

与えられた二次式をf(n)とおくと、f(0) = b、f(1) = a + b + 1です。 f(n)が長さ2以上の素数列を生成するならこれらは素数ですから、次のことがいえます：

1. bは素数である
2. a + b + 1は素数である
3. b = 2のとき、aは偶数である
4. それ以外のとき、aは奇数である

素数判定関数is_primeには同じ引数が与えられることがよくあるのでメモ化しています。

#!/usr/bin/perl

use strict;
use warnings;
use feature qw/say/;

sub prime_seq_len($$) {
  my ($coeff_a, $coeff_b) = @_;
  my $len = 0;
  my $n = 0;
  $len++, $n++ while is_prime($n * ($n + $coeff_a) + $coeff_b);
  return $len;
}

{
  my %primes = ( 2 => 1 );
  sub is_prime(_) {
    my $n = abs shift;

    unless (exists $primes{$n}) {
      return 0 if $n < 2 or $n % 2 == 0;

      my $is_prime = 1;
      for (my $i = 3; $i <= sqrt $n; $i += 2) {
        if ($n % $i == 0) {
          $is_prime = 0;
          last;
        }
      }
      $primes{$n} = $is_prime;
    }

    return $primes{$n};
  }
}

my $longest = 0;
my $product;
for my $coeff_b (grep { is_prime } -999 .. 999) {
  my $modulo = $coeff_b == 2 ? 0 : 1;
  my @a_cands = grep { $_ % 2 == $modulo
                         and is_prime($_ + $coeff_b + 1) } -999 .. 999;
  for my $coeff_a (@a_cands) {
    my $len = prime_seq_len($coeff_a, $coeff_b);
    if ($longest < $len) {
      $longest = $len;
      $product = $coeff_a * $coeff_b;
    }
  }
}

say $product;

追記

コメントで匿名氏に御指摘いただきました。

b = n = 2とするとf(2) = 2² + 2a + 2となり、これは明らかに（2より大きい）偶数です。 つまりb = 2のとき、生成する素数列の長さは高々2であり、無視できることになります。

前掲のコードでb = 2のケースを考慮している部分（$moduloのあたり）は不要になりました：

my $longest = 0;
my $product;
for my $coeff_b (grep { is_prime } -999 .. 999) {
  my @a_cands = grep { $_ % 2 != 0
                         and is_prime($_ + $coeff_b + 1) } -999 .. 999;
  for my $coeff_a (@a_cands) {
    my $len = prime_seq_len($coeff_a, $coeff_b);
    if ($longest < $len) {
      $longest = $len;
      $product = $coeff_a * $coeff_b;
    }
  }
}

say $product;