2010-12-18

Project Euler (Problem 3〜10)

数学

8問ときました。

Problem 3

600851475143の最大素因数を求めよ

「素因数分解せよ」ではなく「最大素因数だけ求めよ」ですから素数のテーブルは作りたくないので、反復により求める事を考えます。
整数Nを素因数分解してみると
$N = p_1^{m_1}p_2^{m_2}\cdots p_n^{m_n}\qquad(p_1\le p_2\le\cdots p_n)$
となりますが、これの最小の素数を全部取り除けば、以下の右辺の最大素因数をもとめることと同じです。
$\frac{N}{p_1^{m_1}} = p_2^{m_2}\cdots p_n^{m_n}$
これをどんどん繰り返していって、最後に残ったものが求める最大素数です。
で、ある素数pを取り去った後次に取り去る素数qですが、これはp+2,p+4,..と試して言って最初に割り切れた奴になります。なぜならば、qより小さい素数はすでに全部なくなっているからです。従ってコードは次の様になります。

n = 600851475143
k = 3
while n > k do
    while n % k == 0 do
        n /= k
    end
    k += 2
end
puts k

Problem 4

３桁の整数の積からなる最大の回文数を求めよ

3桁の整数の積は5 or 6桁ですからまず6桁の時を調べます。
3桁の整数の積を列挙するというのも手ですが、数え上げの基本は「条件の厳しい物を良く調べる」ですので、回文数の方を先に調べるのが筋です。最大の物を求めるので、大きい方から順に列挙すると

999999, 998899, 997799, 996699, 995599, 994499, 993399, ..., 102201, 101101, 100001

となりますので、綺麗に下降している上3桁に注目すると良いでしょう。これをkとおきます。
kの100,10,1の位はそれぞれ
$a = \lfloor k/100 \rfloor \\b = \lfloor k/10 \rfloor - \lfloor k/100 \rfloor 10 \\c = k - \lfloor k/10 \rfloor$
となるので、これから回文数の一般項が作れます。
$a10^5+b10^4+c10^3+c10^2+b10+a\\= 11(100k-90\lfloor k/10 \rfloor-9\lfloor k/100 \rfloor)$
6桁の回文数は11の倍数である事が分かります。一般的には偶数桁の回文数はかならず11の倍数です。
100k-...の部分ををAとおいて、2つの整数を11pとqとすると、桁数の制限から
$pq = A\\10^2 \le 11p < 10^3 \\10^2 \le q < 10^3 \\10^5 \le 11A < 10^6$
となるので、以下を満たせばOKです。
$\frac{A}{10^3} < p \le \frac{10^3}{11}$
以上より解答は以下のようになります。

999.downto(100) do |k|
    a = 100*k-90*(k/10)-9*(k/100)
    ((a+999)/1000..90).each do |p|
        if a % p == 0
            puts "#{11*a} = #{11*p} * #{a/p}"
            exit
        end
    end
end

Problem 5

1 から 20 の整数全てで割り切れる最小の整数を求めよ

つまり最小公倍数の事だから必要な素数を列挙すると

2は16=2^4なので最大4つ
3は9=3^2なので最大2つ
5,7,11,13,17,19は2乗すると20を超えるから最大1つ

という訳で

2^4*3^2*5*7*11*13*17*19

が答え

Problem 6

(1+2+...+100)^2-1^2+2^2+3^2+...+100^2を求めよ

和の公式を使って
$\left( \sum_{0< k\le n}k \right)^2 - \sum_{0< k\le n}k^2 \\= \frac{1}{4}n^2(n+1)^2-\frac{1}{6}n(n+1)(2n+1) \\\frac{1}{12}n(n-1)(n+1)(3n+2)$
となるのでn=100を代入して終わり。

Problem 7

10001番目の素数を求めよ

