Unit Numeriek;
{
Automated Balancing of Chemical Equations
=========================================
This software has been designed and it is
CopyLefted by Han de Bruijn:

                          (===)
                         @-O^O-@
                          #/_\#
                           ###

Solver [0] Fractions Arithmetic (default)
Solver [1] Floating Point Gauss Solver
Solver [2] Cramer's Rule (# molecules < 9)
}
INTERFACE

Uses Algemeen, Breuken;

const
  debug : boolean = FALSE;

procedure Bsolve(A : Bmatrix; var b : Bvektor); { Fractions Arithmetic (default) }
procedure Dsolve(A : Dmatrix; var b : Dvektor; var Det : double); { Floating Point }
procedure DirekteInverse(var onze : matrix; var Det : integer); { Cramer's Rule }
procedure Vermenigvuldig(een : matrix; var twee : vektor); { Matrix times vector }

IMPLEMENTATION

function subdet(ii,jj,ndm : integer; getallen : matrix) : integer;
{        ======
  Subdeterminant of small matrix
}
var
  rij : vektor;
  det, term : integer;
  m,i,j,k, even : integer;

  function teken(L : integer) : integer;
  {
    Plus (+1) or minus sign (-1)
  }
  begin
    teken := 1-2*(L mod 2);
  end;

  function als(prop : boolean) : integer;
  {        ===
    IF like in APL
  }
  begin
    als := 0;
    if prop then als := 1;
  end;

begin
  SetLength(getallen,ndm+1,ndm+1);
  SetLength(rij,ndm-1);
  for k := 0 to ndm-2 do
    rij[k] := k+1;

  det := 0;
  even := +1;
  for m := 1 to Faculteit(ndm-1) do
  begin
    if m > 1 then
      rij := Dijkstra(rij,even);
    term := even;
    for k := 1 to ndm-1 do
    begin
      i := k + als(k >= ii);
      j := rij[k-1] + als(rij[k-1] >= jj);
      term := term*getallen[i,j];
    end;
    det := det + term;
  end;
  subdet := teken(ii+jj) * det;
end;

procedure DirekteInverse(var onze : matrix; var Det : integer);
{
  Direct inversion of small matrices
  ---------------------------------- }
var
  i,j,ndm : integer;
  invert : matrix;
  getallen : matrix;
begin
  ndm := Length(onze);
  if ndm < 1 then 
  begin
    Writeln('DirekteInverse: wrong input');
    Exit;
  end;
  if (ndm > 9) then
  begin
    Writeln(' Error: factorial has become too large!');
    Exit;
  end;

  SetLength(invert,ndm+1,ndm+1);
  SetLength(getallen,ndm+1,ndm+1);

  for i := 1 to ndm do
  for j := 1 to ndm do
     getallen[i,j] := onze[i-1,j-1];

  for i := 1 to ndm do
  for j := 1 to ndm do
    invert[i,j] := subdet(j,i,ndm,getallen);

  det := 0;
  for i := 1 to ndm do
    det := det + getallen[i,1]*invert[1,i];

  for i := 1 to ndm do
  for j := 1 to ndm do
    onze[i-1,j-1] := invert[i,j];
end;

procedure Vermenigvuldig(een : matrix; var twee : vektor);
{
  Matrix Vector Multiply
}
var
  uit : vektor;
  i,j,N : integer;
  h : integer;
begin
  N := Length(een);
  SetLength(uit,N);
  for i := 0 to N-1 do
  begin
    h := 0;
    for j := 0 to N-1 do
      h := h + een[i,j]*twee[j];
    uit[i] := h;
  end;
  twee := uit;
end;

procedure Dsolve(A : Dmatrix; var b : Dvektor; var Det : double);
{
  Simple Gauss Solver with
  all pivots on diagonal
  Floating point, adapted
}
var
  N,i,j,k : integer;
  p,s : double;
begin
  N := Length(b);
{ Elimination }
  Det := 1; p := 0;
  for k := 0 to N-1 do
  begin
    p := A[k,k];
    if p = 0 then Break;
    for i := k+1 to N-1 do
    begin
      if A[i,k] = 0 then Continue;
      A[i,k] := A[i,k]/p;
      for j := k+1 to N-1 do
        A[i,j] := A[i,j] - A[i,k]*A[k,j];
      b[i] := b[i] - A[i,k]*b[k];
    end;
    Det := Det*p;
  end;
{ Backsubstitution }
  if p = 0 then Det := 0;
  if p > 0 then
  for k := N-1 downto 0 do
  begin
    s := b[k];
    for j := k+1 to N-1 do
      s := s - A[k,j]*b[j];
    b[k] := s/A[k,k];
  end;
end;

procedure Bsolve(A : Bmatrix; var b : Bvektor);
{
  Simple Gauss Solver
  -------------------
  Pivots on diagonal
  with true Fractions
}
var
  N,i,j,k : integer;
  p,s,Det : breuk;
begin
  N := Length(b);
{ Elimination }
  for k := 0 to N-1 do
  begin
    p := A[k,k];
    for i := k+1 to N-1 do
    begin
      if A[i,k].teller = 0 then Continue;
      A[i,k] := deel(A[i,k],p);
      for j := k+1 to N-1 do
        A[i,j] := minus(A[i,j],maal(A[i,k],A[k,j]));
      b[i] := minus(b[i],maal(A[i,k],b[k]));
    end;
  end;
  if debug then
  begin
    Det := maak(1,1);
    for i := 0 to N-1 do
    begin
      Det := maal(Det,A[i,i]); 
      for j := 0 to i-1 do
        Write(' ':12);
      for j := i to N-1 do
        Write(afdruk(A[i,j]):12);
      Writeln(' |',afdruk(b[i]):12);
    end;
    Writeln('  Det = ',afdruk(Det));
    Writeln;
  end;
{ Backsubstitution }
  for k := N-1 downto 0 do
  begin
    s := b[k];
    for j := k+1 to N-1 do
      s := minus(s,maal(A[k,j],b[j]));
    b[k] := deel(s,A[k,k]);
  end;
end;

END.