ScalaでProject Euler（68）

Scala

Problem 39

この問題はピタゴラス数の生成を使えば簡単です。

a = 2lmn b = l(m2 - n2) c = l(m2 + n2)
p = 2lm(m + n)

ここで、n' = m + nとおけば、

p / 2 = lmn'
(m, n') = 1 m < n' < 2m

p / 2を素因数分解して3つの積に分解してそれが上の条件を満たす組合せを数えればよいです。
もう少し工夫します。3つの積に分解するとき、まずlを決めて残りをmとn'に分配するのですが、2の因子はmにのみ分配されて、他の因子は同じ因子が全部同じ方に分配されます。すなわち、mn' = 23・32・53なら、

(m, n') = (23・32・53, 1), (23・32, 53), (23・53, 32), (23, 1, 32・53)

の4通りのみです。しかし、このロジックを入れてもN = 100000で2倍の速度にしかなりませんでした。

import scala.testing.Benchmark

def pow(n :Int, e :Int) = (1 to e).foldLeft(1)((x, y) => x * n)

type Facts = List[(Int,Int)]

def value(fs :Facts) =
    fs.foldLeft(1)((x, y) => x * pow(y._1, y._2))

def divide(fs1 :Facts, fs2 :Facts) :Facts = (fs1, fs2) match {
    case (_, Nil) => fs1
    case ((p1, e1) :: t1, (p2, e2) :: t2) if p1 == p2 =>
                if(e1 == e2) divide(t1, t2) else (p1, e1 - e2) :: divide(t1, t2)
    case (h1 :: t1, _) => h1 :: divide(t1, fs2)
    case _ => Nil
}

def sieve(N :Int) = {
    val a = Array.range(0, N)
    for(p <- Iterator.from(2).takeWhile(n => n * n < N).filter(n => a(n) == n);
                    k <- Iterator.range(p * 2, N, p))
        a(k) = p
    a
}

def divs(fs :Facts) :Iterator[Facts] = fs match {
    case Nil => Iterator(Nil)
    case (p, e) :: fst => for(fs2 <- divs(fst); e1 <- Iterator.range(0, e + 1))
                            yield if(e1 == 0) fs2 else (p, e1) :: fs2
}

def divs2(fs :Facts) :Iterator[Facts] = fs match {
    case Nil => Iterator(Nil)
    case (2, e) :: fst => divs2(fst).map((2, e) :: _)
    case f :: fst => divs2(fst).flatMap(x => Iterator(x, f :: x))
}

def divs3(fs :Facts) :Iterator[(Int,Int,Int)] =
    for(fs1 <- divs(fs); fs2 = divide(fs, fs1); fs3 <- divs2(fs2))
        yield (value(fs1), value(fs3), value(divide(fs2, fs3)))

def div_pow(n :Int, d :Int) :(Int,Int) =
    if(n % d != 0)
        (0, n)
    else {
        val (e, m) = div_pow(n / d, d)
        (e + 1, m)
    }

def is_right_triangle(x :(Int,Int,Int)) = {
    val (l, m, n) = x
    m < n && n < m * 2 && n % 2 == 1
}

def solve() = {
    val a = sieve(N / 2 + 1)
    
    def factorize(n :Int) :Facts =
        if(n == 1)
            Nil
        else {
            val p = a(n)
            val (e, m) = div_pow(n, p)
            (p, e) :: factorize(m)
        }
    
    def max(x :(Int,Int), y :(Int, Int)) = (x, y) match {
        case ((n1, m1), (n2, m2)) if m1 > m2 => (n1, m1)
        case ((n1, m1), (n2, m2))            => (n2, m2)
    }
    
    println (Iterator.range(2, N + 1, 2).map(n => (n, divs3(factorize(n / 2)).
                            filter(is_right_triangle).size)).reduceLeft(max))
}

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

val N = 100000
println (Test.runBenchmark(5))