エラトステネスのふるいで良いわけですが、ふるいのサイズを決めなければなりません。こういうときは素数定理が役に立ちます。
今回は素数の値の上限を調べればよいので、以下の式を使います。
n番目の素数をp(n)とするとnが6以上の時は
$p(n) < n(\ln n + \ln \ln n)$
で抑えられます。
またnが20以上の時は
$p(n) < n(\ln n + \ln \ln n - \frac{1}{2})$
で抑えられます。

include Math
n = 10001
N = ( n * (log(n) + log(log(n)) - 1/2) ).floor
sieve = []
count = 0
2.upto(N) do |i|
    next if sieve[i]
    if (count += 1) == n then
        puts i;
        exit
    end
    (2*i).step(N, i) {|p| sieve[p] = true }
end

実行時間は0.15秒くらいだったのでこれで良しとします。
素数定理は計算量を推定する上でも役に立つので是非覚えておくと良いかと思います。

Problem 8

これはおもしろくないので省略

Problem 9

a < b < c かつ a^2+b^2=c^2 かつ a+b+c=1000 を満たす整数a,b,cを求め、abcを答えよ

1文字消去できるのでこれは高校１年くらいで学ぶ2元2次不定方程式の問題。現役高校生ならすぐできるでしょう。
まずはcを消去して定石通り因数分解します。
$a^2+b^2 = (1000-a-b)^2 \\\leftrightarrow (1000-a)(1000-b) = 500000$
今0 < a < bだから1000-a > 1000-bに注意すると
$707 < 1000-a < 1000 \qquad (\sqrt{5000} = \frac{1000\sqrt{2}}{2} = \frac{1414.2...}{2} 707.\cdots)$
が必要です。この範囲での500000の約数は800しか存在しませんからa=200。残りは自動的に出ます。

Problem 10

2000,000以下の素数の和を求めよ

ふつうにエラトステネスのふるいでやります。

N = 2000_000
sieve = []
sum = 0
2.upto(N) do |i|
    next if sieve[i]
    sum += i
    (2*i).step(N, i) {|p| sieve[p] = true }
end
puts sum

nineties% time ruby problem10.rb
142913828922
ruby problem10.rb  3.92s user 0.04s system 99% cpu 3.968 total

よくある高速化としては偶数を最初から取り除いておくという手があります。

N = 2000_000
sieve = []
sum = 2
3.step(N,2) do |i|
    next if sieve[i/2]
    sum += i
    (3*i).step(N, 2*i) {|p| sieve[p/2] = true }
end
puts sum

nineties% time ruby problem10.rb
142913828922
ruby problem10.rb  2.44s user 0.03s system 99% cpu 2.475 total

1.5倍くらい速くなりました。

2010-12-18

Project Euler Problem 12

数学

これはかなりの良門ではないでしょうか。

約数の数が500以上になる最初の三角数を求めよ

