Haskell

2013年11月19日 (火)

フィボナッチ数・いろいろ

《 はじめに … フィボナッチ数の定義 》

またまたフィボナッチ数の話題です(笑)

前回の記事にも書きましたが、フィボナッチ数の定義は次のようになります。

\begin{eqnarray*} \begin{cases} F_{0} &=& 0 \\ F_{1} &=& 1 \\ F_{n+2} &=& F_{n+1} + F_{n} \end{cases} \end{eqnarray*}

フィボナッチ数の求め方は、これまでもいろいろと考えられてきたようです。 今回は、私の知っているいくつかの方法を紹介するとともに、実行速度についても比較してみようと思います。 (速度の比較は "GHCi" 上で行なっているので、コンパイルして実行した場合とは少し違っているかもしれません)

《 再帰的プロセス 》

まずは、フィボナッチ数の定義そのままのもの。

fib1 :: Int -> Integer
fib1 0 = 0
fib1 1 = 1
fib1 n = fib1 (n - 1) + fib1 (n - 2)
*Main> mod (fib1 10) 3
1
(0.00 secs, 0 bytes)
*Main> mod (fib1 20) 3
0
(0.05 secs, 4263488 bytes)
*Main> mod (fib1 30) 3
2
(1.87 secs, 223412808 bytes)
*Main> mod (fib1 35) 3
2
(20.30 secs, 2456928716 bytes)

この関数はとても時間がかかります。\(F_{35}\) を求めるのに 20 秒以上かかってしまいました。

SICP で詳しく解説されいるのですが、関数 fib1 は呼ばれるたびに自分を二回呼び出すので、とても非効率的です。

《 反復的プロセス 》

フィボナッチ数は反復的プロセスでも求めることができます。

fib2 :: Int -> Integer
fib2 n = loop n 0 1 
  where
    loop 0 a _ = a
    loop n a b = loop (n - 1) b (a + b)

fib2 の内部関数である "loop 関数" は初期値として \(F_{0}\) と \(F_{1}\) を持ち、それをもとに次々に \(F_{2}\), \(F_{3}\) … を計算していきます。 求めるフィボナッチ数を \(F_{n}\) とすると、この "loop" は n 回繰り返されます。ですから計算量は Θ(n) ということになると思います。

*Main> mod (fib2 1000) 3
0
(0.02 secs, 2136108 bytes)
*Main> mod (fib2 10000) 3
0
(0.05 secs, 8659128 bytes)
*Main> mod (fib2 100000) 3
0
(1.58 secs, 611890040 bytes)
*Main> mod (fib2 1000000) 3
0
(336.62 secs, 45991518544 bytes)

「再帰的プロセス」に比べれば高速ですが、まだまだ遅いですね。

《 行列を使ったもの 》

Wikipedia にも載っているのですが、フィボナッチ数の漸化式は行列でも表現できます。

\[ \left( \begin{array}{ll} F_{n + 1} & F_{n} \\ F_{n} & F_{n - 1} \end{array}\right) = \left( \begin{array}{ll} 1 & 1 \\ 1 & 0 \end{array} \right) ^ {n} \]

前回の記事で使用した "power 関数" を使用すれば行列の計算方法を "power 関数" に教えてやるだけで OK です。

行列 \( \left( \begin{array}{ll} a & b \\ c & d \end{array} \right) \) を、 Haskell のタプルを使って "(a b c d)" と表現すると、行列の積を計算する「関数 f」は次のように定義できます。

f (a1, b1, c1, d1) (a2, b2, c2, d2) = (a', b', c', d')
  where
    a' = a1 * a2 + b1 * c2
    b' = a1 * b2 + b1 * d2
    c' = c1 * a2 + d1 * c2
    d' = c1 * b2 + d1 * d2

これを使ってフィボナッチ数を求める関数を定義してみます。

power :: Integral a => (t -> t -> t) -> t -> a -> t -> t
power _ _ 0 e = e
power f a n e = power f (f a a) (div n 2) (if odd n then f a e else e)

fib3 :: Int -> Integer
fib3 n = (power f [1,1,1,0] n [1,0,0,1]) !! 1

f [a1, b1, c1, d1] [a2, b2, c2, d2] = [a', b', c', d']
  where
    a' = a1 * a2 + b1 * c2
    b' = a1 * b2 + b1 * d2
    c' = c1 * a2 + d1 * c2
    d' = c1 * b2 + d1 * d2
*Main> mod (fib3 100000) 3
0
(0.03 secs, 1940120 bytes)
*Main> mod (fib3 1000000) 3
0
(0.09 secs, 7170248 bytes)
*Main> mod (fib3 10000000) 3
0
(1.09 secs, 87110184 bytes)
*Main> mod (fib3 100000000) 3
0
(14.41 secs, 917842092 bytes)

"power 関数" は「n 乗」を Θ(log n) で計算してくれるので、これまでの方法と比べてかなり高速です。

《 SICP で紹介されているアルゴリズム 》

こちらのブログ に次のような、 SICP に載っているアルゴリズムを Haskell に移植したものが紹介されていま す(関数名は変えさせてもらってます)。

fib4 :: Int -> Integer
fib4 = iter 1 0 0 1
  where
    iter a b p q count
      | count == 0 = b
      | even count = iter a b (p*p+q*q) (2*p*q+q*q) (count `div` 2)
      | otherwise  = iter (b*q+a*q+a*p) (b*p+a*q) p q (count - 1)

SICP によると、この関数は Θ(log n) で \(F_{n}\) を求めることができるそうです。

