Project Euler (74 〜 80) ピタゴラス数、五角数定理、トポロジカルソートなど - ブートストラッピングでコンパイラを作る日記

今日は80まで。このところ、数学を使う問題が沢山で面白いですね。

ピタゴラス数とブロコット木 Problem75

簡単には

1,500,000以下のLで
a^2 + b^2 = c^2 かつ a + b + c = L
を満たす自然数a,b,cの組がちょうど一組となるようなLの個数を求めよ

という問題です。a^2 + b^2 = c^2を満たす自然数(a,b,c)をピタゴラス数と言います。また、この方程式をピタゴラス方程式と言いますが、ペル方程式などと同じくディオファントス方程式の一種です。
(a,b,c)がピタゴラス数の時kを自然数として(ka,kb,kc)もピタゴラス数なので、(a,b,c)が互いに素であるものについて考えれば良いです。これを原始ピタゴラス数といいます。
さて、原始ピタゴラス数はその生成方法が知られていて

m > n かつ m,nは互いに素かつ片方が偶数のとき(m^2-n^2,2mn,m^2+n^2)は原始ピタゴラス数

となります。
これを列挙する時に厄介なのが「m,nが互いに素」の条件ですが、そこでProblem 71, 73でも使ったブロコット木が役に立ちます。
ブロコット木を使うと既約分数を高速に列挙する事ができるわけですが、このとき分子と分母が必ず互いに素になっている性質が使えます。

L = 1500000
$count = {} # 解の個数
def dfs(la,lb,ra,rb)
    # naとnbは互いに素
    na = la + ra
    nb = lb + rb
    c = 0
    if (na+nb)%2 == 1
        l = 2*(na+nb)*nb
        return if l > L
        # 原始ピタゴラス数の倍数をカウント
        (1..L/l).each do |k|
            $count[k*l] || $count[k*l] = 0
            $count[k*l] += 1
        end
    end
    dfs(la,lb,na,nb)
    dfs(na,nb,ra,rb)
end
dfs(0,1,1,1)
count = 0
$count.each_value {|v| count += 1 if v == 1}
puts count

実行時間は1.11sでした。

整数分割、五角数定理(Problem 76, 77, 78)

この3問は整数分割に関する問題です。整数分割とは整数nを自然数の和で表す場合の数です。例えばn=5の場合は以下の7通りがあります。

5
4 + 1
3 + 2
3 + 1 + 1
2 + 2 + 1
2 + 1 + 1 + 1
1 + 1 + 1 + 1 + 1

とどのつまりこれは、硬貨支払い問題と同じです。5円を1円玉、2円玉、3円玉、4円玉、5円玉の組み合わせで支払う方法は7通りという訳です。従って以前書いた方法と全く同じアルゴリズムで計算できます。

Problem 76: 100を2つ以上の自然数の和で表す方法の数を求めよ

100円を1,2,...,99円玉で表す方法が答えです。よって硬貨支払い問題と同じ漸化式を使うと以下のようになります。a[k][n]はnをk以下の和で表す方法の数です。

class Array
    def [](i)
        i < 0 ? 0 : at(i)
    end
end

N = 100
a = (0..N).collect { [] }
(0..N).each do |n|
    t = a[1][n] = a[1][n-1] + ((n == 0) ? 1 : 0)
    (2..N).each do |k|
        t = a[k][n] = a[k][n-k] + t
    end
end
puts a[99][100]

実行時間は0.01sくらいです。

Problem 77: 素数のみの和で表す方法が5000通り以上となる最初の数は何か？

2,3,5,7,11,13,..円硬貨で支払う方法の数と考えれば全く同じです。金額の上限は適当に決める必要があります。

class Array
    def [](i)
        i < 0 ? 0 : at(i)
    end
end

#上限を適当に決める
N = 100

# 素数のふるいをつくっておく
sieve = []
2.upto(N) do |p|
    next if sieve[p]
    (2*p).step(N,p) {|k| sieve[k] = true}
end

a = (0..N).collect { [] }
(0..N).each do |n|
    t = a[2][n] = a[2][n-2] + ((n == 0) ? 1 : 0)
    (3..N).each do |k|
        next if sieve[k] # kが合成数の時をスキップ
        t = a[k][n] = a[k][n-k] + t
        if a[k][n] > 5000
            puts n
            exit
        end
    end