三角数とは以下の式で書けるものです。
$T(n) = \frac{1}{2}n(n+1)$
Nの素因数分解が下のようになる時
$N = p_1^{m_1}p_2^{m_2}\cdots p_n^{m_n}\qquad(p_1\le p_2\le\cdots p_n)$
約数の個数は
$F(n) = (m_1+1)(m_2+1)(m_3+1)\cdots$
となります。また、互いに素な2数は指数が重複しないですから、
nが奇数なら
F(T(n)) = F(n)F<a href=
で偶数なら
F(T(n *1 = F(n/2)F(n+1)" src=http://www.texify.com/img/%5CLARGE%5C%21F%28T%28n%29%29%20%3D%20F%28n/2%29F%28n%2B1%29.gif align=center border=0>
です。従って、三角数の事はもう忘れて良くて、「自然数nのF(n)を求める方法」を考えれば良いです。
「よしじゃぁ、まず素数のテーブルを作って・・・」というのではダサいです。もっと調べて行きましょう。
まず、上で述べた事よりFは以下の漸化式を満たします。
$F(p_1^{m_1}\cdots p_{n-1}^{m_{n-1}}p_n^{m_n}) = F(p_1^{m_1}\cdots p_{n-1}^{m_{n-1}})(m_n+1)$
こういった素因数分解に基づく漸化式は、エラトステネスのふるいをちょっと変形させれば解けます。
従って後は、m_nを求める方法を考えれば良いです。整数nの素因数に素数pが何個含まれるかをG(n,p)と表すと、
G(np,p) = G(n,p)+1
という漸化式になりますから、小さいnから順番に求める事ができます。
以上をまとめると以下のコードになります。

N = 20000
sieve = Array.new(N+1, 1)
exp = Array.new(N+1)
2.upto(N) do |n|
    if sieve[n] == 1 then
        exp.fill(0)
        n.step(N, n) do |k|
            exp[k] = exp[k/n] + 1
            sieve[k] *= (exp[k] + 1)
        end
    end
    r = (n % 2 == 0) ? sieve[n/2]*sieve[n-1] : sieve[n]*sieve[(n-1)/2]
    if r >= 500 then
        puts n*(n-1)/2
        exit
    end
end

実行時間は0.16秒くらいでした。

*1:n+1)/2)" src=http://www.texify.com/img/%5CLARGE%5C%21F%28T%28n%29%29%20%3D%20F%28n%29F%28%28n%2B1%29/2%29.gif align=center border=0> で偶数なら

2010-12-17

Project Euler始めた

数学

Project Euler解き始めました。自分の数学の勉強も兼ねて、できるだけプログラムでの力技の計算は避けようと思っています。どうしてもプログラミングが必要なやつは、あとでrowlで書こうと思ってます。

Problem 1

1000未満の3または5の倍数の和を求めよ

集合の加法定理が使えます。(3の倍数の和)+(5の倍数の和) - (15の倍数の和)ですから

となります。和の公式使って一発です。

Problem 2

4000,000以下のフィボナッチ数のうち偶数の物の和を求めよ。

書き出してみると

1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, ...

なので3n項目が偶数になる事がわかります。帰納法で証明できます。
フィボナッチ数列の一般項は

ですから偶数の項は

です。等比数列の和の公式を使って

が求める和です。ここで一般項、和ともに第２項の絶対値は0.5未満ですから、第１項だけ計算してもっとも近い整数に丸めれば答えが出ます。
f(3n) <= 4000,000となるnも対数つかって求められます。従ってRubyでの解答は以下のようになります。

include Math
SQRT5 = sqrt(5)
PHI3 = ( (1 + SQRT5)/2 )**3

n = ( log(4000_000 * SQRT5) / log(PHI3) ).floor # 項数
puts ((PHI3**(n+1)-1) / (SQRT5 * (PHI3 - 1))).round  # 和

2010-12-16

母関数を使って括弧のパターンを求めてみた

数学

n個の対応する括弧のパターンの数を母関数を使って求めてみました。
これはid:nobsunさんが指摘なさっているようにカタラン数という数になるんですが、そういう前提知識なしで機械的に解きたいような場合に使えます。
いろんな教科書があると思いますが、私は以下の本で学びました。
参考書籍：

コンピュータの数学

作者: ロナルド・L.グレアム,オーレンパタシュニク,ドナルド・E.クヌース,Ronald L. Graham,Oren Patashnik,Donald E. Knuth,有沢誠,萩野達也,安村通晃,石畑清
出版社/メーカー: 共立出版
発売日: 1993/08
メディア: 単行本
購入: 5人クリック: 224回
この商品を含むブログ (38件) を見る

まず求めるパターンを全て書きだした形式的なべき級数を考えてみます。式の中の1は括弧を一個も使わないパターンのつもり。

T = 1 + () + ()() + (()) + ()()() + ()(()) + (())() + (()()) + ((())) + ()()()() + ...

ここで括弧一組をxと置き換えてみると、

T = 1 + x + x^2 + x^2 + x^3 + x^3 + x^3 + x^3 + x^3 + x^4 + ... = 1 + x + 2 x^2 + 5 x^3 + 14 x^4 + ...