*Main> mod (fib4 100000) 3
0
(0.00 secs, 0 bytes)
*Main> mod (fib4 1000000) 3
0
(0.09 secs, 5390084 bytes)
*Main> mod (fib4 10000000) 3
0
(0.78 secs, 68803364 bytes)
*Main> mod (fib4 100000000) 3
0
(11.31 secs, 766052000 bytes)

計算量が Θ(log n) なので、速度は fib3 とほぼ同様になります。

《 Gosper and Salamin のアイデア 》

ネットでフィボナッチ数についていろいろ調べている時に、『フィボナッチ数とリュカ数の計算法』という文献を見つけました。その中に「Gosper and Salamin のアイデア」というものが載っていました。

詳しくはこの文献を読んでいただくとして、要約すると、「順序対 (a, b) と (c, d) に対して乗法を (a, b)(c, d) = (ac + ad + bc, ac + bd) と定義すると、この乗法は結合律を満し、指数法則が成り立つ。そして f ≡ (1, 0) とすると \((F_{n}, F_{n + 1})\) = f n で表わすことができる ということのようです。これは SICP に載っている fib4 のと同じ方法のように見えます。

これを自作の "power 関数" を使って定義すると、次のようになりました。

power :: Integral a => (t -> t -> t) -> t -> a -> t -> t
power _ _ 0 e = e
power f a n e = power f (f a a) (div n 2) (if odd n then f a e else e)

fib5 :: Int -> Integer
fib5 n = fst $ power gsMulti (1, 0) n (0, 1)
  where gsMulti (a, b) (c, d) = (a * (c + d) + b * c, a * c + b * d)
*main> mod (fib5 100000) 3
0
*Main> :set +s
*Main> mod (fib5 100000) 3
0
(0.02 secs, 1686056 bytes)
*Main> mod (fib5 1000000) 3
0
(0.09 secs, 6295976 bytes)
*Main> mod (fib5 10000000) 3
0
(0.81 secs, 68134876 bytes)
*Main> mod (fib5 100000000) 3
0
(10.89 secs, 727571692 bytes)

"power 関数" を使っているので、計算量はやはり Θ(log n) になります。

《 一般項の式を使ったもの 》

最後に前回の記事に書いた「一般項の式を使ったもの」です。

import Data.Ratio

power :: Integral a => (t -> t -> t) -> t -> a -> t -> t
power _ _ 0 e = e
power f a n e = power f (f a a) (div n 2) (if odd n then f a e else e)

fib6 :: Int -> Integer
fib6 n = truncate . (* 2) . snd $ power f (1 % 2, 1 % 2) n (1, 0)
  where f (a1, b1) (a2, b2) = (a1 * a2 + 5 * b1 * b2, a1 * b2 + a2 * b1)
*Main> mod (fib6 100000) 3
0
(0.03 secs, 2389868 bytes)
*Main> mod (fib6 1000000) 3
0
(0.11 secs, 8976944 bytes)
*Main> mod (fib6 10000000) 3
0
(0.94 secs, 106752416 bytes)
*Main> mod (fib6 100000000) 3
0
(12.31 secs, 1174476528 bytes)

"power 関数" を使っているので、計算量は Θ(log n) のはずですが、有理数の計算をしてる分、"fib4" や "fib5" よりは少し遅くなるようです。

《 まとめ 》

これまでの結果をまとめてみます。

関数n時間 (sec)
fib13520.30
fib21,000,000336.62
fib3100,000,00014.41
fib4100,000,00011.31
fib5100,000,00010.89
fib6100,000,00012.31

今回はメモ化などの手法は考えなかったので、結果は「計算量の少なさ=計算速度」ということになっています。 fib3, fib4, fib5, fib6 はいずれも計算量が Θ(log n) なので、かなり速くフィボナッチ数を求めるこができます。

2013年11月14日 (木)

フィボナッチ数を一般項から求めてみる。

フィボナッチ数を一般項から求めてみる

《 はじめに 》

フィボナッチ数とは、

\begin{eqnarray*} \begin{cases} F_{0} &=& 0 \\ F_{1} &=& 1 \\ F_{n+2} &=& F_{n+1} + F_{n} \end{cases} \end{eqnarray*}

という式で定義される数ですが、次のような一般項を表わす式もわかっています。

\[ F_{n} = \frac{1}{\sqrt{5}} \left\{ \left( \frac{1 + \sqrt{5}}{2} \right) ^ n - \left( \frac{1 - \sqrt{5}}{2} \right) ^ n \right\} \]

今回はこの一般項の式を用いてフィボナッチ数を求める方法を考えてみました。


《 一般項のとおりに計算してみると… 》

「一般項の式があるのなら、それをそのまま使えばいいのでは…」と考える方もいると思います。
では、試してみましょう。

phi  = (1 + sqrt 5) / 2
phi' = (1 - sqrt 5) / 2

