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]

平方根の求め方 その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 でコーディングしている方がいらっしゃいました。でもやっていることは同じはずなのに、処理の仕方が違うのでコードも全然違って見えます。
　こんな風にコーディングに個性が出るのが、プログラミングの面白いところだと思います。

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]]

平方根の求め方

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

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

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

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

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

面白そうだったので、早速 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" を使ったコードの変えてみました。