end

実行時間は0.01sです。

この考え方は分かりやすくて良いのですが、メモリを多く使用するという欠点があります。そこで役に立つのがオイラーの五角数定理です。これは
$\prod_{n=1}^{\infty}(1-x^n) = \sum_{n=-\infty}^{\infty}(-1)^nx^{n(3n-1)/2}$
というものです。
nの整数分割の数p(n)の母関数は
$\sum_{n=0}^{\infty}p(n)x^n = \prod_{n=1}^{\infty}(1 + x^n + x^{2n} + \cdots) \\ = \prod_{n=1}^{\infty}\frac{1}{1-x^n}$
となるので、五角数定理を使えば
$\sum_{n=0}^{\infty}p(n)x^n = \left( \sum_{n=-\infty}^{\infty}(-1)^nx^{n(3n-1)/2} \right)^{-1}$
となります。これを展開して係数比較をすれば、分割数の漸化式が作れます。
Problem 78はこの方法で解けました。まだ計算時間が掛かっているのでもっと良い方法があるかもしれません。

Problem 78: 分割数が初めて1,000,000の倍数となる整数を求めよ

p = [1]
1.upto(1/0.0) do |n|
    p[n] = 0
    1.upto(1/0.0) do |k|
        if (m = n - k*(3*k+1)/2) >= 0
            p[n] += p[m] * (2*(k%2)-1)
        end
        if (m = n - k*(3*k-1)/2) >= 0
            p[n] += p[m] * (2*(k%2)-1)
        else
            break
        end
    end
    if p[n] % 1000_000 == 0
        puts n
        exit
    end
end

実行時間は41.5sでした。

トポロジカルソート(Problem 79)

ある長さ不明の数列がある。その長さ3の部分数列が
319
680
180
690
129
...(略)
となるとき、もとの数列として可能な物で最も短い物を求めよ

これは仮にもとの数列に同じ数が複数回出現しないなら、トポロジカルソートの問題です。
そこで、まずGraphvizでグラフを描いてみました。DOT言語で記述します。

digraph problem79 {
    3 -> 1 -> 9
    6 -> 8 -> 0
    1 -> 8 -> 0
    6 -> 9 -> 0
    1 -> 2 -> 9
...
}

この絵から答えは一目瞭然でした。
一応Rubyでトポロジカルソートを行うコードも載せておきます。もし、グラフにループがあれば無限ループになってしまう非常に雑なコードですが..。

file = open("problem79-keylog.txt")
# グラフを作る
$graph = (0..9).collect {Array.new(10,false)}
$used = []
while line = file.gets
    x, y, z = line.split(//).map {|v| v.to_i}
    $used[x] = $used[y] = $used[z] = true
    $graph[x][y] = $graph[y][z] = true
end

# トポロジカルソート
$visited = []
def visit(v)
    return if $visited[v]
    $visited[v] = true
    0.upto(9).each do |s|
        next if not $used[s] or not $graph[s][v]
        visit(s)
    end
    puts v
end
0.upto(9).each {|v| $used[v] && visit(v)}

実行時間は0.01s。

整数探索に直す(Problem 80)

平方数でない100未満のnについて
√nの最初100桁の数の和
の和を求めよ

開平法という手もありますが実装が面倒です。今回は100桁と分かっているので、もっと簡単な整数の探索問題に直しちゃう事ができます。つまりnが100未満ならば
$k \le 10^{99} \sqrt{n} < k + 1 \leftrightarrow k^2 \le n10^{198} < (k+1)^2$
を満たす整数kがその最初の100桁分です。ということで、以下の様に解いてみました。

sum = 0
2.upto(99) do |a|
    next if Math.sqrt(a).ceil**2 == a
    # 二分探索
    n = a*10**198
    l = 10**99
    r = (a+1)*10**99
    while r - l > 1
        t = (r+l)/2
        t**2 - n > 0 ? r = t : l = t
    end
    # 桁の和を取る
    while l > 0
        sum += l%10
        l /= 10
    end
end
puts sum

実行時間は0.26sでした。