という級数になりますが、このx^nの各係数が求める答えであるという事が分かるかと思います。このように求める数列を係数としてもつ関数を母関数といいます。この母関数は式の再帰的な構造に注目すると簡単に求める事ができます。

では順番にやっていきます。
とりあえず、xに置き換える前の式に戻ります。

T = 1 + () + ()() + (()) + ()()() + ()(()) + (())() + (()()) + ((())) + ()()()() + ...

ここで途中で複数に分割できるパターン("()()"とか"()(())"とか)を(交換出来ない)掛け算だと思って、左端をくくり出すと

T = 1 + () {1 + () + ()() + (()) + ()()() + ...} + (()) {1 + () + ()() + (()) + ()()() + ...} + (()()) {1 + () + ()() + ...} + ...

くくった結果再びTが現れたので整理すると

T = 1 + () T + (()) T + (()()) T + ...

となります。今度はTをくくると

T = 1 + {() + (()) + (()()) + ...} T

となります。ここで残った中括弧の中身はTの各項を括弧で囲んだ物であることがわかります。つまり

() + (()) + (()()) + ... = '(' {1 + () + ()() + (()) + ...} ')'

となってますから、以下のシンプルな式を得ます。

T = 1 + '(' T ')' T

'...'がなくなったら括弧をxに置き換えます。この置き換えでx^nの係数は変わらない事が分かると思います。

T = 1 + x T^2

これは2次方程式ですから解の公式で解けて

T = (1 - sqrt(1- 4 x)) / (2 x)

と母関数が求まります。コイツのx^nの係数が求める答えです。
一応こいつの無限級数展開をWolfram Alphaで聞いてみると

確かに各係数(1,1,2,5,14,42,..)が求める答えになっている事がわかります。この係数の一般項を具体的求めると、それがカタラン数である事が分かるのですがWolfram Alphaはそこまでは教えてくれませんでした。

一般項を求めないで、漸化式を導くという事をやっても良いかと思います。
一つ前の式に戻ります。

T = 1 + x T^2

これから漸化式が作れます。n組の括弧からなるパターンの数をT(n)とおくと

T(n) = T(0) + T(1) x + T(2) x^2 + T(3) x^3 + T(4) x^4 + ...

となってる訳ですから、両辺を展開してx^nの係数だけ比較してみると、n > 0の時は

T(n) = T(0) T(n-1) + T(1) T(n-2) + ... + T(n-1) T(0)

となります。カタラン数の漸化式を機械的に求める事が出来ました。
参考までにですが、右辺の形は畳込みと言って「母関数の積(T^2)の係数 = 母関数の係数の畳込み」という有名な性質です。他にも微分・積分やシフト演算などを用いる母関数についての公式がありますので、それを覚えとくとこういった計算が素早くできます。

さて、「漸化式」があるならば「動的計画法」が定石です。rubyで書いてみるとこんな感じになります。(n=100まで求めてみます)

t = []
t[0] = 1
(1..100).each do |n|
    x = 0
    (0..n-1).each do |k|
        x += t[k] * t[n-k-1]
    end
    t[n] = x
    puts t[n]
end

%nineties% time ruby kakko.rb
1
2
5
14
42
132
429
1430
4862
16796
58786
208012
742900
2674440
... 中略 ...
896519947090131496687170070074100632420837521538745909320
ruby kakko.rb  0.03s user 0.01s system 77% cpu 0.057 total

文字列集合を得る

ここでも上で求めた式が使えます

T = 1 + '(' T ')' T

Haskellでやってみようと思いますが、右辺に無限リストが2つあるので漸化式に直しとく必要があります。

T(0) = 1
T(n) = '(' T(0) ')' T(n-1) + '(' T(1) ')' T(n-2) + ... + '(' T(n-1) ')' T(0)

これを立式すると次のようになります。

pattern 0 = [""]
pattern n = [concat["(",l,")",r] | i <- [0..n-1], l <- pattern i, r <- pattern (n-i-1)]

