ScalaでProject Euler（51） - inamori’s diary

Problem 29

もう一度考え直してみましょう。N = 100で底8の重複を考えてみましょう。底2に重複するのは公差1で33まで、底4は公差2で66まで。とりあえず8¹も含めて両方とも33個あります。ただし、ダブっている部分があって、公差2で33までだから16個あるので33 + 33 - 16 = 50で、8¹の分を引いて49個です。包除原理ですね。
包除原理は、次のように考えるとわかりやすいかもしれません。全体集合Uと部分集合S₁, S₂, ... , S_nとして、

U - S₁∪S₂∪...∪S_n = ( U - S₁ )×( U - S₂ )×...×( U - S_n )

ここで、×は

|A|×|B| = |A∩B|

です。明らかに記号として間違っていますが、あくまでイメージです。
ここで、S₁を先ほどの底2に重複する集合、S₂を底4に重複する集合として、

S₁ = (1, 33, 1)
S₂ = (2, 66, 1)

と書きます。(公差, 上限, 交わる集合の個数)です。そうすると、

S₁×S₂ = (2, 33, 2)

となります。一般には、

(d₁, u₁, n₁)×(d₂, u₂, n₂) = (lcm(d₁, d₂), min(u₁, u₂), n₁ + n₂)

という簡単な計算になります。符号は交わる集合の個数をたよりにして底8を計算すると、

100 - 33 / 1 - 66 / 2 + 33 / 2 - 1 = 49

となります。

min

小さいほうを選ぶには、

min(1, 2)

ではなく、

1.min(2)

なんですね。オブジェクト指向は対称性の犠牲の上に成り立っています。

さて、以上の考えのもとに下のように組んでみました。たぶん、(d, u, n)で、d > uとなったら無視するというのが高速化のポイントです。これで、メモリの制限で7×10⁸が限界でした。4秒弱です。Iteratorを使っているのになぜメモリを食うのでしょうね。

（追記）
少しコードを間違えていて、修正したコードで5×10⁸までで1.3秒程度。6×10⁸にすると30秒くらいかかっていました。原因不明です。

（追記）
単に素数乗のとき遅くて、29乗のときの計算をしなければならないと遅い、ということのようです。

import scala.math
import scala.testing.Benchmark

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

def lcm(n :Int, m :Int) = n / gcd(n, m) * m

def pow(n :Int, e :Int) :Int =
    if(e == 0)
        1
    else {
        val m = pow(n, e / 2)
        if(e % 2 == 1)
            m * m * n
        else
            m * m
    }

// [n^(1/e)]
def int_root(n :Int, e :Int) = {
    def dec(m :Int) :Int =
        if(pow(m - 1, e) < n) m else dec(m - 1)
    
    def inc(m :Int) :Int =
        if(pow(m + 1, e) > n) m else inc(m + 1)
    
    val m = math.pow(n, 1. / e).toInt
    if(pow(m, e) > n) dec(m) else inc(m)
}

def divs(n :Int) =
    Iterator.range(1, n).filter(n % _ == 0)

def num_pows() = {
    val a = Iterator.from(2).map(int_root(N, _) - 1).takeWhile(0 <).toArray
    for(e <- Iterator.range(a.size + 1, 3, -1); e2 <- divs(e).drop(1))
        a(e2 - 2) -= a(e - 2)
    Iterator.range(0, a.size).map(a(_))
}

def fraction(num :Int, den :Int) = {
    val d = gcd(num, den)
    (num / d, den / d)
}

def count_overlap_old(e :Int) = {
    val a = new Array[Boolean](N + 1)
    val fs = List.range(1, e).map(fraction(_, e)).toMap.toList
    for((n, d) <- fs; k <- Iterator.range(n, N / d * n + 1, n))
        a(k) = true
    Iterator.range(2, N + 1).count(k => a(k))
}

def count_overlap(e :Int) :Int = {
    def mul(x :(Int,Int,Int), y :(Int,Int,Int)) = {
        val ((d1, u1, n1), (d2, u2, n2)) = (x, y)
        (lcm(d1, d2), u1.min(u2), n1 + n2)
    }
    
    def mul_iter(it_x :Iterator[(Int,Int,Int)], y :(Int,Int)) =
        it_x.flatMap(x => Iterator((1, N, 0), (y._1, N / y._2 * y._1, 1)).
                                map(mul(x, _)).filter(x => x._1 <= x._2))
    
    val fs = List.range(1, e).map(fraction(_, e)).toMap.toIterator
    val it_n = fs.foldLeft(Iterator((1, N, 0)))(mul_iter)
    N - it_n.map(x => if(x._3 % 2 == 0) x._2 / x._1 else -x._2 / x._1).sum - 1
}

val N = 5e8.toInt
def solve() = {
    val a = num_pows()
    val it = a.zip(Iterator.from(2).map(count_overlap))
    println ((N.toLong - 1) * (N - 1) - it.map(x => x._1.toLong * x._2).sum)
}

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

println (Test.runBenchmark(10))