ガウスの消去法

Haskell

2chでガウスの消去法をHaskellでどう書くかという話題があったので、書いてみた:

module Gauss where

type E a = (a, [a])

forward :: Fractional a => [E a] -> [E a]
forward [] = []
forward (h@(a, k:kt):t) = h : forward (map f t)
  where
    f (a', k':kt') =
      let g x x' = x' - x * k' / k
          a'' = g a a'
          kt'' = zipWith g kt kt'
      in (a'', kt'')

back :: Fractional a => [E a] -> [a]
back [] = []
back ((a, k:[]):t) = x : back (map f t)
  where
    x = a / k
    f (b, k':t') = (b - x * k', t')

gauss :: Fractional a => [E a] -> [a]
gauss = reverse . back . reverse . map f . forward
  where
    f (a, ks) = (a, reverse ks)

動作確認として

  • 2x - 3y + z = 1
  • x + 2y - 3z = 4
  • 3x + 2y - z = 5

の連立方程式を解いてみる。

*Gauss> :m + Data.Ratio
*Gauss Data.Ratio> gauss [(1, [2,-3,1]), (4, [1,2,-3]),(5,[3,2,-1])]  :: [Ratio Integer]
[5 % 4,1 % 4,(-3) % 4]
*Gauss Data.Ratio>

x = 5/4, y = 1/4, z = -3/4の答えがちゃんと出ている。

事前のチェックをしていないので、不完全なところもある。たぶん、解が定まらないときにはゼロ除算例外を出すと思う。まぁ、それは構わない。そういうときに例外を出すのは、一つのありうる設計。

それはそれとして、reverseを3回、それも1回はmapの中で使っていることが気に入らない。そのおかげでnを変数の数として、O(n^2)になるはず…　だけど、よく考えたらガウスの消去法はもともとO(n^2)かな？　それなら、まぁ…　いやよくない。