もしくは無限リストとして自分自身を参照するようにするとこんな感じになります。

pattern = [""] : [[concat["(",l,")",r] | i <- [0..n-1], l <- pattern !! i, r <- pattern !! (n-i-1)] | n <- [1..]]

patternは括弧のつけ方の無限リストになっています。

%ghci
GHCi, version 6.12.1: http://www.haskell.org/ghc/  :? for help
Loading package ghc-prim ... linking ... done.
Loading package integer-gmp ... linking ... done.
Loading package base ... linking ... done.
Prelude> :l kakko.hs
[1 of 1] Compiling Main             ( kakko.hs, interpreted )
Ok, modules loaded: Main.
*Main> take 7 pattern
[[""],["()"],["()()","(())"],["()()()","()(())","(())()","(()())","((()))"],["()()()()","()()(())","()(())()","()(()())","()((()))","(())()()","(())(())","(()())()","((()))()","(()()())","(()(()))","((())())","((()()))","(((())))"],["()()()()()","()()()(())","()()(())()","()()(()())","()()((()))","()(())()()","()(())(())","()(()())()","()((()))()","()(()()())","()(()(()))","()((())())","()((()()))","()(((())))","(())()()()","(())()(())","(())(())()","(())(()())","(())((()))","(()())()()","(()())(())","((()))()()","((()))(())","(()()())()","(()(()))()","((())())()","((()()))()","(((())))()","(()()()())","(()()(()))","(()(())())","(()(()()))","(()((())))","((())()())","((())(()))","((()())())","(((()))())","((()()()))","((()(())))","(((())()))","(((()())))","((((()))))"],["()()()()()()","()()()()(())","()()()(())()","()()()(()())","()()()((()))","()()(())()()","()()(())(())","()()(()())()","()()((()))()","()()(()()())","()()(()(()))","()()((())())","()()((()()))","()()(((())))","()(())()()()","()(())()(())","()(())(())()","()(())(()())","()(())((()))","()(()())()()","()(()())(())","()((()))()()","()((()))(())","()(()()())()","()(()(()))()","()((())())()","()((()()))()","()(((())))()","()(()()()())","()(()()(()))","()(()(())())","()(()(()()))","()(()((())))","()((())()())","()((())(()))","()((()())())","()(((()))())","()((()()()))","()((()(())))","()(((())()))","()(((()())))","()((((()))))","(())()()()()","(())()()(())","(())()(())()","(())()(()())","(())()((()))","(())(())()()","(())(())(())","(())(()())()","(())((()))()","(())(()()())","(())(()(()))","(())((())())","(())((()()))","(())(((())))","(()())()()()","(()())()(())","(()())(())()","(()())(()())","(()())((()))","((()))()()()","((()))()(())","((()))(())()","((()))(()())","((()))((()))","(()()())()()","(()()())(())","(()(()))()()","(()(()))(())","((())())()()","((())())(())","((()()))()()","((()()))(())","(((())))()()","(((())))(())","(()()()())()","(()()(()))()","(()(())())()","(()(()()))()","(()((())))()","((())()())()","((())(()))()","((()())())()","(((()))())()","((()()()))()","((()(())))()","(((())()))()","(((()())))()","((((()))))()","(()()()()())","(()()()(()))","(()()(())())","(()()(()()))","(()()((())))","(()(())()())","(()(())(()))","(()(()())())","(()((()))())","(()(()()()))","(()(()(())))","(()((())()))","(()((()())))","(()(((()))))","((())()()())","((())()(()))","((())(())())","((())(()()))","((())((())))","((()())()())","((()())(()))","(((()))()())","(((()))(()))","((()()())())","((()(()))())","(((())())())","(((()()))())","((((())))())","((()()()()))","((()()(())))","((()(())()))","((()(()())))","((()((()))))","(((())()()))","(((())(())))","(((()())()))","((((()))()))","(((()()())))","(((()(()))))","((((())())))","((((()()))))","(((((())))))"]]

