« Project Euler : Problem 58 | トップページ | 基数ソート »

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]

« Project Euler : Problem 58 | トップページ | 基数ソート »

Haskell」カテゴリの記事

数学」カテゴリの記事

コメント

コメントを書く

コメントは記事投稿者が公開するまで表示されません。

(ウェブ上には掲載しません)

トラックバック

この記事のトラックバックURL:
http://app.cocolog-nifty.com/t/trackback/112020/50859380

この記事へのトラックバック一覧です: Haskell で行列の計算をしてみる:

« Project Euler : Problem 58 | トップページ | 基数ソート »

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

最近のトラックバック

無料ブログはココログ