ScalaでProject Euler（98）

Problem 64

ループの長さが奇数になる数を具体的に並べてみましょう。

2, 5, 10, 13, 17, 26, 29, 37, 41, 50, 53, 58, 61, 65, 73, ...

これらを素因数分解すると、因数は2か4で割って1余る素数になっています。
ただ、逆は違うようです。すなわち、すべて2か4で割って1余る素数から成っていてもループの長さが奇数になるとは限りません。すぐにわかるのは、2は一つしか使ってはいけないということです。例えば、20などはこの中に入りません。それ以外では、例えば、34は偶数のようです。これらの例外を並べてみると、

34, 146, 178, 194, 205, 221, 305, 377, 386, 410, ...

素数は例外にならないようです。それ以外は例外に成り得ます。なので、素数以外の場合は前回行った判定をしなければなりません。このロジックで問題無いかは1000万まで確認しました。

これで組んでみたところ、前回のコードより10倍以上速くなりました。

import scala.math.sqrt
import scala.testing.Benchmark

def gcd(m :Int, n :Int) :Int = if(n == 0) m else gcd(n, m % n)

def num_steps(N :Int) :Int = {
    val n = sqrt(N.toDouble).toInt
    val x0 = (1, n, 1, 0)
    
    def next(x :(Int,Int,Int,Int)) :(Int,Int,Int,Int) = {
        val (a, b, c, m) = x
        val m1 = (a * n + b) / c
        val e = m1 * c - b
        val (a1, b1, c1) = (a * c, c * e, a * a * N - e * e)
        val d = gcd(gcd(a1, b1), c1)
        (a1 / d, b1 / d, c1 / d, m1)
    }
    
    def equalto(x :(Int,Int,Int,Int)) =
        x._1 == x0._1 && x._2 == x0._2 && x._3 == x0._3
    
    Iterator.iterate(next(x0))(next).takeWhile(x => !equalto(x)).size + 1
}

// 250, 5 -> (3, 2)
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)
    }

// [ 0, 1, 2, 3, 2, 5, 2, 7, 2, 3, ... ]
def sieve(M :Int) :Array[Int] = {
    val a = new Array[Int](M + 1)
    for(k <- 2 to M) a(k) = k
    for(p <- Iterator.from(2).takeWhile(n => n * n <= M) if a(p) == p;
        k <- Iterator.range(p * 2, M + 1, p) if a(k) == k)
            a(k) = p
    a
}

def is_odd_loop_length(n :Int, a :Array[Int]) = {
    def f(m :Int, bsq :Boolean) :Boolean = {
        val p = a(m)
        if(m == 1)
            !bsq
        else if(p % 4 == 3)
            false
        else if(p == m)
            true
        else {
            val (e, q) = div_pow(m, p)
            if(p == 2 && e > 1)
                false
            else
                f(q, bsq && (e % 2 == 0))
        }
    }
    
    if(a(n) == n)       // prime
        (n & 3) != 3
    else
        f(n, true) && num_steps(n) % 2 == 1
}

def solve() = {
    val M = 1000000
    val a = sieve(M)
    println ((2 to M).filter(n => is_odd_loop_length(n, a)).size)
}

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

println (Test.runBenchmark(10))