fib1 :: Int -> Integer
fib1 n = truncate $ (phi ^ n - phi' ^ n) / (sqrt 5)

この関数を "GHCi" に読み込んで実行してみます。

*Main> map fib1 [1..10]
[1,1,2,3,5,8,13,21,34,54]

F10 (10番目のフィボナッチ数) は "55" のはずなのに、"fib1 10" の結果は "54" になってしまいました。
これは fib1 が浮動小数点数で計算した値を整数に変換しているため、計算誤差が生じていると考えられます。

正しいフィボナッチ数を求めるためには、他の方法を考えないといけないようです。


《 一般項の式をよく見てみると… 》

一般項の式を使って、実際にいくつかのフィボナッチ数を計算してみましょう。

\begin{eqnarray*} & F_{1} = & \frac{1}{\sqrt{5}} \left\{ \left( \frac{1}{2} + \frac{1}{2} \sqrt{5} \right) - \left( \frac{1}{2} - \frac{1}{2} \sqrt{5} \right) \right\} & = & \frac{1}{\sqrt{5}} ・ \frac{2}{2} \sqrt{5} & = & 1 \\ & F_{2} = & \frac{1}{\sqrt{5}} \left\{ \left( \frac{6}{4} + \frac{2}{4} \sqrt{5} \right) - \left( \frac{6}{4} - \frac{2}{4} \sqrt{5} \right) \right\} & = & \frac{1}{\sqrt{5}} ・ \frac{4}{4} \sqrt{5} & = & 1 \\ & F_{3} = & \frac{1}{\sqrt{5}} \left\{ \left( \frac{16}{8} + \frac{8}{8} \sqrt{5} \right) - \left( \frac{16}{8} - \frac{8}{8} \sqrt{5} \right) \right\} & = & \frac{1}{\sqrt{5}} ・ \frac{16}{8} \sqrt{5} & = & 2 \\ & F_{4} = & \frac{1}{\sqrt{5}} \left\{ \left( \frac{56}{16} + \frac{24}{16} \sqrt{5} \right) - \left( \frac{56}{16} - \frac{24}{16} \sqrt{5} \right) \right\} & = & \frac{1}{\sqrt{5}} ・ \frac{48}{16} \sqrt{5} & = & 3 \\ & F_{5} = & \frac{1}{\sqrt{5}} \left\{ \left( \frac{176}{32} + \frac{80}{32} \sqrt{5} \right) - \left( \frac{176}{32} - \frac{80}{32} \sqrt{5} \right) \right\} & = & \frac{1}{\sqrt{5}} ・ \frac{160}{32} \sqrt{5} & = & 5 \end{eqnarray*}

こうして見ていくと、二項定理からも明らかなように \(( 1 / 2 + \sqrt{5} / 2 )^{n} - ( 1 / 2 - \sqrt{5} / 2 )^{n}\) を計算すると、有理数は打ち消しあって無理数( \(\sqrt{5}\) を含む項 )だけが残ります。
従って、\(( \frac{1}{2} + \frac{\sqrt{5}}{2} ) ^ {n}\) の \(\sqrt{5}\) の係数を2倍すると \(F_{n}\) と等しくなります。

では、どうしたら \(( \frac{1}{2} + \frac{\sqrt{5}}{2} ) ^ {n}\) の \(\sqrt{5}\) の係数を求められるのでしょうか?

\(「(有理数 + 有理数・\sqrt{5})^{n}」\) は必ず \(「有理数 + 有理数・\sqrt{5}」\) の形になるので、 こちらの記事 にも書いてあるように、虚数の計算のように「有理数」と \(「有理数・\sqrt{5}」\) をわけて計算すれば良さそうです。 例えば \(「5 + 3\sqrt{5}」\) を (5, 3) と表現することにします。


《 下準備 》

ここで、ちょっと寄り道をします。

SICP の中で、 bn を高速に計算する方法として、次のような手続きが紹介されています。

(define (fast-expt b n)
  (cond ((= n 0) 1)
        ((even? n) (square (fast-expt b (/ n 2))))
        (else (* b (fast-expt b (- n 1))))))

でも、この方法は「乗算」だけに限定する必要はないんです。結合則を満す演算なら何でも使えそうです。
そこで『素因数分解と素数判定』という本に載っていたアルゴリズムを参考に次のような関数を作ってみました。

-----------------------------------------
-- 冪乗法
-----------------------------------------
-- e の初期値 : 関数 f の単位元
-- ex : 2^100 == power (*) 2 100 1
--
power :: Integral a => (t -> t -> t) -> t -> a -> t -> t
power _ _ 0 e = e
power f a n e = power f (f a a) (div n 2) (if odd n then f a e else e)

この関数を使うと \(「a^{n}」\) だけでなく、 素数判定などにでてくる \(「a^{n} \pmod{m}」\) のような計算も高速に行えます。

--
-- a^n (mod m) を計算する
--
powerMod :: (Integral a, Integral t) => t -> a -> t -> t
powerMod a n m = power (\x y -> mod (x * y) m) a n 1




《 一般項からフィボナッチ数を求める関数 》

話をフィボナッチに戻します。

\(a_{1}, a_{2}, b_{1}, b_{2}\) を有理数とすると、

\[( a_{1} + b_{1} \sqrt{5} ) ( a_{2} + b_{2} \sqrt{5} ) = ( a_{1}a_{2} + 5b_{1}b_{2} ) + ( a_{1}b_{2} + a_{2}b_{1} ) \sqrt{5} \]

となることから、次のような関数を定義します。

f (a1, b1) (a2, b2) = (a1 * a2 + 5 * b1 * b2, a1 * b2 + a2 * b1)

この関数と前述の「power 関数」を使えば、次のような「一般項からフィボナッチ数を求める関数」を定義することができます(有理数を扱うので「import Data.Ratio」が必要です)。

import Data.Ratio

power :: Integral a => (t -> t -> t) -> t -> a -> t -> t
power _ _ 0 e = e
power f a n e = power f (f a a) (div n 2) (if odd n then f a e else e)

fibonacci :: Int -> Integer
fibonacci 0 = 0
fibonacci n = truncate . (* 2) . snd $ power f (1 % 2, 1 % 2) n (1, 0)
  where f (a1, b1) (a2, b2) = (a1 * a2 + 5 * b1 * b2, a1 * b2 + a2 * b1)

「power 関数」がかなり強力なので、この「fibonacci 関数」も高速です。

GHCi 上で実行したら、次のような結果になりました。 (実際に計算させると 20898764 桁の数になるので、結果の表示に時間がかからないように "mod" を付け加えました)。

> mod (fib9 100000000) 3
0
(12.32 secs, 1175839488 bytes)

100,000,000 番目のフィボナッチ数が 12 秒ほどで計算できるのだから、なかなかのものではないでしょうか。

2012年10月19日 (金)

セクシー素数

先日、「Haskellで見る双子素数, いとこ素数, セクシー素数」というブログを見つけました。

以前から素数に興味があり、プログラミング言語を覚える度に、素数を求めるメソッドや関数を作ったりしてました。
「双子素数」や「三つ子素数」は知ってたんですけどね。「セクシー素数」は知りませんでした。

こちらのブログには「双子素数」「いとこ素数」「セクシー素数」などを求めるHaskell のコードが載っています。
例えば「双子素数」や「いとこ素数」は

twins = [(p1, p2) | (p1, p2) <- zip primes (tail primes), p1 + 2 == p2]
cousins = [(p, p+4) | p <- primes, isPrime (p+4)]

といった具合です。(primes と isPrime は Data.Numbers.Primes というパッケージが使用されているようです)

「双子素数」のコードに異論はないのですが、「いとこ素数」のコードはちょっと気になるところがありました。
せっかく primes という素数列が手元にあるのだから、「いとこ素数」のペアを探すのにわざわざ isPrime 使う必要はないと思うんです。
言葉だけでは分かりにくいので具体的なコードを書いてみます。(以下のコードではすべて、私が考えた primes と isPrime のコードを使用しています)

-- 素数列
-- http://d.hatena.ne.jp/notogawa/20110114/1295006865 を参考に
-- gcd の使い方が秀逸
primes :: Integral a => [a]
primes = map fromIntegral primes'
    where
      primes' = [2, 3, 5] ++ f 5 7 (drop 2 primes')
      f m s (p : ps) = [n | n <- ns, gcd m n == 1] ++ f (m * p) (p * p) ps
          where ns = [x + y | x <- [s, s + 6 .. p * p - 2], y <- [0, 4]]

