Решил выложить код вычисления определителей. Код рабочий, хотя и не претендует на виртуозность. Просто было интересно решить эту задачу именно на Хаскелле. Рассмотрены два подхода к решению задачи: простая рекурсия и метод Гаусса.

Немного теории


Как известно, определитель квадратной матрицы n*n — это сумма n! слагаемых, каждое из которых есть произведение, содержащее ровно по одному элементу матрицы из каждого столбца и ровно по одному из каждой строки. Знак очередного произведения:

${a}_{1,i1}*{a}_{2,i2}*...{a}_{n,in}$


определяется чётностью подстановки:

$\begin{pmatrix}1 & 2 & ... & n \\ {i}_{1} & {i}_{2} & ... & {i}_{n} \end{pmatrix}$



Прямой метод вычисления определителя состоит в разложении его по элементам строки или столбца в сумму произведений элементов какой-либо строки или столбца на их алгебраические дополнения. В свою очередь, алгебраическое дополнение элемента матрицы

${a}_{i,j}$

есть

${(-1)}^{i+j}*{M}_{i,j}$

при этом

${M}_{i,j}$

— есть минор элемента (i,j), т.е. определитель, получающийся из исходного определителя вычеркиванием i-й строки и j-го столбца.

Такой метод порождает рекурсивный процесс, позволяющий вычислить любой определитель. Но производительность этого алгоритма оставляет желать лучшего — O(n!). Поэтому применяется такое прямое вычисление разве что при символьных выкладках (и с определителями не слишком высокого порядка).

Гораздо производительнее оказывается метод Гаусса. Его суть основывается на следующих положениях:

1. Определитель верхней треугольной матрицы \begin{pmatrix}{a}_{1,1} & {a}_{1,2} &… & {a}_{1,n} \\ 0 & {a}_{2,2} &… & {a}_{2,n} \\ 0 & 0 &… & ...\\ 0 & 0 &… & {a}_{n,n} \\\end{pmatrix} равен произведению ее диагональных элементов. Этот факт сразу же следует из разложения определителя по элементам первой строки или первого столбца.

2. Если в матрице к элементам одной строки прибавить элементы другой строки, умноженные на одно и то же число, то значение определителя не изменится.

3. Если в матрице поменять местами две строки (или два столбца), то значение определителя изменит знак на противоположный.

Мы можем, подбирая коэффициенты, складывать первую строку матрицы со всеми остальными и получать в первом столбце нули во всех позициях, кроме первой. Для получения нуля во второй строке, нужно прибавить ко второй строке первую, умноженную на

$-{a}_{2,1}/{a}_{1,1}$

Для получения нуля в третьей строке, нужно к третьей строке прибавить первую строку, умноженную на

$-{a}_{3,1}/{a}_{1,1}$

и т.д. В конечном итоге, матрица приведется к виду, в котором все элементы

${a}_{n,1}$

при n>1 будут равны нулю.

Если же в матрице элемент

${a}_{1,1}$

оказался равным нулю, то можно найти в первом столбце ненулевой элемент (предположим, он оказался на k-м месте) и обменять местами первую и k-ю строки. При этом преобразовании определитель просто поменяет знак, что можно учесть. Если же в первом столбце нет ненулевых элементов, то определитель равен нулю.

Далее, действуя аналогично, можно получить нули во втором столбце, затем в третьем и т.п. Важно, что при сложении строк полученные ранее нули не изменятся. Если для какой-либо строки не удастся найти ненулевой элемент для знаменателя, то определитель равен нулю и процесс можно остановить. Нормальное завершение процесса Гаусса порождает матрицу, у которой все элементы, расположенные ниже главной диагонали, равны нулю. Как говорилось выше, определитель такой матрицы равен произведению диагональных элементов.

Перейдем к программированию.


Мы работаем с данными с плавающей точкой. Матрицы представляем списками строк. Для начала определим два типа:

type Row    = [Double]
type Matrix = [Row]

Простая рекурсия


Ничтоже сумняшеся, мы будем раскладывать определитель по элементам первой (т.е. нулевой) строки. Нам понадобится программа построения минора, получающегося вычеркиванием первой строки и k-го столбца. 

-- Удаление k-го элемента изо всех строк матрицы
deln :: Matrix -> Int -> Matrix
deln matrix k = map (\ r -> (take (k) r)++(drop (k+1) r)) matrix

А вот и минор:

-- Минор k-го элемента нулевой строки
minor :: Matrix -> Int -> Double
minor matrix k = det $ deln (drop 1 matrix) k

Обратите внимание: минор — это определитель. Мы вызываем функцию det, которую еще не реализовали. Для реализации det, нам придется сформировать знакочередующуюся сумму произведений очередного элемента первой строки на определитель очередного минора. Ч��обы избежать громоздких выражений, создадим для формирования знака суммы отдельную функцию:

sgn :: Int -> Double
sgn n = if n `rem` 2 == 0 then 1.0 else (-1.0)

Теперь можно вычислить определитель:

