数学

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

2011年2月13日 (日)

Haskell で行列の計算をしてみる

 みなさんは『オイラーの贈物』という本を御存じですか?
 「オイラーの公式」を導くために必要な、代数、幾何、解析の知識がこの一冊で学べるという本です。高校から大学で習う数学の知識があれば、素晴しい数学の世界を垣間見ることができます。いろいろな定義や公式が導き方とともに載っているので、数学的な知識の確認をしたい時に重宝します。

 ここ数日、この本の「ベクトルと行列」の章を読み返していました。
 読んでいて、ふと「Haskell で行列の計算をするためのモジュールを作ってみよう」と思ました。
 Haskell の行列演算モジュールなど、絶対誰かが作っているはずですが、ちょっと頭の体操代わりに自分で作ってみたくなったんです。
 で、できたのがこれです。

{- Matrix.hs -- ベクトル及び行列の計算 coding by Tsumuji -} module Matrix where import Data.List (transpose) -- ベクトルの演算 -- * ベクトルは数値のリストの形で表す。 -- ex : 3次元のベクトル [1,2,3] type Vector a = [a] -- Scalar * Vector (*/) :: Num a => a -> Vector a -> Vector a (*/) n xs = map (* n) xs infixl 7 */ -- Vector + Vector (/+/) :: Num a => Vector a -> Vector a -> Vector a (/+/) xs ys = zipWith (+) xs ys infixl 6 /+/ -- Vector - Vector (/-/) :: Num a => Vector a -> Vector a -> Vector a (/-/) xs ys = zipWith (-) xs ys infixl 6 /-/ -- Vector * Vector (inner product) (/*/) :: Num a => Vector a -> Vector a -> a (/*/) xs ys = sum $ zipWith (*) xs ys infixl 7 /*/ -- Vector の絶対値 vectorAbs :: Floating a => Vector a -> a vectorAbs xs = sqrt (xs /*/ xs) -- n 次元のゼロベクトル zeroVector :: Num a => Int -> Vector a zeroVector n = take n $ repeat 0 -- 行列の演算 -- * 行列は行ベクトルのリストの形で表す。 -- ex : 2行3列の行列 [[1,2,3],[4,5,6]] -- * 行と列の整合性は調べていないので注意!! type Matrix a = [Vector a] -- 行列の m 行目 n 列目の要素(通常のリストと同様に 0 から数え始める) matrixRef :: Num a => Int -> Int -> Matrix a -> a matrixRef m n xss = (xss !! m) !! n -- Scalar * Matrix (*|) :: Num a => a -> Matrix a -> Matrix a (*|) n xss = map (n */) xss infixl 7 *| -- Matrix + Matrix (|+|) :: Num a => Matrix a -> Matrix a -> Matrix a (|+|) xss yss = zipWith (/+/) xss yss infixl 6 |+| -- Matrix - Matrix (|-|) :: Num a => Matrix a -> Matrix a -> Matrix a (|-|) xss yss = zipWith (/-/) xss yss infixl 6 |-| -- Matrix * Matrix (|*|) :: Num a => Matrix a -> Matrix a -> Matrix a (|*|) xss yss = [[xs /*/ ys | ys <- yss'] | xs <- xss] where yss' = transpose yss infixl 7 |*| -- m 行 n 列のゼロ行列 zeroMatrix :: Num a => Int -> Int -> Matrix a zeroMatrix m n = take m $ repeat $ zeroVector n -- n 行 n 列の単位行列 unitMatrix :: Num a => Int -> Matrix a unitMatrix n = [take n $ f x | x <- [0 .. n - 1]] where f x = replicate x 0 ++ [1] ++ repeat 0 -- Matrix ^ n (正方行列のみ可) (|^) :: Num a => Matrix a -> Int -> Matrix a (|^) xss 0 = unitMatrix (length xss) (|^) xss n = power (|*|) xss n infixr 8 |^ -- 高速累乗計算 -- ex : 2^100 == power (*) 2 100 power :: (t -> t -> t) -> t -> Int -> t power op a n = loop (n - 1) a a where loop 0 _ r = r loop n x r = loop (div n 2) (op x x) r' where r' = if odd n then op r x else r
 ベクトルや行列の演算では、要素にランダムにアクセスすることは無いので、データの型をベクトルは「数字のリスト」、行列は「ベクトルのリスト」としました。
 ベクトルや行列のサイズや整合性をチェックしていない手抜きバージョンなので、使うときにはちょっと注意が必要ですが……。

 さらに、『オイラーの贈物』には行列の計算でフィボナッチ数列を求める方法も載っていたので、これもコーディングしてみました。

{- fib_matrix.hs coding by Tsumuji -} import Matrix -- フィボナッチ数列の第n項を求める (n > 0) fibonacci :: Int -> Integer fibonacci n = matrixRef 1 0 $ [[1,1],[1,0]] |^ (n - 1) |*| [[1],[1]] main :: IO () main = print $ map fibonacci [1..20]
>ghc -O -o fib fib_matrix.hs --make ghc -O -o fib fib_matrix.hs --make Linking fib.exe ... >fib fib [1,1,2,3,5,8,13,21,34,55,89,144,233,377,610,987,1597,2584,4181,6765]

2011年1月29日 (土)

平方根の求め方 その2 〜 今度は Haskell で

 以前、Ruby で平方根を求めたことがありましたが、今度は Haskell でやってみました。
 今回は入力(引数)も出力(戻値)も桁数の制限を受けないように、どちらの型も String にしてみました。
 また、今回は「筆算による開平法」をシュミレートしたものも作ってみました。

 

《 筆算による開平法 》

 「筆算による開平法」の具体な方法については、こちらを参照してください。

import Data.Char {- -- 筆算による開平法 -- -} squareRoot :: String -> String squareRoot cs = root' (int' ++ frac' ++ repeat '0') 0 0 where (int, frac) = break (== '.') cs int' = if odd (length int) then ('0' : int) else int frac' = if null frac then "." else frac root' ('.' : xs) r d = '.' : root' xs r d root' (a : b : xs) r d = (intToDigit . fromIntegral) n' : root' xs r' d' where n = r * 100 + read [a, b] n' = last $ takeWhile (\x -> (d + x) * x <= n) [0 .. 9] r' = n - (d + n') * n' d' = 10 * (d + n' + n')

 このコードは、まず始めに整数部分と少数部分を整形しています。
 整数部分は偶数桁になるように、必要があれば最上位の桁に '0' を補います。
 少数部分は引数に少数点が含まれていなければ補い、さらに '0' の無限リストを付け加えています。
 その後は、root' 関数で「筆算による開平法」をシミュレートしているだけです。

 実行すると、次のようになります。

*Main> take 50 $ squareRoot "2" "1.414213562373095048801688724209698078569671875376" *Main> take 50 $ squareRoot $ show 3.0 "1.732050807568877293527446341505872366942805253810" *Main> take 50 $ squareRoot "0.000000000000000000000000123" "0.000000000000350713558335003638336349349661310276"
 squareRoot 関数をそのまま実行すると、メモりが枯渇するまで文字列を出力し続けるので、take などを使って適当に切ってください。(その場合、少数点も一桁分に数えられるのでご注意を)
 また、大きな桁の数を show を使って文字列に変換すると、浮動少数点表記になってしまい、そのまま squareRoot 関数に渡すとエラーになります。
*Main> show 0.000000000000000000000000123 "1.23e-25" *Main> take 50 $ squareRoot $ show 0.000000000000000000000000123 "1.1*** Exception: Prelude.read: no parse"

 

《 タイガー計算機による開平法 》

 「タイガー計算機による開平法」の具体的な方法については、こちらを参照してください。

import Data.Char {- -- タイガー計算機による開平法 -- -} squareRoot2 :: String -> String squareRoot2 cs = root' 0 (int' ++ frac' ++ repeat '0') 1 where (int, frac) = break (== '.') cs int' = if odd (length int) then ('0' : int) else int frac' = if null frac then "." else frac root' rest ('.' : xs) n = '.' : root' rest xs n root' rest (a : b : xs) n = intToDigit c : root' rest' xs (n' * 10 - 9) where (rest', n', c) = count (rest * 100 + read [a, b]) n -- count : -- m から 奇数 n, n + 2 .. を順に引いていった時、 -- (残った数, 次に引くはずだった奇数, 奇数を引いた回数) を返す。 count m n = head $ dropWhile (\(x, y, _) -> x >= y) $ iterate g (m, n, 0) where g (a, b, c) = (a - b, b + 2, c + 1)

 このコードの内容は root' 関数の中身が違う以外は「筆算による開平法」の squareRoot と同じです。ですから、注意事項も同じです。
 ただし、squareRoot が内部でかけ算を多用しているのに対し、こちらは内部で引き算を多用しているせいか、若干こちらのほうが速い気がします。

 

 ネットを見ていたら、私と同じように「筆算による開平法」を Haskell でコーディングしている方がいらっしゃいました。でもやっていることは同じはずなのに、処理の仕方が違うのでコードも全然違って見えます。
 こんな風にコーディングに個性が出るのが、プログラミングの面白いところだと思います。

2010年2月20日 (土)

Haskell で「順列」と「組み合わせ」

 "Project Euler" などの数学的な問題を解くときに「辞書順に生成された順列」を使いたい時があります。( ここでいう「辞書順」とは、例えば "[1,2,3]" を関数に渡したときに "[[1,2,3],[1,3,2],[2,1,3],[2,3,1],[3,1,2],[3,2,1]]" が返ってくることを言います )
 Hasekell には Data.List モジュールに "permutations" という順列を生成する関数があるのですが、これが辞書順の答えを返してくれないんです。"[1,2,3]" を渡すと "[[1,2,3],[2,1,3],[3,2,1],[2,3,1],[3,1,2],[1,3,2]]" が返ってきます。

 そこで、Haskell の練習も兼ねて自分で作ってみることにしました。
 紙と鉛筆で順列を考える場合のように、樹形図を作っていくイメージです。

permutation :: [a] -> Int -> [[a]] permutation [] _ = [[]] permutation ns size = perm size (length ns - 1) [(ns, [])] where perm 0 _ xs = [a | (_, a) <- xs] perm c n xs = perm (c - 1) (n - 1) $ concatMap (f n) xs f n (xs, ys) = [(as ++ bs, ys ++ [b]) | x <- [0 .. n], let (as, b : bs) = splitAt x xs]
*Main> permutation [1 .. 4] 2 [[1,2],[1,3],[1,4],[2,1],[2,3],[2,4],[3,1],[3,2],[3,4],[4,1],[4,2],[4,3]]

 ちょっと行数が多い気はしますが、ネットで見つけたいくつかの「辞書順に順列を生成する関数」と比べても 4 〜 5 倍の速さで答えを出せるようなので自分としては満足かなと……。

 

 さらにこの関数をもとに「組み合わせを生成する関数」もつくってみました。

combination :: [a] -> Int -> [[a]] combination [] _ = [[]] combination ns size = comb size [(ns, [])] where comb 0 xs = [a | (_, a) <- xs] comb c xs = comb (c - 1) $ concatMap comb' xs comb' (x : xs, ys) = (xs, ys ++ [x]) : comb' (xs, ys) comb' _ = []
*Main> combination [1 .. 4] 2 [[1,2],[1,3],[1,4],[2,3],[2,4],[3,4]]

2009年9月 3日 (木)

平方根の求め方

みずぴーさんのブログを見ていたら、面白い記事を見つけました。

"奇数を順に引いていくだけで平方根が求められる" というのです。


リンクをたどって行って、『平方根の求め方』というホームページにたどり着きました。

なんでも、歯車式の計算機で平方根を求める方法らしいです。

ホームページの説明をよく見てみると……平方数の考え方を使って平方根を求めているのですね。こんな方法を考えつくなんてすごいです。(ホームページの説明は正方形の面積を使って平方数をイメージできるようになっていて、非常にわかりやすかったです)


面白そうだったので、早速 Ruby でコードを書いてみました。

# -*- coding: utf-8 -*- # # = root.rb # # from 『平方根の求めかた』 (http://www.wizforest.com/gear/tiger/sqrt/) # class Numeric # call-seq: # n.root -> float # # * n の平方根を返す。 # # === Example # 2.root #=> 1.4142135623731 # def root n = self d1 = 0 # 整数部が 2 桁以下になるように繰り下げる。 n, d1 = n / 100.0, d1 - 1 while n >= 100 # 整数部が 1 桁以上になるように繰り上げる。 n, d1 = n * 100.0, d1 + 1 while n < 1 ans = 0 odd = 1 # 引いていく奇数 0.upto(1/0.0) do |d2| # d2 : ans の桁数 0.upto(1/0.0) do |c| if n < odd ans = ans * 10 + c break end n = n - odd odd = odd + 2 end if (d2 > 16 || n == 0) return ans / 10.0 ** (d1 + d2) end n = n * 100 odd = (odd - 1) * 10 + 1 end end end if __FILE__ == $0 n = ARGV[0].to_f puts n.root end

他のメソッドからも呼び出せるように最後の答えを float 型で返すようにしたので、有効桁数に制限がありますが、"Math#sqrt" や "** 0.5" の結果と同じ値が出せます。(ただし、float 型なので計算上の微妙な誤差が出るようで、"==" で評価すると "false" が帰ってくることはあるようです)

コードの最初の部分で桁数を調整してから計算を始めているので、かなり大きな数からかなり小さな数までしっかり計算できます。

$ ruby root.rb 2 1.4142135623731 $ ruby root.rb 3 1.73205080756888 $ ruby root.rb 12345678901234567890 3513641828.82014 $ ruby root.rb 0.00000000000000000000000123 1.10905365064094e-12

今回は float 型で値を返すようにしたため有効桁数に制限がありますが、"BigDecimal" を使えば有効桁数の制限を気にせずに計算できると思います。まあ、これは今後の課題ということで……。

 

追記 (2010/05/04) : 以前の "root.rb" が "while" を多用していて、あまり Ruby っぽくないように思えたので、"Integer#upto" を使ったコードの変えてみました。

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            
フォト

最近のトラックバック

無料ブログはココログ