E127. Вычисление определителя методом Краута за O (N3)
emaxx algorithm
#algorithm#emaxx#linear-algebra#math#matrix
Task
Источник: e-maxx.ru/algo, страница PDF 417.
Здесь будет рассмотрена модификация метода Краута (Crout), позволяющая вычислить определитель матрицы за O (N3).
Собственно алгоритм Краута находит разложение матрицы A в виде A = L U, где L - нижняя, а U - верхняя треугольная матрицы. Без ограничения общности можно считать, что в L все диагональные элементы равны 1. Но, зная эти матрицы, легко вычислить определитель A: он будет равен произведению всех элементов, стоящих на главной диагонали матрицы U. Имеется теорема, согласно которой любая обратимая матрица обладает LU-разложением, и притом единственным, тогда и только тогда, когда все её главные миноры отличны от нуля. Следует напомнить, что мы рассматриваем только такие разложения, в которых диагональ L состоит только из единиц; иначе же, вообще говоря, разложение не единственно.
Пусть A - матрица, N - её размер. Мы найдём элементы матриц L и U. Сам алгоритм состоит из следующих шагов:
1. Положим Li i = 1 для i = 1, 2, ..., N
2. Для каждого j = 1, 2, ..., N выполним: 1. Для i = 1, 2, ..., j найдём значение Ui j:
Ui j = Ai j - SUM Li k Uk j ,
где сумма по всем k = 1, 2, ..., i-1. 2. Далее, для i = j+1, j+2, ..., N имеем:
Li j = (Ai j - SUM Li k Uk j) / Uj j ,
где сумма берётся по всем k = 1, 2, ..., j-1.
Реализация
Код на Java (с использованием дробной длинной арифметики):
static BigInteger det (BigDecimal a [][], int n)
{
try {
for (int i=0; i<n; i++){
boolean nonzero = false;
for (int j=0; j<n; j++)if (a[i][j].compareTo (new BigDecimal(BigInteger.ZERO)) > 0)
nonzero = true;
if (!nonzero)return BigInteger.ZERO;}
BigDecimal scaling [] = new BigDecimal [n];
for (int i=0; i<n; i++){
BigDecimal big = new BigDecimal (BigInteger.ZERO);
for (int j=0; j<n; j++)if (a[i][j].abs().compareTo (big) > 0)big = a[i][j].abs();
scaling[i] = (new BigDecimal (BigInteger.ONE)) .divide
(big, 100, BigDecimal.ROUND_HALF_EVEN);
}
int sign = 1;for (int j=0; j<n; j++){
for (int i=0; i<j; i++){
BigDecimal sum = a[i][j];
for (int k=0; k<i; k++)sum = sum.subtract (a[i][k].multiply (a[k][j]));
a[i][j] = sum;
}
BigDecimal big = new BigDecimal (BigInteger.ZERO);
int imax = -1;for (int i=j; i<n; i++){
BigDecimal sum = a[i][j];
for (int k=0; k<j; k++)sum = sum.subtract (a[i][k].multiply (a[k][j]));
a[i][j] = sum;
BigDecimal cur = sum.abs();
cur = cur.multiply (scaling[i]);
if (cur.compareTo (big) >= 0){
big = cur;
imax = i;
} }
if (j != imax){
for (int k=0; k<n; k++){
BigDecimal t = a[j][k];
a[j][k] = a[imax][k];
a[imax][k] = t;
}
BigDecimal t = scaling[imax];
scaling[imax] = scaling[j];
scaling[j] = t;
sign = -sign;
}
if (j != n-1)for (int i=j+1; i<n; i++)a[i][j] = a[i][j].divide
(a[j][j], 100,
BigDecimal.ROUND_HALF_EVEN);
}
BigDecimal result = new BigDecimal (1);
if (sign == -1)result = result.negate();
for (int i=0; i<n; i++)result = result.multiply (a[i][i]);
return result.divide(BigDecimal.valueOf(1), 0, BigDecimal.
ROUND_HALF_EVEN).toBigInteger();
}
catch (Exception e)
{
return BigInteger.ZERO;} }
C# solution
auto-draft, review before submitusing System;
using System.Collections.Generic;
using System.Linq;
public static class AlgorithmDraft
{
// Auto-generated C# draft from the original e-maxx C/C++ listing. Review before production use.
static BigInteger det (BigDecimal a [][], int n)
{
try {
for (int i=0; i<n; i++)
{
boolean nonzero = false;
for (int j=0; j<n; j++)
if (a[i][j].compareTo (new BigDecimal
(BigInteger.ZERO)) > 0)
nonzero = true;
if (!nonzero)
return BigInteger.ZERO;
}
BigDecimal scaling [] = new BigDecimal [n];
for (int i=0; i<n; i++)
{
BigDecimal big = new BigDecimal (BigInteger.ZERO);
for (int j=0; j<n; j++)
if (a[i][j].abs().compareTo (big) > 0)
big = a[i][j].abs();
scaling[i] = (new BigDecimal (BigInteger.ONE)) .divide
(big, 100, BigDecimal.ROUND_HALF_EVEN);
}
int sign = 1;
for (int j=0; j<n; j++)
{
for (int i=0; i<j; i++)
{
BigDecimal sum = a[i][j];
for (int k=0; k<i; k++)
sum = sum.subtract (a[i][k].multiply (a[k][j]));
a[i][j] = sum;
}
BigDecimal big = new BigDecimal (BigInteger.ZERO);
int imax = -1;
for (int i=j; i<n; i++)
{
BigDecimal sum = a[i][j];
for (int k=0; k<j; k++)
sum = sum.subtract (a[i][k].multiply (a[k][j]));
a[i][j] = sum;
BigDecimal cur = sum.abs();
cur = cur.multiply (scaling[i]);
if (cur.compareTo (big) >= 0)
{
big = cur;
imax = i;
}
}
if (j != imax)
{
for (int k=0; k<n; k++)
{
BigDecimal t = a[j][k];
a[j][k] = a[imax][k];
a[imax][k] = t;
}
BigDecimal t = scaling[imax];
scaling[imax] = scaling[j];
scaling[j] = t;
sign = -sign;
}
if (j != n-1)
for (int i=j+1; i<n; i++)
a[i][j] = a[i][j].divide
(a[j][j], 100,
BigDecimal.ROUND_HALF_EVEN);
}
BigDecimal result = new BigDecimal (1);
if (sign == -1)
result = result.negate();
for (int i=0; i<n; i++)
result = result.multiply (a[i][i]);
return result.divide
(BigDecimal.valueOf(1), 0, BigDecimal.
ROUND_HALF_EVEN).toBigInteger();
}
catch (Exception e)
{
return BigInteger.ZERO;
}
}
}C++ solution
matched/originalstatic BigInteger det (BigDecimal a [][], int n)
{
try {
for (int i=0; i<n; i++)
{
boolean nonzero = false;
for (int j=0; j<n; j++)
if (a[i][j].compareTo (new BigDecimal
(BigInteger.ZERO)) > 0)
nonzero = true;
if (!nonzero)
return BigInteger.ZERO;
}
BigDecimal scaling [] = new BigDecimal [n];
for (int i=0; i<n; i++)
{
BigDecimal big = new BigDecimal (BigInteger.ZERO);
for (int j=0; j<n; j++)
if (a[i][j].abs().compareTo (big) > 0)
big = a[i][j].abs();
scaling[i] = (new BigDecimal (BigInteger.ONE)) .divide
(big, 100, BigDecimal.ROUND_HALF_EVEN);
}
int sign = 1;
for (int j=0; j<n; j++)
{
for (int i=0; i<j; i++)
{
BigDecimal sum = a[i][j];
for (int k=0; k<i; k++)
sum = sum.subtract (a[i][k].multiply (a[k][j]));
a[i][j] = sum;
}
BigDecimal big = new BigDecimal (BigInteger.ZERO);
int imax = -1;
for (int i=j; i<n; i++)
{
BigDecimal sum = a[i][j];
for (int k=0; k<j; k++)
sum = sum.subtract (a[i][k].multiply (a[k][j]));
a[i][j] = sum;
BigDecimal cur = sum.abs();
cur = cur.multiply (scaling[i]);
if (cur.compareTo (big) >= 0)
{
big = cur;
imax = i;
}
}
if (j != imax)
{
for (int k=0; k<n; k++)
{
BigDecimal t = a[j][k];
a[j][k] = a[imax][k];
a[imax][k] = t;
}
BigDecimal t = scaling[imax];
scaling[imax] = scaling[j];
scaling[j] = t;
sign = -sign;
}
if (j != n-1)
for (int i=j+1; i<n; i++)
a[i][j] = a[i][j].divide
(a[j][j], 100,
BigDecimal.ROUND_HALF_EVEN);
}
BigDecimal result = new BigDecimal (1);
if (sign == -1)
result = result.negate();
for (int i=0; i<n; i++)
result = result.multiply (a[i][i]);
return result.divide
(BigDecimal.valueOf(1), 0, BigDecimal.
ROUND_HALF_EVEN).toBigInteger();
}
catch (Exception e)
{
return BigInteger.ZERO;
}
}Java solution
auto-draft, review before submitimport java.util.*;
import java.math.*;
public class AlgorithmDraft {
// Auto-generated Java draft from the original e-maxx C/C++ listing. Review before production use.
static BigInteger det (BigDecimal a [][], int n)
{
try {
for (int i=0; i<n; i++)
{
booleanean nonzero = false;
for (int j=0; j<n; j++)
if (a[i][j].compareTo (new BigDecimal
(BigInteger.ZERO)) > 0)
nonzero = true;
if (!nonzero)
return BigInteger.ZERO;
}
BigDecimal scaling [] = new BigDecimal [n];
for (int i=0; i<n; i++)
{
BigDecimal big = new BigDecimal (BigInteger.ZERO);
for (int j=0; j<n; j++)
if (a[i][j].abs().compareTo (big) > 0)
big = a[i][j].abs();
scaling[i] = (new BigDecimal (BigInteger.ONE)) .divide
(big, 100, BigDecimal.ROUND_HALF_EVEN);
}
int sign = 1;
for (int j=0; j<n; j++)
{
for (int i=0; i<j; i++)
{
BigDecimal sum = a[i][j];
for (int k=0; k<i; k++)
sum = sum.subtract (a[i][k].multiply (a[k][j]));
a[i][j] = sum;
}
BigDecimal big = new BigDecimal (BigInteger.ZERO);
int imax = -1;
for (int i=j; i<n; i++)
{
BigDecimal sum = a[i][j];
for (int k=0; k<j; k++)
sum = sum.subtract (a[i][k].multiply (a[k][j]));
a[i][j] = sum;
BigDecimal cur = sum.abs();
cur = cur.multiply (scaling[i]);
if (cur.compareTo (big) >= 0)
{
big = cur;
imax = i;
}
}
if (j != imax)
{
for (int k=0; k<n; k++)
{
BigDecimal t = a[j][k];
a[j][k] = a[imax][k];
a[imax][k] = t;
}
BigDecimal t = scaling[imax];
scaling[imax] = scaling[j];
scaling[j] = t;
sign = -sign;
}
if (j != n-1)
for (int i=j+1; i<n; i++)
a[i][j] = a[i][j].divide
(a[j][j], 100,
BigDecimal.ROUND_HALF_EVEN);
}
BigDecimal result = new BigDecimal (1);
if (sign == -1)
result = result.negate();
for (int i=0; i<n; i++)
result = result.multiply (a[i][i]);
return result.divide
(BigDecimal.valueOf(1), 0, BigDecimal.
ROUND_HALF_EVEN).toBigInteger();
}
catch (Exception e)
{
return BigInteger.ZERO;
}
}
}Explanation
Материал разбит как алгоритмическая задача: изучить постановку, понять асимптотику и реализовать алгоритм на выбранном языке.