「日本の硬貨を使ってn円払う払い方の総数は？」

母関数の応用として有名な例を一つ紹介。

1 + x + x^2 + x^3 + ...

のx^nの係数が1円玉のみでn円払う払い方の総数(つまり1通り)。

1 + x^5 + x^10 + x^15 + x^20 + ...

のx^nの係数が5円玉のみでn円払う払い方の総数(つまり1通り)。
...
これらの組み合わせを考えれば良いので、母関数は

T = (1 + x + x^2 + x^3 + ..) (1 + x^5 + x^10 + x^15 + x^20 + ...) (1 + x^10 + x^20 + ...) (1 + x^50 + x^100 + ...) (1 + x^100 + x^200 + ...) (1 + x^500 + x^1000 + ...)

となります。それぞれの形式的級数(無限等比級数の和の公式を使って)を閉形式にすると

T = 1/(1-x) 1/(1-x^5) 1/(1-x^10) 1/(1-x^50) 1/(1-x^100) 1/(1-x^500)

となります(ここでは発散とか考えなくてよし)。
これから一般解を求める事ができますが、めんどくさいので省略。
漸化式を作ります。
分母を払うと

(1-x)(1-x^5)(1-x^10)(1-x^50)(1-x^100)(1-x^500)T = 1 .... (1)

となります。
さてTのx^nの係数T(n)の漸化式を求めます。
まず(1-x^500)Tのx^nの係数は

A(n) = T(n)-T(n-500)

です。(但し、X(n)の中身が負になったときはX(n)=0とします。)
よって(1-x^100)(1-x^500)Tのx^nの係数は

B(n) = A(n) - A(n-100)

です。これを繰り返すと

A(n) = T(n)-T(n-500)
B(n) = A(n) - A(n-100)
C(n) = B(n) - B(n-50)
D(n) = C(n) - C(n-10)
E(n) = D(n) - D(n-5)
F(n) = E(n) - E(n-1)

となります。
これを移項して整理すると

T(n) = T(n-500) + A(n)
A(n) = A(n-100) + B(n)
B(n) = B(n-50) + C(n)
C(n) = C(n-10) + D(n)
D(n) = D(n-5) + E(n)
E(n) = E(n-1) + F(n)

となって下から上に向かって動的計画法で解ける形になります。
F(n)が(1)式の左辺のx^nの係数なのでF(0) = 1で、n > 0の時はF(n) = 0となる事に注意してください。
rubyで書くと以下の様になります。せっかくのRubyなので「インデックスが負なら0」をArrayの[]を書き換えて実装してみました。

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

t = []
a = []
b = []
c = []
d = []
e = []

(0..10000).each do |n|
    f = (n == 0) ? 1 : 0
    e[n] = e[n-1] + f
    d[n] = d[n-5] + e[n]
    c[n] = c[n-10] + d[n]
    b[n] = b[n-50] + c[n]
    a[n] = a[n-100] + b[n]
    t[n] = t[n-500] + a[n]
    puts "#{n}円 : #{t[n]}通り"
end

nineties% time ruby kouka.rb
0円 : 1通り
1円 : 1通り
2円 : 1通り
3円 : 1通り
4円 : 1通り
5円 : 2通り
6円 : 2通り
7円 : 2通り
8円 : 2通り
9円 : 2通り
10円 : 4通り
11円 : 4通り
...

9996円 : 7825615000通り
9997円 : 7825615000通り
9998円 : 7825615000通り
9999円 : 7825615000通り
10000円 : 7844606371通り
ruby kouka.rb  0.19s user 0.06s system 63% cpu 0.395 total

この例では問題の難易度に比べて大げさすぎる気がしますが、まぁおもしろいんではないかと思います。他にもいろんな応用がありますので、是非調べてみて下さい。。

2010-12-13

rowlの式

rowl

rowlの式として以下の定義を使う方針にしました。その為、いろいろと実装の修正を行っています。

