Unit Algemeen;
{
  PunchCard Reader
  ================
  This program has been designed and is CopyLefted by:

* Han de Bruijn; Systems for Research "A little bit of Physics would be  (===)
* and Education (DTO/SOO), Mekelweg 6  NO Idleness in Mathematics"(HdB) @-O^O-@
* 2628 CD Delft, The Netherlands  http://huizen.dto.tudelft.nl/deBruijn  #/_\#
* E-mail: Han.deBruijn@DTO.TUDelft.NL  Tel: +31 15 27 82751. Fax: 81722   ###

All I ask is to be credited when it is appropriate.

  Mathematics in General
  ----------------------
}
INTERFACE

Uses SysUtils;

type
  integers = array of integer;
  vektor = record
    x,y : double;
  end;
  Jordan = array of vektor;
  punt = record
    x,y : integer;
  end;
  punten = array of punt;
  tensor = record
    xx,yy,xy : double;
  end;
  kolom = array of double;
  matrix = array of array of double;

procedure LogFile(naam : string); { Name of file }
procedure Schijf(bericht : string); { for Logging }
procedure TijdMeting(afdruk : string); { & Timing }

function Hoek(p,q : vektor) : double; { Angle between Vectors }

function TelSort(rij : array of double) : integers; { Sort by Counting }

function Verplaats(p,dd : vektor; pm : integer) : vektor; { Translate }
function Roteer(p,dd : vektor; pm : integer) : vektor; { Rotate }
function Schaal(p,dd : vektor) : vektor; { Scale }

function Faculteit(n : integer) : integer; { Factorial }
{ Permutations in lexographic order. Algorithm by Edsger W. Dijkstra }
function Dijkstra(waarde : integers; var teken : integer) : integers;
procedure Edsger(N : integer); { Just a Test }

type
  VektorRuimte = class(TObject) { Vector Space }
  private
    NDM : integer; { Dimension }
  public
    procedure Dimensie(grootte : integer); { Define Dimension }
    function DirekteInverse(var getallen : matrix) : boolean; { Direct Inverse }
    function Vermenigvuldig(een : matrix; twee : kolom) : kolom; { Multiply }
  end;

  Schuine = class(TObject) { Skewed 2-D Gaussian Distribution }
  private
    sigma : tensor; { 2nd order Moments }
    mu : vektor; { 1st order Moments }
  public
    procedure Definitie(S : tensor; M : vektor);
    function Gewogen : tensor;
    function Gauss(r : vektor) : double;
  end;

procedure DirectDemo(NDM : integer); { Just a Test }

IMPLEMENTATION

var
  genoemd : string;

procedure LogFile;
{
  Name of LogFile
}
begin
  genoemd := naam;
end;

procedure Schijf;
{
  (Debugging and) Logging
}
var
  bug : text;
  eerst : boolean;
begin
  eerst := not FileExists(genoemd);
  AssignFile(bug,genoemd);
  if eerst then Rewrite(bug)
           else Append(bug);
  Writeln(bug, bericht);
  Closefile(bug);
  Writeln(bericht);
end;

procedure TijdMeting;
{
  StopWatch
}
const
  eerder : integer = 0;
var
  heden : integer;
begin
  if afdruk = '' then
  begin
    eerder := DateTimeToTimeStamp(Now).Time;
    Exit;
  end;
  heden := DateTimeToTimeStamp(Now).Time;
  Schijf(afdruk+' = '+IntToStr(heden-eerder)+' ms');
end;

function Hoek;
{
  Calculates the Angle between two Vectors
}
var
  uit,inp,len1,len2,len12 : double;
  a,b,s : vektor;
begin
  Hoek := 0;

{ Outer product (area parallellogram): }
  uit := p.x*q.y - q.x*p.y;
{ Area parallellogram also equals two times
  the distance to the boundary, times half
  the base of the triangle (= dist*len12) }
  len12 := sqrt(sqr(q.x-p.x) + sqr(q.y-p.y));
{ Distance to boundary < 1 pixel too close: }
  if abs(uit) < len12 * 1 then Exit;