-- Определитель квадратной матрицы
det :: Matrix -> Double
det [[a,b],[c,d]] = a*d-b*c
det matrix = sum $ map (\c -> ((matrix !! 0)!!c)*(sgn c)*(minor matrix c))  [0..n]
             where n = length matrix - 1

Код очень прост и не требует особых комментариев. Чтобы проверить работоспособность наших функций, напишем функцию main:

main = print $ det [[1,2,3],[4,5,6],[7,8,(-9)]]

Значение этого определителя равно 54, в чем можно убедиться.

Метод Гаусса


Нам понадобится несколько служебных функций (которые можно будет использовать и в других местах). Первая из них — взаимный обмен двух строк в матрице:

-- Обмен двух строк матрицы
swap :: Matrix -> Int -> Int -> Matrix
swap matrix n1 n2 = map row [0..n]
                    where n=length matrix - 1
                          row k | k==n1 = matrix !! n2
                                | k==n2 = matrix !! n1
                                | otherwise = matrix !! k

Как можно понять по приведенному выше коду, функция проходит строку за строкой. При этом, если встретилась строка с номером n1, принудительно подставляется строка n2 (и наоборот). Остальные строки остаются на месте.

Следующая функция вычисляет строку r1 сложенную со строкой r2, умноженной поэлементно на число f:

-- Вычислить строку r1+f*r2
comb :: Row -> Row -> Double -> Row
comb r1 r2 f = zipWith (\ x y -> x+f*y) r1 r2

Здесь все предельно прозрачно: действия выполняются над строками матрицы (т.е. над списками [Double]). А вот следующая функция в��полняет это преобразование над матрицей (и, естественно, получает новую матрицу):

-- прибавить к строке r1 строку r2, умноженную на f
trans :: Matrix -> Int -> Int -> Double -> Matrix
trans matrix n1 n2 f = map row [0..n]
                       where n=length matrix - 1
                             row k | k==n1 = comb (matrix !! n1) (matrix !! n2) f
                                   | otherwise = matrix !! k

Функция getNz ищет номер первого ненулевого элемента в списке. Она нужна в случае, когда очередной диагональный элемент оказался равным нулю.

-- Номер первого ненулевого в списке
getNz :: Row -> Int
getNz xs = if length tmp == 0 then (-1) else snd $ head tmp
           where tmp=dropWhile (\ (x,k) -> (abs x) <= 1.0e-10) $ zip xs [0..]

Если все элементы списка равны нулю, функция вернет -1.

Функция search проверяет, подходит ли матрица для очередного преобразования (у нее должен быть ненулевым очередной диагональный элемент). Если это не так, матрица преобразовывается перестановкой строк.

-- Поиск ведущего элемента и перестановка строк при необходимости
search :: Matrix -> Int -> Matrix
search matrix k | (abs ((matrix !! k) !! k)) > 1.0e-10 = matrix
                | nz < 0 = matrix  -- матрица вырождена    
                | otherwise = swap matrix k p 
                           where n   = length matrix
                                 lst = map (\ r -> r !! k) $ drop k matrix
                                 nz  = getNz lst
                                 p   = k + nz

Если ведущий (ненулевой) элемент найти невозможно (матрица вырождена), то функция вернет ее без изменений. Функция mkzero формирует нули в очередном столбце матрицы:

-- получение нулей в нужном столбце
mkzero :: Matrix -> Int -> Int -> Matrix
mkzero matrix k p | p>n-1 = matrix
                  | otherwise = mkzero (trans matrix p k (-f)) k (p+1)
                    where n = length matrix
                          f = ((matrix !! p) !! k)/((matrix !! k) !! k)

Функция triangle формирует верхнюю треугольную форму матрицы:

-- Получение верхней треугольной формы матрицы
triangle :: Matrix -> Int -> Matrix
triangle matrix k | k>=n = matrix
                  | (abs v) <= 1.0e-10 = [[0.0]] -- матрица вырождена
                  | otherwise = triangle (mkzero tmp k k1) k1 
                    where n   = length matrix
                          tmp = search matrix k
                          v   = (tmp !! k) !! k -- диагональный элемент
                          k1  = k+1

Если на очередном этапе не удалось найти ведущий элемент, функция возвращает нулевую матрицу 1-го порядка. Теперь можно составить парадную функцию приведения матрицы к верхней треугольной форме:

-- Парадная функция
gauss :: Matrix -> Matrix
gauss matrix = triangle matrix 0

Для вычисления определителя нам нужно перемножить диагональные элементы. Для этого составим отдельную функцию:

-- Произведение диагональных элементов
proddiag :: Matrix -> Double
proddiag matrix = product $ map (\ (r,k) -> r !!k) $ zip matrix [0,1..]

Ну, и «бантик» — собственно вычисление определителя:

-- Вычисление определителя
det :: Matrix -> Double
det matrix = proddiag $ triangle matrix 0

Проверим, как работает эта функция:

main = print $ det  [[1,2,3],[4,5,6],[7,8,-9]]

[1 of 1] Compiling Main             ( main.hs, main.o ) 
Linking a.out ...                                                                                 
54.0

Спасибо тем, кто дочитал до конца!

Код можно скачать здесь