isPrime :: Integral a => a -> Bool
isPrime n = (n > 1) && isPrime' primes
    where
      isPrime' (p : ps)
          | p * p > n    = True
          | rem n p == 0 = False
          | otherwise    = isPrime' ps

cousinPrimes :: Integral a => [(a, a)]
cousinPrimes = loop primes
  where
    loop (p1 : ps1)
      | p2 - p1 == 4 = (p1, p2) : loop ps1
      | otherwise     = loop ps1
      where (p2 : _) = dropWhile (< p1 + 4) ps1

見てもらいたいのは、 cousinPrimes の

(p2 : _) = dropWhile (< p1 + 4) ps1

の部分です。

dropWhile で (p1 + 4) 未満の素数を落していくと、p2 には当然 (p1 + 4) 以上の素数のうち一番小さい素数が入ります。(p1 + 4) が素数なら、その数が必ず p2 になります。
もし (p1 + 4) が素数でなければ p2 には (p1 + 4) より大きい素数が入るので、p1 と p2 の差が 4 であることが確認できれば、「いとこ素数」のペアを一つ見つけたことになります。

ということで、時間のかかる素数判定は使わなくても「いとこ素数」は見つけられます。
ghci で確認してみると正しく動いているようです。素数判定を行なっていない分、オリジナルよりちょっぴり速くなっています。

*Main> take 10 cousinPrimes
[(3,7),(7,11),(13,17),(19,23),(37,41),(43,47),(67,71),(79,83),(97,101),(103,107)]
(0.00 secs, 514988 bytes)
*Main> take 10 cousins
[(3,7),(7,11),(13,17),(19,23),(37,41),(43,47),(67,71),(79,83),(97,101),(103,107)]
(0.00 secs, 0 bytes)
*Main> cousinPrimes !! 2000
(185869,185873)
(0.41 secs, 39850840 bytes)
*Main> cousins !! 2000
(185869,185873)
(1.51 secs, 137477412 bytes)

 

こちらのブログを読み進めていくと最後の方に「prime k-tuple」という言葉が出てきます。
これを見ていたら自作の cousinPrimes を改良してもっと一般化できそうな気がしてきました。

それで、素数の二つ組、三つ組、四つ組を求めるコードを作ってみました。引数には一つ目の素数との差を渡すことにしました。

-- 素数の二つ組
--
-- 双子素数 : prime2Tuples 2
-- いとこ素数 : prime2Tuples 4
-- セクシー素数 : prime2Tuples 6
--
prime2Tuples :: Integral a => a -> [(a, a)]
prime2Tuples d1 = loop primes
  where
    loop (p1 : ps1)
      | p2 - p1 == d1 = (p1, p2) : loop ps1
      | otherwise     = loop ps1
      where (p2 : _) = dropWhile (< p1 + d1) ps1


-- 素数の三つ組
--
-- 三つ子素数(1) : prime3Tuples 2 6
-- 三つ子素数(2) : prime3Tuples 4 6
-- セクシー素数の三つ組 : 
--   filter (\(p, _, _) -> (not . isPrime) (p + 18)) $ prime3Tuples 6 12
--
prime3Tuples :: Integral a => a -> a -> [(a, a, a)]
prime3Tuples d1 d2 = loop primes
  where
    loop (p1 : ps1)
      | checkDiff = (p1, p2, p3) : loop ps1
      | otherwise = loop ps1
      where
        (p2 : ps2) = dropWhile (< p1 + d1) ps1
        (p3 : _  ) = dropWhile (< p1 + d2) ps2
        checkDiff = map (subtract p1) [p2, p3] == [d1, d2]


