Project Euler(51 〜 73) ブロコット木など

数学

久々にProject Eulerやりました。まだrowlは使える段階にないので、Rubyつかってます。
このあたりから数論の問題がいろいろ出てきて面白かったです。

ぺル方程式(Problem 66)

簡単には

Dが平方数でない整数の時
x^2 - Dy^2 = 1
の整数解を求めよ

という問題。これはぺル方程式といって、ディオファントス方程式の一種です。
以下のページに詳しく理論と解法が紹介されています。連分数展開を使います。
http://aozoragakuen.sakura.ne.jp/suuron/suuron.html

余談ですが、ディオファントス方程式(特にベズー方程式 a x + b y = d)はコンパイラの教科書にも出てきまして、その関係で私は前からぺル方程式の事を知っていました。
二つの配列アクセスa[p*i+q]とb[r*j+s]が同じ場所を参照するかどうかは

a + sizeof(a[0]) * (p * i + q) == b + sizeof(b[0]) * (r * j + s)

が解を持つかどうかで決まります。これはi,jを変数とするベズー方程式になっています。この解を調べる事でループの並列性検査などが出来ます。

オイラーのトーティエント関数(Problem 69, 70)

オイラーのトーティエント関数φ(n)とはnと互いに素なn以下の自然数の数として定義されています。
これはnの素因数分解が
n = p_1^{m_1}p_2^{m_2}\cdots p_n^{m_n}
となるとき、
\varphi(n) = n\frac{p_1-1}{p_1}\frac{p_2-1}{p_2}\cdots\frac{p_n-1}{p_n}
を満たします。nの素因数分解は非常に重いですから、トーティエント関数に関する問題は

  1. 素数のリストの上で探索する
  2. エラトステネスの篩で計算

といった方法を取ることになると思います。例えば後者に関しては以下のようなコードで2〜Nについて関数の値を求める事が出来ます。

phi = (0..N).to_a
2.upto(N) do |n|
    if phi[n] == n # n is a prime number
        n.step(N, n) {|p| phi[p] = phi[p]*(n-1)/n}
    end
end

ファレイ数列、ブロコット木(Problem 71,72,73)

この3問はファレイ数列についての問題です。
ファレイ数列というのは0以上1以下の既約な分数を並べたものとして定義されます。
例えば分母が8以下の時は

F8 = 0/1, 1/8, 1/7, 1/6, 1/5, 1/4, 2/7, 1/3, 3/8, 2/5, 3/7, 1/2, 4/7, 3/5, 5/8, 2/3, 5/7, 3/4, 4/5, 5/6, 6/7, 7/8, 1/1

となります。

ファレイ数列の長さ

Problem 72はファレイ数列の長さを求める問題ですが、これはトーティエント関数を用いて
|F_n| = 1 + \sum_{k=1}^n\varphi(k)
と表せます。分母がkである既約な1未満の分数は、kと互いに疎なk以下の自然数の数と一致するのでこれが成立します。

ブロコット木

Problem 71と73はファレイ数列の一部だけを考える問題ですが、こういう問題ではブロコット木(Stern Brocot tree)が便利です。
ブロコット木というのは
既約な分数の分子同士、分母同士を加える事を繰り返して作られる木で図のようになります(wikipediaより引用)。
\frac{a}{b},\quad\frac{c}{d}\rightarrow\frac{a+c}{b+d}
http://upload.wikimedia.org/wikipedia/commons/thumb/3/37/SternBrocotTree.svg/2000px-SternBrocotTree.svg.png
このブロコット木は

  • 生成される分数は必ず既約になっている
  • 全ての分数を網羅する

という非常に良い性質を持っています。分母分子の約分の問題を考えなくて良いので、既約分数に関する探索問題が高速に解けます。

Problem 71: 分母が1,000,000以下の分数で3/7の未満の最大の物は何か？

3/7を挟む区間でブロコット木を辿って行けば良いです。

la, lb = 0, 1
ra, rb = 1, 2
while true
    ma, mb = la + ra, lb + rb
    if 7*ma < 3*mb
        if mb > 1000000
            puts "#{la}/#{lb}"
            exit
        end
        la, lb = ma, mb
    else
        ra, rb = ma, mb
    end
end

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

Problem 73: 分母が12,000以下の分数で1/2と1/3の間の数は何個か？

木のノード数を深さ優先探索で求めるだけ。また、分子の数字は考える必要ないです。

def dfs(l, r)
    m = l + r
    return 0 if m > 12000
    dfs(l, m) + dfs(m, r) + 1
end
puts dfs(3, 2)

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

ここら辺の話はKnuthの本に書いてあったと思います。

コンピュータの数学

コンピュータの数学