原始式(シンボル，整数リテラル，浮動小数点数リテラル，文字列)は式である。
hがシンボル，e1,...,enが式の時 h(e1, .., en)は式である。hを式のheadと言う。

この表現もLISPと同じくSymbol expression、つまりS式と呼ばれます。各々を区別する呼称があった気がするけど忘れたので、ここでは仮にσ式と呼ぶ事にしときます。
σ式とS式は本質的には同じで、

S式のドット対(a . b)はσ式のCons(a, b)で表現できる

と考えても良いし

σ式のh(e1, ..., en)はS式の(h e1 ... en)と表現できる

と考えても良いです。
σ式では例えば，1+2*3はPlus(1, Times(2, 3))の様な自然な表現になります。このようなヘッドを明に持った表現は構文木表現としてはS式よりも一般的で、それを表層のレベルでも採用している処理系としてはMathematicaが有名です。

S式にしろσ式にしろ、構文木もデータも同じ形で表現されるからいじくりまわすのが簡単であるという利点はそんなに変わりません。S式の方がより表現がシンプルなので扱い安いのですが、LISPの記法が嫌なのと、式の意味を表すヘッドが明示的である方が好みであることからσ式を採用する事にしました。

構文糖衣

組み込みの構文糖衣

[e1, ..., en] → List(e1, ..., en)
(e1, ..., en) → Tuple(e1, ..., en)
{e1 ... en} → Block(e1, ..., en)
e1(e2, ..., en)
- e1がシンボルの時はそのまま
- それ以外はApply(e1, List(e2, ..., en))
e1[e2, ..., en] → Subscr(e1, List(e2, ..., en))

前置・後置・中置演算子とコンストラクタは組み込みは無しでユーザ側・ライブラリ側で定義。

2010-11-28

パターンマッチの実装(3)

rowl

型でのマッチングの実装を一部実装しました。
これによって関数のオーバーロードが出来るようになりました↓

print(10)            # print_intが呼ばれる
print("Hello World") # print_stringが呼ばれる

関数定義時にはこんな感じで書きます↓

fib(n!Int): fib(n-1) + fib(n-2)
fib(0): 0
fib(1): 1

これでInt以外でfibを呼ぶとパターンマッチエラー(型エラーではなく)になります。これまでと同じく、型のパターンマッチのコードもJIT コンパイルで生成しています。

静的型付け言語と違って、rowlでは型によるマッチングであっても実行時検査で行うので、実行時オーバヘッドがあります。
例えば、上の定義でfib(36)を走らせた場合は、4.66秒だったのが6.91秒と遅くなりました。今回整数演算でもパターンマッチを行うようしたので、だいぶ大きなオーバヘッドになっています。

一方、型パターン使用時はJIT コンパイラが型の情報を使って最適化をかける事が出来るので、最終的に性能的にもプラスにはたらくことを期待しています。rowlはRuby等と違ってグローバルに関数定義情報を持っていて、高速化が楽になるように作っています。

2010-11-23

最適化コンパイラ屋がHaskellの副作用について考えてみる

Haskellに副作用があるのか？というのは難しいテーマだと思いますが、少なくとも最適化理論での一般的な「副作用」の定義ではHaskellは全く副作用がない言語と言えると思います。
理論的に美しいという点がHaskellの設計の一番重要なところだと思うのですが、バックエンドの泥臭い話も面白いかもと思ったので書いてみます。
私は主に手続き型言語の最適化をやっていて関数型言語は詳しくありませんので、見当違いな事を書いていたら指摘して下さい。

まずは副作用の定義の為に「データ依存関係」(Wikipedia:依存関係)の説明をします。
ある変数xを使用・定義する二つの文を考える時、その評価順序関係には以下の4種類があり、それぞれの関係を以下のように呼びます。

真の依存

x = ...
...
... = x

逆依存

... = x
...
x = ...

出力依存

x = ...
...
x = ...

