Wyznacznik macierzy

W jaki sposób mogę obliczyć wyznacznik macierzy kwadratowej (n x n). Powiedzmy, ze jest ona zawarta w tablicy dynamicznej (n x n). Mam na mysli oczywiscie sytuacje ktora zarowno bedzie dzialac na macierzy 3x3, (ktora mozna rozwalic schematem Sarrusa swoja droga) jak i dowolnie wiekszej 7x7, 11x11...

type { you can modify this dimension }
arr = array[1..4, 1..4] of real;

procedure printarr(new_a: arr; dim: integer);
var i, j: integer;
begin
for i := 1 to dim do
begin
for j := 1 to dim do
write(new_a[i,j]3,' ');
writeln;
end;
end;

{ Gauss-Jordan Elimination }
function gauss(a: arr; dim: integer): real;
var
new_a : arr;
i, j, k : integer;
factor, temp, det : real;
begin
{ copy the array }
for i := 1 to dim do
for j := 1 to dim do
new_a[i,j] := a[i,j];

det := 1.0;
{ do the elimination }
for i := 1 to dim-1 do
begin
{ if the main diagonal value is zero }
{ re-sort the array }
if (new_a[i,i] = 0) then
begin
for j := i+1 to dim do
begin
if (new_a[j,i] 0) then
begin
for k := 1 to dim do
begin
temp := new_a[i,k];
new_a[i,k] := new_a[j,k];
new_a[j,k] := temp;
end;
{ For Gauss-Jordan Elimination, }
{ if we do a switch, the determinant }
{ switches sign. }
det := -det;
break;
end;
end;
end;

 { if after the resorting, the value is still zero }
 { then the determinant is definitely zero }
 if (new_a[i,i] = 0) then
 begin
   gauss := 0; exit;
 end;

 { eliminate the lower rows to achieve triangular zeroes }
 for j := i+1 to dim do
 begin
   if (new_a[j,i]  0) then
   begin
     factor := (new_a[j,i] * 1.0) / new_a[i,i];
     for k := i to dim do
     begin
       new_a[j,k] := new_a[j,k] - factor * new_a[i,k];
       {2}
     end;
     {1}
   end;
 end;

end;

{ calculate the main diagonal }
for i := 1 to dim do
det := det * new_a[i,i];

gauss := det;
end;

[źródło: http://http://www.geocities.com/SiliconValley/Park/3230/misc/misc20011214-0000.html]

Kod ten znalazłem w ciągu kilku sekund używając http://google.pl :-[

{ copy the array }
for i := 1 to dim do
for j := 1 to dim do
new_a[i,j] := a[i,j];

new_a:=a;

Szybciej :)

Mozecie looknac na moj kod?
Probowalem wykorzystac powyzszy algorytm, ale sypie mi bledami z overflow
i wartosc wyznacznika nie zgadza sie z obliczonym wyznacznikiem w excelu.

program Project1;

{$APPTYPE CONSOLE}

uses
  SysUtils;

//begin
  { TODO -oUser -cConsole Main : Insert code here }


type { wymiar macierzy }
arr = array[1..4, 1..4] of real;

var
macierz : arr;
x:string;

procedure setarr(new_a: arr; dim: integer);
var i, j: integer;
begin
   for i := 1 to dim do
   begin
      for j := 1 to dim do
         macierz[i,j]:=j*random(2);
   end;
end;

procedure printarr(new_a: arr; dim: integer);
var i, j: integer;
begin
   for i := 1 to dim do
   begin
      for j := 1 to dim do
         write(new_a[i,j]:5:3,' ');
      writeln;
   end;
end;

{ eliminacja Gaussa-Jordana }
function gauss(a: arr; dim: integer): real;
var
  new_a : arr;
  i, j, k : integer;
  factor, temp, det : real;

begin
   { kopiuj tablice }
   for i := 1 to dim do
      for j := 1 to dim do
         new_a[i,j] := a[i,j];

   det := 1.0;
   { wykonaj eliminacje }
   for i := 1 to dim-1 do
   begin

     if (new_a[i,i] = 0) then
     begin
       for j := i+1 to dim do
       begin
         if (new_a[j,i] = 0) then
         begin
           for k := 1 to dim do
           begin
             temp := new_a[i,k];
             new_a[i,k] := new_a[j,k];
             new_a[j,k] := temp;
           end;
           { Eliminacja Gaussa-Jordana, }
           { jesli zamieni sie kolumny/wiersze }
           { to wyznacznik zmienia znak. }
           det := -det;
           break;
         end;
       end;
     end;

     { jesli po przestawieniu wartosc jest dalej zero }
     { wtedy wyznacznik tez jest zero }
     if (new_a[i,i] = 0) then
     begin
       gauss := 0;
       x:=floattostr(det);
     end;

     { eliminacja wierszy ponizej by uzyskac zera }
     for j := i+1 to dim do
     begin
       if (new_a[j,i] = 0) then
       begin
         factor := (new_a[j,i] * 1.0) / new_a[i,i];
         for k := i to dim do
         begin
           new_a[j,k] := new_a[j,k] - factor * new_a[i,k];
           {2}
         end;
         {1}
       end;
     end;
   end;

   { obliczenie glownego wyznacznika }
   for i := 1 to dim do
      det := det * new_a[i,i];

   gauss := det;

   x:=floattostr(det);

end;

begin
setarr(macierz,4);
printarr(macierz,4);
gauss(macierz,4);
writeln('wyznacznik = ',x);
readln;

end.