ScalaでProject Euler（93）

Scala

Problem 60

まず、素数列を用意します。これにどの向きに結合しても素数になるような素数を探します。例えば(3, 7)ですね。これを結合素数と呼ぶことにします。こうして2つの素数の組を列挙します。これにさらにもう一つ素数を加えて3つの素数による結合素数を列挙します。こうして5個の素数による結合素数を探します。実際には、素数の通し番号の組と和のタプルを持つことにします。(3, 7)なら((1, 3), 10)です。なぜなら、(3, 7)に結合する素数は7より大きいものだけ考えればよいからです。

さて、列挙すると書きましたが、上限がわからないので難しいですよね。小さい方から列挙することができればいいのですがそれも難しいです。なぜなら、例えば(43, 97)と(7, 127)は総和は後者の方が小さいですが、素数をさらに一つ加えたときは前者の方が小さくなる可能性があるからです。

そこで上限を適当に設定します。そしてその上限で計算します。解が無ければ上限を2倍します。これを解が現れるまで繰り返します。13秒でした。まあまあですね。
ところで、これを実際に計算してみると上限2倍で6、7倍の時間がかかります。6倍とすると、最悪の場合は最良の場合に比べて6倍の時間がかかることになります。こういうときは、21/log26（≒1.3）倍ずつしていくとよいです。こうすると最悪の場合でもそんなに悪くならないです（あとで詳しく書くと思います）

import scala.collection.mutable.ArrayBuffer
import scala.testing.Benchmark

object Prime {
    val primes = new ArrayBuffer[Int]
    primes += 2
    
    def apply(k :Int) =
        if(k < primes.size)
            primes(k)
        else {
            for(n <- Stream.from(primes(primes.size-1) + 1).
                    takeWhile(n => k >= primes.size) if is_prime(n))
                primes += n
            primes(primes.size-1)
        }
    
    def is_prime(n :Int, k :Int = 0) :Boolean = {
        val p = Prime(k)
        if(p * p > n)       true
        else if(n % p == 0) false
        else                is_prime(n, k + 1)
    }
}

def cat(n :Int, m :Int) :Int = {
    def pow10(d :Int) :Int =
        if(d > m) d else pow10(d * 10)
    
    n * pow10(1) + m
}

def is_cat_prime(n :Int, m :Int) :Boolean =
    Prime.is_prime(cat(n, m)) && Prime.is_prime(cat(m, n))

def is_all_cat_prime(ps :List[Int], q :Int) :Boolean = ps match {
    case Nil => true
    case p :: tail if is_cat_prime(p, q) => is_all_cat_prime(tail, q)
    case p :: tail => false
}

def primes(k :Int, limit :Int) :List[Int] = {
    val p = Prime(k)
    if(p >= limit) Nil else p :: primes(k + 1, limit)
}

def make_sets(N :Int, num :Int, a :List[(Int,List[Int])]) =
    if(num == 1)
        primes(0, N / M).zipWithIndex.map(x => (x._1, List(x._2)))
    else
        for((sum_p, ks) <- a;
            limit = (N - sum_p) / (M - ks.size);
            (p, k) <- primes(ks.head + 1, limit).zipWithIndex;
            ps = ks.map(k => Prime(k))
            if(is_all_cat_prime(ps, p)))
                yield (sum_p + p, (k + ks.head + 1) :: ks)

def solutions(N :Int) :List[(Int,List[Int])] = {
    val a = new Array[List[(Int,List[Int])]](M + 1)
    for(n <- 1 to M)
        a(n) = make_sets(N, n, a(n-1))
    a(M)
}

def solve() = {
    val it = Iterator.iterate(1000)(2 *).map(solutions).dropWhile(Nil ==)
    println (it.next.map(_._1).min)
}

object Test extends Benchmark {
    def run() = solve
}

val M = 5
println (Test.runBenchmark(1))