入力依存

... = x
...
... = x

入力依存は今回の話と関係ありません。
前者３つの依存関係は「順序を入れ替えるとプログラムの意味が変わっていまう」関係です。
最適化はプログラムの意味を変えてはならないというのが大前提なので、これら依存関係を解析した上でプログラムの変換をしなければなりません。

ここで、上の変数xの使用・代入からなる依存は分かりやすいので「自明な依存」と呼ぶ事にします。
さて、「自明でない依存」というのもあります。
例えばポインタを使うケース：

*p = ...;
... = x;

これはもし(p == &x)なら依存関係にあります。
次は、集成体を使うケース(例えば配列の場合):

for (i = 0; i < n; i++) {
    a[i+d] = a[i];
}

これは0 < d < nの時のみ並列化できないループです。
そして、グローバルな状態を変更するケース:

f();
g();

これはf()とg()が同じグローバル変数にアクセスするなどの場合に依存関係にあります。

上で述べた自明な依存に加えてこれらの自明でない依存も考慮しなければ正しいプログラム変換ができません。
最適化理論ではこのような「自明でない依存」を引き起こす可能性のある命令が「副作用」を持つと言う事が多いと思います。
ただし、「参照透過的である」のような統一的な定義があるわけでもないですし、自明でない依存もポインタ解析・シンボリック解析・関数間解析などで頑張って調べれば部分的に判別できる場合もあります。

ここでは、仮に「データ依存以外に評価順序を決める要因」全てを「副作用」と定義した場合の話をします。

さて、入出力、グローバル変数への破壊的代入、書き換え可能な集成体・・・等の機能があると*普通は*この手の自明でない依存が生じる事が分かると思います。
Haskellでも入出力は可能だし、グローバル変数への破壊的代入や集成体の書き換えも(メモリの上書きという意味では)可能です。しかし、自明でない依存は存在しません。
例えば以下のコードを見て下さい。

main = do x <- getLine
          y <- getLine
          putStrLn x
          putStrLn y

このレベルで見ると、x <- getLineとy <- getLineの間にはデータ依存がないのにこれらが入れ替えられない事から、一見自明でない依存があるように見えます。
しかしコンパイラが実際に扱うのはこれを脱糖したコードで、おおざっぱに書くと下の様な感じの物になります。(getline,putstrはgetLineやputStrの実際の仕事をするプリミティブのつもり)。

main = IO (\s1 ->
    let (s2, x) = getline s1 in
    let (s3, y) = getline s2 in
    let (s4, _) = putstr x s3 in
    putstr y s4
    )

前で定義されたsNを次の命令に渡している事から、これらの間には真の依存があります。これは「自明な依存」です。

つまり、Haskellは自明でない依存を、上のように変数を引き渡す事で自明な依存に直しているのです。
do記法という構文糖衣の上から見るとこの変数は見えないので、人間側からはあたかも副作用が存在しているかのように見える訳です。

というわけで、Haskellは「真の依存」のみによって評価順序が決定されるので、例えばラムダ計算の様な体系にそのまま乗っける事が出来るし、難しい事考えなくても項書き換えシステムで正しく動作させる事が出来るし、自由自在にプログラムを変換してもその意味が決して変わらないのでアグレッシブな最適化が出来るわけです。
副作用がないというのはユーザだけでなくコンパイラにとっても大きなメリットなのです。

一方で、デメリットは

ユーザにとってはプログラムの字面から評価順序を判別する事が難しくなる
本来依存関係にない命令間に依存が生じる(上の例だとy <- getLineとputStrLn xの間には本来は依存がないのに、データ依存が生じている)

などでしょうか。特に後者の為に、思ったほどはHaskellの並列化は楽ではなかったりするのだと思います。

まとめると

入出力や破壊的代入などは普通は自明でない依存を生む
Haskellではそれらが自明な依存としてコードに現れる
Haskellにはコード変換時に考慮すべき副作用が存在しない

という話でした。