На самом деле, помощь понадобилась четверокурсникам. Они уже давно забыли, как решать системы линейных уравнений, да и не заставишь их (в том числе меня). Но что поделаешь, если преподавателю по численным методам захотелось народ помучить? Что касается меня, меня решать что-то уже не заставишь, мне проще написать программу, которая делает это автоматически. Собственно, я и написал программу, которая решает линейную систему по формулам крамера, выводя все сделанные выкладки в TeX, при этом определитель находится не методом Гаусса, а "наивным" методом, как обычно и делают в институте. Собственно, это первое (размером более 50 строк), что я написал на Хаскелле. Прошу оценить и помочь советом, что делается по-другому. Собственно, сам код (компилится последним GHC с опциями -XMultiParamTypeClasses -XTypeSynonymInstances, причём вторую опцию поставил по интуиции, зачем она нужна, толком не понял):
import qualified Array
import Data.Ratio
import Data.Ix
-------------------
-- StringBuilder
--
infixl 0 :++
data StringBuilder = SB String | StringBuilder :++ StringBuilder
buildString :: StringBuilder -> String
buildString builder =
let
build' (SB content) acc = content ++ acc
build' (l :++ r) acc = build' l (build' r acc)
in
build' builder []
emptyBuilder = SB []
join :: StringBuilder -> [StringBuilder] -> StringBuilder
join sep [] = emptyBuilder
join sep (h:t) = h :++ foldr (\e -> \acc -> sep :++ e :++ acc) emptyBuilder t
-------------------
-- Type classes
infix 9 !
class (Ix i) => Indiced t i c where
(!) :: t -> i -> c
class ShowTex t where
showTex :: t -> StringBuilder
-------------------
-- Matrix
type MatrixIndex = (Int, Int)
data Matrix = Matrix (Array.Array MatrixIndex Rational)
instance Indiced Matrix MatrixIndex Rational where
(Matrix content) ! index = content Array.! index
mkMatrix :: Int -> [((Int, Int), Rational)] -> Matrix
mkMatrix size content = Matrix (Array.array ((1, 1), (size, size)) content)
(//) :: Matrix -> [((Int, Int), Rational)] -> Matrix
(Matrix content) // newContent = Matrix (content Array.// newContent)
matrixSize :: Matrix -> Int
matrixSize (Matrix content) = let ((_, _), (size, _)) = Array.bounds content in size
instance ShowTex Matrix where
showTex matrix =
let
size = matrixSize matrix
rowToTex :: Int -> StringBuilder
rowToTex row = join (SB "&") [showTex ((matrix!(row, i :: Int)) :: Rational) | i <- [1..size]]
in
SB "\\begin{array}{" :++ (SB . take size . repeat) 'c' :++ SB "}" :++
join (SB "\\\\") [rowToTex i | i <- [1..size]] :++
SB "\\end{array}"
instance ShowTex Rational where
showTex value =
let
denom = abs (denominator value)
numer = abs (numerator value)
absTex =
if denom == 1 then
SB (show numer)
else
SB "\\frac{" :++ SB (show numer) :++ SB "}{" :++ SB (show denom) :++ SB "}"
in
if value < 0 then SB "(-" :++ absTex :++ SB ")" else absTex
simpleMatrix :: Int -> [[Rational]] -> Matrix
simpleMatrix size content = mkMatrix size
(foldr (++) [] [extractRow row rowIndex | (row, rowIndex) <- zip content [1..size]])
where
extractRow :: [Rational] -> Int -> [((Int, Int), Rational)]
extractRow row rowIndex = foldr (:) []
[((rowIndex, colIndex), v) | (v, colIndex) <- zip row [1..size]]
minor :: (Int, Int) -> Matrix -> Matrix
minor (row, col) matrix =
let
size = matrixSize matrix
fixIndices r c = (if r >= row then r + 1 else r, if c >= col then c + 1 else c)
in
mkMatrix (size - 1) [((row, col), matrix!(fixIndices row col)) | row <- [1..size-1], col <- [1..size-1]]
-------------------
-- Vector
data Vector = Vector (Array.Array Int Rational)
mkVector :: Int -> [(Int, Rational)] -> Vector
mkVector size content = Vector (Array.array (1, size) content)
updateVector :: Vector -> [(Int, Rational)] -> Vector
updateVector (Vector content) newContent = Vector (content Array.// newContent)
vectorSize :: Vector -> Int
vectorSize (Vector content) = let (_, size) = Array.bounds content in size
instance Indiced Vector Int Rational where
(Vector content) ! index = content Array.! index
instance ShowTex Vector where
showTex vector =
let
size = vectorSize vector
in
SB "\\begin{array}{c}" :++
join (SB "\\\\") [showTex ((vector!(i)) :: Rational) | i <- [1..size]] :++
SB "\\end{array}"
-------------------
-- Utilities
splitString :: String -> [String]
splitString s =
let
split' (h:t) =
let
(p:q) = split' t
in
if h == ' ' || h == '\r' || h == '\n' then
if null p then p : q else [] : p : q
else
(h : p) : q
split' [] = [[]]
in
split' s
readIntegers :: String -> [Int]
readIntegers = map read . splitString
readExtMatrix :: String -> (Matrix, Vector)
readExtMatrix text =
let
(size : content) = readIntegers text
indices = [(r, c) | r <- [1..size], c <- [1..size]]
vectorContent = drop (size * size) content
in
(mkMatrix size (zip indices [(toInteger x) % 1 | x <- content]),
mkVector size (zip [1..size] [(toInteger x) % 1 | x <- vectorContent]))
makeVarVectorTex :: StringBuilder -> Int -> StringBuilder
makeVarVectorTex var size =
SB "\\left(\\begin{array}{c}" :++
join (SB "\\\\") [var :++ SB "_" :++ SB (show i) | i <- [1..size]] :++
SB "\\end{array}\\right)"
replaceColumn :: Int -> Vector -> Matrix -> Matrix
replaceColumn index v m = m // [((i, index), (v!i) :: Rational) | i <- [1..matrixSize m]]
-------------------
-- Solving
data Solution a = Solution StringBuilder a
determ :: Matrix -> Solution Rational
determ matrix =
let
size = matrixSize matrix
in
case size of
1 ->
Solution emptyBuilder (matrix!(1 :: Int, 1 :: Int))
2 ->
let
value =
matrix!(1 :: Int, 1 :: Int) * matrix!(2 :: Int, 2 :: Int) -
matrix!(2 :: Int, 1 :: Int) * matrix!(1 :: Int, 2 :: Int)
builder' =
SB "$$\\left|" :++ showTex matrix :++ SB "\\right| = " :++
showTex ((matrix!(1 :: Int, 1 :: Int)) :: Rational) :++ SB " \\cdot " :++
showTex ((matrix!(2 :: Int, 2 :: Int)) :: Rational) :++ SB " - " :++
showTex ((matrix!(1 :: Int, 2 :: Int)) :: Rational) :++ SB " \\cdot " :++
showTex ((matrix!(2 :: Int, 1 :: Int)) :: Rational) :++ SB " = " :++
showTex value :++ SB "$$\r"
in
Solution builder' value
n ->
let
addSolution index (solutions, builder) =
let sol@(Solution builder' _) = determ (minor (1, index) matrix)
in (sol : solutions, builder' :++ builder)
(solutions', builder) = foldr addSolution ([], emptyBuilder) [1..n]
solutions = reverse solutions'
addResult (Solution _ v, index) acc =
acc + (if index `mod` 2 == 0 then -1 else 1) * v * matrix!(1 :: Int, index :: Int)
result = foldr addResult 0 (zip solutions [1..n])
appendMinor index builder =
let
v = (matrix!(1 :: Int, index :: Int)) :: Rational
minorTex = showTex (minor(1, index) matrix)
signTex = if index == 1 then SB "" else if index `mod` 2 == 0 then SB "-" else SB "+"
in
signTex :++ showTex v :++ SB "\\cdot\\left|" :++ minorTex :++ SB "\\right|" :++ builder
appendMinor' (Solution _ m, index) builder =
let
v = (matrix!(1 :: Int, index :: Int)) :: Rational
minorTex = showTex m
signTex = if index == 1 then SB "" else if index `mod` 2 == 0 then SB " - " else SB " + "
in
signTex :++ showTex v :++ SB " \\cdot " :++ minorTex :++ builder
builder' =
builder :++
SB "$$\\left|" :++ showTex matrix :++ SB "\\right| = " :++
foldr appendMinor emptyBuilder [1..n] :++ SB " = " :++
foldr appendMinor' emptyBuilder (zip solutions' [1..n]) :++ SB " = "
in
Solution (builder' :++ showTex result :++ SB "$$\r") result
kramerSolve :: (Matrix, Vector) -> Solution Vector
kramerSolve (m, v) =
let
size = matrixSize m
task =
SB "$$\\left(" :++ showTex m :++ SB "\\right)" :++ makeVarVectorTex (SB "x") size :++
SB " = \\left(" :++ showTex v :++ SB "\\right)$$\r"
makeBindingTex var matrix det =
SB "$$" :++ var :++ SB " = \\left|" :++ showTex matrix :++
SB "\\right| = " :++ showTex det :++ SB "$$\r"
makeIndexBindingTex var index matrix det =
makeBindingTex (var :++ SB "_" :++ SB (show index)) matrix det
(Solution prologue det) = determ m
matrices = [replaceColumn i v m | i <- [1..size]]
solutions = [determ x | x <- matrices]
prologue' =
makeBindingTex (SB "\\Delta") m det :++
join (SB "") [s :++ makeIndexBindingTex (SB "\\Delta") i m v | ((Solution s v), m, i) <- zip3 solutions matrices [1..size]]
resultList = [d / det | (Solution _ d) <- solutions]
result = mkVector size (zip [1..size] resultList)
epilogue =
SB "$$" :++ makeVarVectorTex (SB "x") size :++
SB " = \\left(" :++ showTex result :++ SB "\\right)$$\r"
in
Solution (task :++ prologue :++ prologue' :++ epilogue) result
-------------------
-- Main
main = do
input <- getContents
m <- do return (readExtMatrix input)
putStrLn (let (Solution b _) = kramerSolve m in buildString b)
Меня одолевают смутные сомнения, что в случае со StringBuilder я изобрёл велосипед, но когда мне там разбираться в неудобной справке GHC? Кроме того, я уверен, что существуют готовые модули для работы с линейной алгеброй, однако, как мне кажется, от них мало толку, так как нужно не просто получить ответ, а расписать выкладки.
... << RSDN@Home 1.2.0 alpha rev. 672>>