-- 素数の四つ組
--
-- 四つ子素数 : prime4Tuples 2 6 8
--
prime4Tuples :: Integral a => a -> a -> a -> [(a, a, a, a)]
prime4Tuples d1 d2 d3 = loop primes
  where
    loop (p1 : ps1)
      | checkDiff = (p1, p2, p3, p4) : loop ps1
      | otherwise = loop ps1
      where
        (p2 : ps2) = dropWhile (< p1 + d1) ps1
        (p3 : ps3) = dropWhile (< p1 + d2) ps2
        (p4 : _  ) = dropWhile (< p1 + d3) ps3
        checkDiff = map (subtract p1) [p2, p3, p4] == [d1, d2, d3]

と、ここまで作って、「それぞれの関数で同じような処理を繰り返してるな」と気付きました。それで、さらにこれらの関数を一つにまとめてみました。
可変長引数をリストで与えることによって、二つ組でも三つ組でもこの関数一つで求められます。

-- Prime k-tuple の一般化
--
-- 双子素数 : primeKTuples [2]
-- いとこ素数 : primeKTuples [4]
-- セクシー素数 : primeKTuples [6]
--
-- 三つ子素数(1) : primeKTuples [2, 6]
-- 三つ子素数(2) : primeKTuples [4, 6]
-- セクシー素数の三つ組 :
--     filter (\(p : _) -> (not . isPrime) (p + 18)) $ primeKTuples [6, 12]
--
-- 四つ子素数 : primeKTuples [2, 6, 8]
-- セクシー素数の四つ組 : primeKTuples [6, 12, 18]
--
-- 五つ子素数(1) : primeKTuples [2, 6, 8, 12]
-- 五つ子素数(2) : primeKTuples [4, 6, 10, 12]
--
-- 六つ子素数 : primeKTuples [4, 6, 10, 12, 16]
--
primeKTuples :: Integral a => [a] -> [[a]]
primeKTuples ds = loop primes
  where
    loop (p : ps)
      | checkDiff = (p : ps') : loop ps
      | otherwise = loop ps
      where
        checkDiff = map (subtract p) ps' == ds
        ps' = map (head . (\d -> dropWhile (< p + d) ps)) ds

さらに、「三つ子素数」や「五つ子素数」が二つに分れているのが気になったので、マージソートで使われる merge 関数を導入してみました。

-- merge を使用すると「三つ子素数」や「五つ子素数」を一つにまとめられる。
-- 三つ子素数 : merge (primeKTuples [2, 6]) (primeKTuples [4, 6])
-- 五つ子素数 : merge (primeKTuples [2, 6, 8, 12]) (primeKTuples [4, 6, 10, 12])
merge :: Ord a => [a] -> [a] -> [a]
merge xs [] = xs
merge [] ys = ys
merge xs@(x : xs') ys@(y : ys')
  | x <= y = x : merge xs' ys
  | otherwise = y : merge xs ys'

例えば「三つ子素数」ならこんな感じで、まとめて小さい順に表示できます。

*Main> take 10 $ merge (primeKTuples [2, 6]) (primeKTuples [4, 6])
[[5,7,11],[7,11,13],[11,13,17],[13,17,19],[17,19,23],[37,41,43],[41,43,47],[67,71,73],[97,101,103],[101,103,107]]

HTML generated by org-mode 6.33x in emacs 23

2012年8月 7日 (火)

数当て問題 - Hit & Blow -

またしても人様のブログからネタを拝借。

こちらのブログに Haskell で書かれた「数当てゲーム」のコードがありました。
ルールは

  • コンピュータが生成した、重複のない 4 桁の数字を当てる。
  • コンピュータは、予想した数に対して、以下のヒントを返す。
    • Hit: 数字はあっているが、場所が異なる数字の数
    • Blow:数字も場所もあっている数字の数

というものです。

私も自分なりに考えてみました。

コンソールから数字を入力して、その度に Hit と Blow の数を返すようにしたのですが、入力が「文字列」になるので、Hit も Blow も文字列同士を比較することにしました。それによってコードが簡単になったと思います。

入力された文字列に対するエラーチェックはしていない、かなり手抜きのコードですが…。

import Data.List
import Random

countTrue :: [Bool] -> Int
countTrue = length . filter (== True)

countHit :: String -> String -> Int
countHit target = countTrue . map (`elem` target)

countBlow :: String -> String -> Int
countBlow target = countTrue . zipWith (==) target

checkNum :: String -> IO ()
checkNum target = do
  putStrLn "input number : "  
  input <- getLine
  let b = countBlow target input
  if b == 4
    then print "good !!"
    else do let h = countHit target input
            putStrLn $ show h ++ " hit(s), " ++ show b ++ "blow(s)\n"
            checkNum target

-- xs からランダムに重複しない n 個の要素を取り出す
randomPicUp :: Eq a => Int -> [a] -> IO [a]
randomPicUp n xs = loop n (length xs - 1) xs []
  where
    loop 0 _ _ prd = return prd
    loop c m xs prd = do
      g <- newStdGen
      let (i, _) = randomR (0, m) g
      let x = xs !! i
      loop (c - 1) (m - 1) (delete x xs) (x : prd)

main :: IO ()
main = do
  target <- randomPicUp 4 ['0'..'9']
  print target                  -- カンニング
  checkNum target

HTML generated by org-mode 6.33x in emacs 23

2012年4月 5日 (木)

オートマトン

ネットをウロウロしてたら、こちらのブログを発見しました。

「オートマトン」って、詳しくは知らないけど何だか面白そうです。 参考に挙げてあったPDFファイルを見ながら、自分なりにコードを書いてみました。

{-
   << DFA(決定性有限オートマトン) >>

   参考:
   http://d.hatena.ne.jp/cranebird/20120327/1332857936
   http://kurt.scitec.kobe-u.ac.jp/~kikyo/lec/07/automaton/k2.pdf
-}

type Q     = String
type Delta = ((Q, Char), Q)

-- isAccept 遷移表 受理状態の集合 状態 文字列 => 真偽
isAccept :: [Delta] -> [Q] -> Q -> String -> Bool
isAccept qs f st []        = elem st f
isAccept qs f st (c : cs)  = isAccept qs f st' cs
  where Just st' = lookup (st, c) qs -- 遷移表に抜けがないことが前提

-- makeDfa 遷移表 受理状態の集合 初期状態 => dfa
makeDfa :: [Delta] -> [Q] -> Q -> String -> String
makeDfa qs f q0 cs = if (isAccept qs f q0 cs)
                     then "Accept"
                     else "Not Accept"


{-
   例 1
   Σ = {a}, 奇数個の 'a' からなる文字列を受理する
-}
dfa1 :: String -> String
dfa1 = makeDfa qs f "Even"
  where
    qs = [(("Even", 'a'), "Odd"), (("Odd", 'a'), "Even")]
    f  = ["Odd"]

{-
   例 2
   Σ = {a, b}, 奇数個の 'a' からなる文字列を受理する
-}
dfa2 :: String -> String
dfa2 = makeDfa qs f "Even"
  where
    qs = [(("Even", 'a'), "Odd"), (("Even", 'b'), "Even"),
          (("Odd", 'a'), "Even"), (("Odd", 'b'), "Odd")]
    f  = ["Odd"]

{-
   例 3
   Σ = {a, b}, 奇数個の a と偶数個の b からなる文字列を受理する
   (0 は偶数とする)
-}
dfa3 :: String -> String
dfa3 = makeDfa qs f "00"
  where
    qs = [(("00", 'a'), "10"), (("00", 'b'), "01"),
          (("10", 'a'), "00"), (("10", 'b'), "11"),
          (("01", 'a'), "11"), (("01", 'b'), "00"),
          (("11", 'a'), "01"), (("11", 'b'), "10")]
    f  = ["10"]

{-
   問 4
   Σ = {a, b, c}, 奇数個の a と偶数個の b からなる文字列を受理する
   (0 は偶数とする)
-}
dfa4 :: String -> String
dfa4 = makeDfa qs f "00"
  where
    qs = [(("00", 'a'), "10"), (("00", 'b'), "01"), (("00", 'c'), "00"),
          (("10", 'a'), "00"), (("10", 'b'), "11"), (("10", 'c'), "10"),
          (("01", 'a'), "11"), (("01", 'b'), "00"), (("01", 'c'), "01"),
          (("11", 'a'), "01"), (("11", 'b'), "10"), (("11", 'c'), "11")]
    f  = ["10"]

 

とりあえず PDFファイル に書いてあった「例 1」〜「例 3」と「問 4」に対応する dfa を定義してみました。
ghci 上で確認してみたら、思ったように動いてくれているようです。

*Main> dfa1 "aaa"
"Accept"
*Main> dfa1 "aaaa"
"Not Accept"
*Main> dfa2 "ababa"
"Accept"
*Main> dfa2 "abababa"
"Not Accept"
*Main> dfa3 "ababa"
"Accept"
*Main> dfa3 "abababa"
"Not Accept"
*Main> dfa4 "abcacba"
"Accept"
*Main> dfa4 "abcacbacba"
"Not Accept"

 

でも PDFファイル の最後に書いてある「問 5」が未だにわからないんですよね。
御本家に聞いてみようかな?

 

 

HTML generated by org-mode 6.33x in emacs 23

2012年2月24日 (金)

リストモナド

こちらのブログ経由でこの記事にたどり着きました。
一つめのコードは「指定した文字列の大文字小文字の組み合わせ全パターンの リストを取得する」というコードなんですが…確かに少し冗長な気がします。
私なら次のようなコードを書くかな…。

import Data.Char

getAllUpperLowerPattern :: String -> [String]
getAllUpperLowerPattern cs = foldr f [""] cs
  where f c css = concat [[toUpper c : cs, toLower c : cs] | cs <- css]

 

と思った後にこちらのブログをよく読んだら、リストモナドに「sequence」と いう関数があったんですね。
この関数を使えばコードがさらに簡潔になります。

import Data.Char

getAllUpperLowerPattern :: String -> [String]
getAllUpperLowerPattern cs = sequence [[toUpper c, toLower c] | c <- cs]

 

やっぱり、もっとモナドの勉強をしなければいけませんね…。

 

 

Date: 2012-02-24 12:48:56 JST

HTML generated by org-mode 6.33x in emacs 23

2012年2月16日 (木)

Codeforces 35

最近、「Project Euler」を休んで、「Codeforces」の過去問を Haskell で解いて遊んでます。
ずっと「モナド」とか「IO アクション」から逃げてたんですけど、「Codeforces」では「入力」と「出力」が決まってるので「IO」の練習になるかな…と思いまして。

今回、この問題を解いてみました。
一般的な手続型言語では配列にデータを入れておいて、後はスワップを繰り返せばいいのでしょうけど、Haskell だとリストや配列のスワップはコストが高いよな…と思って、次のようなコードを書きました。

test :: String -> [String] -> Int
test "1" xs = loop 1 0 0 xs
test "2" xs = loop 0 1 0 xs
test "3" xs = loop 0 0 1 xs

loop :: Int -> Int -> Int -> [String] -> Int
loop 1 _ _ [] = 1
loop _ 1 _ [] = 2
loop _ _ 1 [] = 3
loop a b c (x1 : x2 : xs)
  | (x1, x2) == ("1", "2") || (x1, x2) == ("2", "1") = loop b a c xs
  | (x1, x2) == ("1", "3") || (x1, x2) == ("3", "1") = loop c b a xs
  | (x1, x2) == ("2", "3") || (x1, x2) == ("3", "2") = loop a c b xs

main :: IO ()
main = do input <- readFile "input.txt"
          let ([i], xs) = splitAt 1 $ words input
          writeFile "output.txt" $ show $ test i xs

一応、結果は「Accepted」でしたが、力技で解いた感があってあまり美しいと言えません。

 

今後の参考にと、他の人の解放を見ていたら、次のようなコードを見つけました。私とはちょっと発想が違っています。

import IO

gao :: [Int] -> Int
gao [a] = a
gao (a:b:c:d) = gao $ (if a == b then c else if a == c then b else a):d

main = do
        hi <- openFile "input.txt" ReadMode
        ho <- openFile "output.txt" WriteMode
        hGetContents hi >>= hPutStr ho . show . gao . map read . words
        hClose ho

このコードを参考に(パクッて?)、自分の好きなスタイルでコードを書いてみました。 (私は "if 〜 then 〜 else 〜 " があまり好きではありません。特に参考にしたコードのように if が連なっていると、直感的に分りにくい気がします)

test :: String -> [String] -> String
test i (a : b : xs)
  | i == a    = test b xs
  | i == b    = test a xs
  | otherwise = test i xs
test i [] = i

main :: IO ()
main = do input <- readFile "input.txt"
          let (i : xs) = words input
          writeFile "output.txt" $ test i xs

最初に自分で書いたコードより、こちらのコードのほうがずっと良い気がします。

Date: 2012-02-16 11:48:22 JST

HTML generated by org-mode 6.33x in emacs 23

2012年1月 4日 (水)

n Queen 問題

 先日、Haskell で「N-Queens 問題」を扱ったブログを見つけました。
 そういえば「N-Queens 問題」を自分で解いたことがなかったなぁ…ということで、自分でもコードを書いてみました。

import System
import Data.List

queen :: Int -> [[Int]]
queen n = loop n [([], [1 .. n])]
    where
      loop 0 xs = [x | (x, _) <- xs]
      loop n xs = loop (n - 1) $ concatMap addLine xs

-- 斜めに重複しない行を付け加える
addLine :: ([Int], [Int]) -> [([Int], [Int])]
addLine (xs, ys) = [(y : xs, delete y ys) | y <- ys, check y 1 xs]

-- 斜めの重複の判定する
check :: Int -> Int -> [Int] -> Bool
check y d (x : xs)
    | y - d == x || y + d == x = False
    | otherwise = check y (d + 1) xs
check _ _ _ = True


main :: IO ()
main = do n <- getArgs
          print $ length $ queen (read $ head n)

 クイーンの配置は「その行の左から何列目にクイーンがあるか」を整数で表わすことにしました。
 例えば

 ---------------
| Q |   |   |   |
|---+---+---+---|
|   | Q |   |   |
|---+---+---+---|
|   |   | Q |   |
|---+---+---+---|
|   |   |   | Q |
 ---------------

という配置は [1,2,3,4] というリストで表わされます。
 このようなデータ構造を使うと、n × n の盤面の場合 [1 .. n] の順列を考えれば縦と横の重複を調べる必要がありません。

 こちらこちらのコードと比べて処理時間が速かったので、自分では満足してたのですが、「ビット演算を使うともう少し速くなる」といったブログの記事を散見したので、時間があれば改良してみようかと思っています。

Date: 2012-01-04 11:45:13 JST

HTML generated by org-mode 6.33x in emacs 23

2011年12月15日 (木)

Haskell の配列

 こちらのブログで Haskell の二次元配列の話題が出ていました。
 コードの例が示してあり、「やはりCなどに比べると配列の準備が面倒ですね」と書いてあります。

 私は Haskell の配列に詳しいわけではないのですが、「listArray 関数」を使えば簡単に配列が作れるはずです。
 The Haskell 98 Reportには

配列は関数 listArray を使って境界対とインデックス順に ならべた値のリストから構築することもできる。
と書いてあります。

 

Prelude> :m +Data.Array Prelude Data.Array> listArray (0, 3) [0 .. 3] array (0,3) [(0,0),(1,1),(2,2),(3,3)] Prelude Data.Array> listArray ((0, 0), (1, 1)) [1 .. 4] array ((0,0),(1,1)) [((0,0),1),((0,1),2),((1,0),3),((1,1),4)]
 このように listArray 関数を使えば多次元配列も簡単に作れます。

 「listArray 関数」を使って作った配列には、(!) でインデックスを指定することでアクセスできます(多次元配列にも一発でアクセスできます)。

Prelude Data.Array> let a = listArray ((0, 0), (1, 1)) [1 .. 4] Prelude Data.Array> a array ((0,0),(1,1)) [((0,0),1),((0,1),2),((1,0),3),((1,1),4)] Prelude Data.Array> a ! (0,1) 2

 

 ということで、こちらのブログのコードと同様なコードを「listArray 関数」を使って書いてみました。
 "(!)" で直接アクセスできるので、「seeValue 関数」は定義しませんでした(「makeArray 関数」も定義しなくてよさそうでしたが、これは一応残しました)。

import Data.Array makeArray :: (Int, Int) -> [a] -> Array (Int, Int) a makeArray (x, y) xs = listArray ((1, 1), (x, y)) xs makeRange :: String -> (Int, Int) makeRange cs = (x, y) where [x, y] = map read $ take 2 $ words cs toIntList :: String -> [Int] toIntList cs = map read $ words cs main :: IO () main = do xy <- getLine xs <- getContents let arr = makeArray (makeRange xy) (toIntList xs) print arr print $ arr ! (2, 3)
$ ./array 4 4 1 2 3 4 5 6 7 8 1 2 4 8 7 7 7 7 array ((1,1),(4,4)) [((1,1),1),((1,2),2),((1,3),3),((1,4),4),((2,1),5),((2,2),6),((2,3),7),((2,4),8),((3,1),1),((3,2),2),((3,3),4),((3,4),8),((4,1),7),((4,2),7),((4,3),7),((4,4),7)] 7

2011年12月 1日 (木)

またまた「FizzBuzz」

 かなり久しぶりのブログ更新です。
 最近、また「FizzBuzz問題」が気になりまして……。

 

 一般的な「FizzBuzz問題」といえばこちらのルールですよね。
 Haskell だとこんな感じでしょうか?

fizzbuzz :: [String] fizzbuzz = map f [1 .. ] where f x | rem x 15 == 0 = "FizzBuzz" | rem x 3 == 0 = "Fizz" | rem x 5 == 0 = "Buzz" | otherwise = show x main :: IO () main = putStr $ unlines $ take 100 fizzbuzz

 

 そういえば、「剰余関数を使わない」という制限を加えたものもありましたね。ネット上では次のようなコードが有名(?)だったと思います。

fizzbuzz :: [String] fizzbuzz = zipWith f [1 .. ] $ zipWith (++) fizz buzz where f x "" = show x f _ y = y fizz = cycle ["", "", "Fizz"] buzz = cycle ["", "", "", "", "Buzz"] main :: IO () main = putStr $ unlines $ take 100 fizzbuzz

 

 この方法とは別に、リストの要素に n 個置きに関数を適用する step 関数を Haskell で定義すると、「剰余関数を使わない」で FizzBuzz ができます。

-- リストの先頭から n 個目ごとに f を適用 -- ex : step 3 (\x -> 0) [1 .. 10] => [0,2,3,0,5,6,0,8,9,0] step :: Int -> (a -> a) -> [a] -> [a] step _ _ [] = [] step n f xs = f a : as ++ step n f bs where (a : as, bs) = splitAt n xs fizzbuzz :: [String] fizzbuzz = tail $ foldl conv ss factors where ss = map show [0 ..] factors = [(3, "Fizz"), (5, "Buzz"), (15, "FizzBuzz")] conv ss (a, b) = step a (\ x -> b) ss main :: IO () main = putStrLn $ unlines $ take 100 fizzbuzz

 初めのうちは、ss = map show [0 .. 100] としていたのですが、試しに無限数列にしてみたら、ちゃんと動きました。
 なぜちゃんと動くかが自分でも分っていないのですが、遅延評価ってすごいなぁと思います。

 

 さて、ここからが本題なのですが、「3 と 5 意外の数で FizzBuzz してもいいじゃないか」ということで、任意の数の倍数を任意の文字列で置き換える「拡張版 FizzBuzz」を考えてみました。

{- 任意の数と文字列を指定できる FizzBuzz -} import System.Environment (getArgs) -- ex : fizzbuzz [(2, "Fizz"), (3, "Buzz")] -- -> ["1","Fizz","Buzz","Fizz","5","FizzBuzz","7","Fizz" ..] fizzbuzz :: [(Int, String)] -> [String] fizzbuzz xs = map f [1 ..] where f n = if null str then show n else str where str = foldr g [] xs g (x, cs) tmp = if rem n x == 0 then cs ++ tmp else tmp -- ex : fizzbuzz 100 3 fizz 5 buzz main :: IO () main = do args <- getArgs let n = read (head args) xs = format (tail args) in mapM_ putStrLn $ take n $ fizzbuzz xs where format [] = [] format (n : x : xs) = (read n, x) : format xs

 実行結果は次のようになります。

$ runghc fizzbuzz8.hs 50 2 Fizz 3 Buzz 5 Bar 1 Fizz Buzz Fizz Bar FizzBuzz 7 Fizz Buzz FizzBar 11 FizzBuzz 13 Fizz BuzzBar Fizz 17 FizzBuzz 19 FizzBar Buzz Fizz 23 FizzBuzz Bar Fizz Buzz Fizz 29 FizzBuzzBar 31 Fizz Buzz Fizz Bar FizzBuzz 37 Fizz Buzz FizzBar 41 FizzBuzz 43 Fizz BuzzBar Fizz 47 FizzBuzz 49 FizzBar

より以前の記事一覧

2016年7月
          1 2
3 4 5 6 7 8 9
10 11 12 13 14 15 16
17 18 19 20 21 22 23
24 25 26 27 28 29 30
31            
フォト

最近のトラックバック

無料ブログはココログ