{ Use sum of normed Vectors: }
  len1 := sqrt(sqr(p.x) + sqr(p.y));
  len2 := sqrt(sqr(q.x) + sqr(q.y));
  if (len1 = 0) or (len2 = 0) then Exit;

  a.x := p.x/len1 ; a.y := p.y/len1;
  b.x := q.x/len2 ; b.y := q.y/len2;
  s.x := a.x + b.x;
  s.y := a.y + b.y;
{ The triangle spanned by (a) and (b) is equilateral.
  Hence the sum-vector (a+b) divides the top angle in
  two equal pieces. Take one of them and multiply the
  angle found by two. But first determine outer/inner
  products: }
  uit := a.x*s.y - s.x*a.y;
  inp := a.x*s.x + a.y*s.y;
{ The quotient of these is sin(Hoek)/cos(Hoek). Hence }
  Hoek := 2 * arctan(uit/inp);
end;

function TelSort;
{
  Sorting by counting.
  Algorithm C on page 76 of "The Art of Computer Programming"
  volume 3 / "Sorting and Searching" by Donald E. Knuth
}
var
  veel, i,j,k,L : integer;
  tel,let : integers;
begin
  veel := Length(rij)-1;
  SetLength(tel,veel+1);
  for k := 1 to veel do
    tel[k] := 0;

  for k := 0 to (veel-2) do
  begin
    i := veel - k;
    for L := 1 to (i-1) do
    begin
      j := i - L;
      if rij[i] < rij[j] then
         tel[j] := tel[j] + 1
      else
         tel[i] := tel[i] + 1;
    end;
  end;

  SetLength(let,veel+1);
  for k := 1 to veel do
    let[tel[k]+1] := k;
  SetLength(tel,0);

  TelSort := let;
  SetLength(let,0);
end;

function Verplaats;
{
  Translation
}
var
  r : vektor;
begin
  r.x := p.x + pm*dd.x;
  r.y := p.y + pm*dd.y;
  Verplaats := r;
end;

function Roteer;
{
  Rotation
}
var
  c,s : double;
  q,r : vektor;
begin
  c := dd.x;
  s := dd.y;
  q := p;
  r.x := c*q.x - pm*s*q.y;
  r.y := pm*s*q.x + c*q.y;
  Roteer := r;
end;

function Schaal;
{
  Scaling
}
var
  r : vektor;
begin
  r.x := p.x*dd.x;
  r.y := p.y*dd.y;
  Schaal := r;
end;

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

function Faculteit;
{
  Factorial n!
  ------------ }
var
  k,m : integer;
begin
  m := 1;
  for k := 2 to n do
    m := m*k;
  Faculteit := m;
end;

function Dijkstra;
{
  Should produce a Permutation (waarde)
  in alphabetical (lexicographic) order.
  Start with the identical permutation.

  According to: E.W.Dijkstra,
  "A discipline of Programming",
  Prentice Hall, 1976
}
var
  i,j,N : integer;

  procedure Wissel(k,L : integer);
  var
    h : integer;
  begin
    h := waarde[L];
    waarde[L] := waarde[k];
    waarde[k] := h;
    teken := - teken;
  end;

begin
  N := Length(waarde);

  i := N-1;
  while (waarde[i-1] >= waarde[i]) do i := i-1;
  j := N;
  while (waarde[j-1] <= waarde[i-1]) do j := j-1;

  Wissel(i-1,j-1);
  i := i+1; j := N;

  while (i < j) do
  begin
    Wissel(i-1,j-1);
    i := i+1;
    j := j-1;
  end;

  Dijkstra := waarde;
end;

procedure Edsger;
{
  Test Permutations
}
var
  orde,F,k,i,teken : integer;
  rij : integers;

  function Extra : string;
  var
    s : string;
  begin
    s := '+ ';
    if teken < 0 then s := '- ';
    Extra := s;
  end;

begin
  orde := N;
  SetLength(rij,orde);
  F := Faculteit(orde);
  for i := 0 to orde-1 do
    rij[i] := i+1;
  teken := +1;
  for k := 1 to F do
  begin
    if k > 1 then
      rij := Dijkstra(rij,teken);
    Write(Extra);
    for i := 0 to orde-1 do
      Write(rij[i]:3);
    Writeln;
  end;
