問題146

問題 146 です。数学的解法により高速化するというより、プログラムの高速化技法を用いて解く問題だと思いました。
(いつも使用している Scheme ではなくて、高速で動作する c で書いたのはヒミツです)

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdbool.h>
#include <math.h>

bool is_prime(int* prime_array, int prime_cnt, long n)
{
  long p = (long)sqrt(n);

  for (int i = 0; i < prime_cnt && p >= prime_array[i]; i++) {
    if (n % prime_array[i] == 0) {
      return false;
    }
  }

  return true;
}

bool is_consecutive_primes(int* prime_array, int prime_cnt, long nn)
{
  long p = sqrt(nn + 27);

  for (int i = 0; i < prime_cnt && p >= prime_array[i]; i++) {
    if ((nn + 1) % prime_array[i] == 0 ||
      (nn + 3) % prime_array[i] == 0 ||
      (nn + 7) % prime_array[i] == 0 ||
      (nn + 9) % prime_array[i] == 0 ||
      (nn + 13) % prime_array[i] == 0 ||
      (nn + 27) % prime_array[i] == 0) {
      return false;
    }
  }

  return true;
}

int main(int argc, char* argv[])
{
  const int limit = 150000000L;
  //const long limit = 1000000L;
  int prime_cnt = limit;

  bool* primes = (bool*)calloc(limit, sizeof(bool));
  memset(primes, true, limit);

  for (int i = 2L; i < limit; i++) {
    if (primes[i]) {
      for (int j = i + i; j < limit; j += i) {
        if (primes[j]) {
          prime_cnt--;
          primes[j] = false;
        }
      }
    }
  }

  prime_cnt -= 2;
  int* prime_array = (int*)calloc(prime_cnt, sizeof(int));
  for (int i = 2, j = 0; i < limit; i++) {
    if (primes[i]) {
      prime_array[j] = i;
      j++;
    }
  }
  free(primes);

  long ans = 0L;
  for (long n = 2L; n < limit; n++) {
    long nn = (long)n * (long)n;
    if (nn % 2L == 0 && is_consecutive_primes(prime_array, prime_cnt, nn)) {
      if (!is_prime(prime_array, prime_cnt, nn + 5) &&
        !is_prime(prime_array, prime_cnt, nn + 11) &&
        !is_prime(prime_array, prime_cnt, nn + 15) &&
        !is_prime(prime_array, prime_cnt, nn + 17) &&
        !is_prime(prime_array, prime_cnt, nn + 19) &&
        !is_prime(prime_array, prime_cnt, nn + 21) &&
        !is_prime(prime_array, prime_cnt, nn + 23) &&
        !is_prime(prime_array, prime_cnt, nn + 25)) {
        ans = ans + n;
      }
    }
  }

  printf("%ld\n", ans);
  free(prime_array);
  return 0;
}

速度を早くするために工夫した点は、
1)素数の生成方法
2)連続する素数の判定
の部分です。
それぞれ説明していきます。

素数の生成方法

純粋に素数を探し出してしまうと処理が遅くなるため素数のテーブルを生成します。
まずエラトステネスの篩を使用して素数か非素数かの配列 (プログラム中では primes) を作り、この配列から素数配列(primes_array)を生成しています。メモリを浪費することで処理速度を上げる方法です。とはいえ、n の素数判定を行うには sqrt(n) まで調べれば良いので、素数配列の要素数についてもコンパクトになるように工夫しています。

calloc や free なんて今どき使用していて大丈夫なのかという気もしますがw、素数配列を求めるコードは他の素数が関係する問題を解く時にも流用できると思います。

連続する素数の判定

連続する素数の判定は、is_consecutive_primes() で行っていますが、調べなければいけない素数 n^2 + k (k=1,3,7,9,13,27) を毎回それぞれ素数配列を使用して検索するのは処理時間が大幅に無駄になるため、連続する素数の最大値 n^2 + 27 まで素数配列を検索する際に、すべての素数判定 n^2 + k (k=1,3,7,9,13,27) をまとめて行うようにしています。これにより素数配列を検索する回数が 1/6 になります。

また、n^2 + k (k=5,11,15,17,19,21,23,25) については素数でないことを確認しなければなりませんが。これはそれぞれ別個に素数配列を検索しています。この処理は候補が見つかった場合にしか実行しないため処理時間は大して増えませんから、もうちょっとどうにかしたかったですがこれで良しとしました。

コンパイラの違いによる注意点

また、gccコンパイルしているため long は 8 バイト長となり正常動作しますが、Visual Studio など他のコンパイラを使用する場合は long は 4 バイト長になったりするため、long long (8バイト長) に修正しなければ正常動作しない事があります。注意してください。
(コンパイラによって正常動作したりしなかったりというコードはどうなんだ と思います。はいw)

課題
ポインタ持ちにすると配列のサイズが求められなくなるため、素数配列生成部分はモジュール分割をしておらず、ひどいコードになっています。

実行結果

$ gcc -O2 -o pe146 pe146.c

Emacs からコンパイルする場合は C-c c

$ time ./pe146.exe 
676333270

real	0m15.796s
user	0m15.500s
sys	0m0.078s

最初に書いたプログラムは数時間経っても解答できなかったため、15 秒強かかっているものの大幅に処理速度は改善したと思います。
(最初に書いたプログラムが駄目すぎるという説もありますw)