end;

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

  function teken(L : integer) : integer;
  begin
    teken := 1-2*(L mod 2);
  end;

begin
  subdet := 0;
  fout := (ndm < 1) or (ndm > 8);
  if fout then 
  begin
    Writeln('subdet: input Error');
    Exit;
  end;

  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 VektorRuimte.Dimensie;
{
  Dimension of Vector Space
}
begin
  NDM := grootte;
end;

function VektorRuimte.DirekteInverse;
{
  Direct inversion of small matrices
  ---------------------------------- }
var
  i,j : integer;
  invert : matrix;
  det : double;
begin
  DirekteInverse := false;
  if ndm < 1 then 
  begin
    Writeln('DirekteInverse: wrong input');
    Exit;
  end;
  SetLength(invert,ndm+1,ndm+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];
  if det = 0 then 
  begin
    Writeln('DirekteInverse: singular matrix');
    Exit;
  end;

  for i := 1 to ndm do
  for j := 1 to ndm do
    invert[i,j] := invert[i,j]/det;

  getallen := invert;
  DirekteInverse := true;
end;

procedure Scherm(ndm : integer; beeld : matrix);
{
  Print
}
var
  k,L : integer;
begin
  for L := 1 to ndm do
  begin
    for k := 1 to ndm do
      Write(' ',beeld[k][L]:8:5);
    Writeln;
  end;
  Writeln;
end;

procedure DirectDemo;
{         ==========
  Direct inversion
  ----------------
  demonstration program
  --------------------- }
var
  a,b,c : matrix;
  i,j,k : integer;
  h : double;
  ruimte : VektorRuimte;
begin
  SetLength(a,NDM+1,NDM+1);
  SetLength(b,NDM+1,NDM+1);
  SetLength(c,NDM+1,NDM+1);

{ Random matrix }
  for i := 1 to NDM do
  for j := 1 to NDM do
    a[i,j] := Random;
  Scherm(NDM,a);
  b := a;

  ruimte := VektorRuimte.Create;
  ruimte.Dimensie(NDM);
{ Direct inversion }
  ruimte.DirekteInverse(b);
  Scherm(NDM,b);

{ Multiply original with inverse }
  for i := 1 to NDM do
  begin
    for j := 1 to NDM do
    begin
      h := 0.;
      for k := 1 to NDM do
        h := h + a[i,k]*b[k,j];
      c[i,j] := h;
    end;
  end;

{ Result must be unity }
  Scherm(NDM,c);
end;

function VektorRuimte.Vermenigvuldig;
{
  Matrix Vector Multiply
}
var
  uit : kolom;
  i,j,N : integer;
  h : double;
begin
  N := ndm;
  SetLength(uit,N+1);
  for i := 1 to N do
  begin
    h := 0;
    for j := 1 to N do
      h := h + een[i,j]*twee[j];
    uit[i] := h;
  end;
  Vermenigvuldig := uit;
end;

procedure Schuine.Definitie;
{
  Definition
}
var
  Det : double;
begin
  sigma := S;
  mu := M;
  Det := sigma.xx*sigma.yy 
       - sqr(sigma.xy);
  sigma.xx := sigma.xx/Det;
  sigma.yy := sigma.yy/Det;
  sigma.xy := sigma.xy/Det;
end;

function Schuine.Gewogen;
begin
  Gewogen := sigma;
end;

function Schuine.Gauss;
{
  Skewed 2-D Gaussian
}
var
  vgl : double;
begin
  Gauss := 0;
  vgl := sigma.yy*sqr(r.x - mu.x)
       - 2*sigma.xy*(r.x - mu.x)*(r.y - mu.y)
       + sigma.xx*sqr(r.y - mu.y);
  if vgl < 0 then Writeln('ALARM');
{ If exponential function near zero }
  if vgl > sqr(2*pi) then Exit;
  Gauss := exp(-vgl/2);
end;

END.