Unit Rekenen;

INTERFACE

type
  integers = array of integer;

function decimaal(rij : integers) : string; { sorted array -> decimal string }
function bitmap(getal : string) : integers; { decimal string -> sorted array }
function plus(a,b : integers) : integers; { a + b }
function maal(a,b : integers) : integers; { a * b }
function macht(x,n : integer) : string; { x^n }
function kleiner(a,b : integers) : boolean; { a < b }
function tien(n : integer) : string; { 10^n as string }
procedure wortel_twee(noemer : string); { sqrt(2) estimates }
function faculteit(N : integer) : string;

IMPLEMENTATION

function TwoMerge(rij1,rij2 : integers) : integers;
{
  Two-way merge.
  Described on page 160 of "The Art of Computer Programming"
  volume 3 / "Sorting and Searching" by Donald E. Knuth.
}
var
  LIM1,LIM2,LIM3 : integer;
  max1,max2,sentinel : integer;
  k, k1, k2, L : integer;
  rij3 : integers;
begin
  LIM1 := Length(rij1);
  LIM2 := Length(rij2);
  LIM3 := LIM1 + LIM2;
  SetLength(rij3,LIM3);
  if (LIM2 = 0) and (LIM1 > 0 )then
  for k := 0 to LIM1-1 do
    rij3[k] := rij1[k];
  if (LIM1 = 0) and (LIM2 > 0 )then
  for k := 0 to LIM2-1 do
    rij3[k] := rij2[k];
  if (LIM1 > 0) and (LIM2 > 0) then
  begin
    max1 := rij1[LIM1-1];
    max2 := rij2[LIM2-1];
    sentinel := max1;
    if max2 > max1 then sentinel := max2;
    sentinel := sentinel + 1;
  
    k1 := 0;
    k2 := 0;
    for L := 0 to LIM3-1 do
    begin
      if rij1[k1] <= rij2[k2] then
      begin
        rij3[L] := rij1[k1];
        rij1[k1] := sentinel;
        if k1 < LIM1-1 then k1 := k1+1;
      end else begin
        rij3[L] := rij2[k2];
        rij2[k2] := sentinel;
        if k2 < LIM2-1 then k2 := k2+1;
      end;
    end;
  end;
  TwoMerge := rij3;
end;

function plus(a,b : integers) : integers;
{
  Summation with sorted arrays
}
var
  q,r,rij,som : integers;
  L,k,P : integer;
begin
  L := Length(a);
  SetLength(q,L);
  for k := 0 to L-1 do
    q[k] := a[k];
  L := Length(b);
  SetLength(r,L);
  for k := 0 to L-1 do
    r[k] := b[k];
  rij := TwoMerge(q,r);
  L := Length(rij);
  if L = 0 then
  begin
    plus := rij;
    Exit;
  end;

  SetLength(som,L);
  P := 1;
  som[0] := rij[0];
  for k := 1 to L-1 do
  begin
    if rij[k] > som[P-1] then
    begin
      som[P] := rij[k]; P := P+1; Continue;
    end;

    if rij[k] = som[P-1] then
    begin
      som[P-1] := rij[k]+1; Continue;
    end;

    if rij[k] < som[P-1] then
    begin
      som[P] := som[P-1]; som[P-1] := rij[k]; P := P+1; Continue;
    end;
  end;
  SetLength(som,P);
  plus := som;
end;

function maal(a,b : integers) : integers;
{
  Multiplication with sorted arrays
}
var
  M,N,i,j,links : integer;
  h,produkt : integers;
begin
  M := Length(a);
  N := Length(b);
  SetLength(h,N);
  SetLength(produkt,0);
  for i := 0 to M-1 do
  begin
    links := a[i];
    for j := 0 to N-1 do
      h[j] := b[j] + links;
    produkt := plus(produkt,h);
  end;
  maal := produkt;
end;

function dubbel(getal : string; carry : boolean) : string;
{
  Make decimal number representation twice as large
}
var
  G,H : string;
  rest : boolean;
  c : byte;
  k,L : integer;
begin
  G := getal;
  L := Length(G);
  rest := carry;
  for k := L downto 1 do
  begin
    c := byte(G[k])-48;
    c := c*2;
    if rest then c := c+1;
    rest := (c > 9);
    if rest then c := c-10;
    G[k] := char(c+48);
   end;
   H := G;
   if rest then H := '1' + G;
   dubbel := H;
end;

function decimaal(rij : integers) : string;
{
  Convert sorted array to decimal string
}
var
  G : string;
  L,p,k : integer;
  rest : boolean;
begin
  decimaal := '0';
  L := Length(rij);
  if L = 0 then Exit;
  G := '0';
  p := L-1;
  k := rij[p];
  while (p >= 0) do
  begin
    rest := false;
    if rij[p] = k then
    begin
       rest := true;
       p := p-1;
    end;
    G := dubbel(G,rest);
    k := k-1;
  end;
  while k >= 0 do
  begin
    G := dubbel(G,false);
    k := k-1;
  end;
  decimaal := G;
end;

function let(n : integer) : integers;
{
  Convert integer to sorted array
}
var
  k,m,L : integer;
  uit : integers;
begin
  m := n;
  SetLength(uit,32);
  for k := 0 to 31 do
    uit[k] := 0;
  L := 0;
  for k := 0 to 31 do
  begin
    if (m and 1) = 1 then
    begin
      uit[L] := k;
      L := L + 1;
    end; 
    m := m shr 1;
  end; 
  SetLength(uit,L);
  let := uit;
end;

function macht(x,n : integer) : string;
{
  Power x^n
}
var
  m : integer;
  p, y : integers;
begin
  m := n;
  SetLength(p,1);
  p[0] := 0;
  SetLength(y,0);
  y := let(x);
  while m > 0 do begin
    if (m and 1) > 0 then p := maal(p,y);
    m := m shr 1;
    y := maal(y,y);
  end;
  macht := decimaal(p);
end;

function kleiner(a,b : integers) : boolean;
{
  ( a < b ) ?
}
var
  M,N,i,j : integer;
  minder : boolean;
  P,Q,R : boolean;
begin
  M := Length(a);
  N := Length(b);
{ Special cases }
  minder := false;
  P := (M = 0) and (N > 0);
  if P then minder := true;
  Q := (M > 0) and (N = 0);
  R := (M = 0) and (N = 0);
  kleiner := minder;
  if P or Q or R then Exit;
{ General case }
  i := M; j := N;
  minder := false;
  while true do
  begin
    i := i-1; j := j-1;
    if (i < 0) and (j < 0) then
    begin
      minder := false;
      Break;
    end;
    if (i < 0) and (j >= 0) then
    begin
      minder := true;
      Break;
    end;
    if (j < 0) and (i >= 0) then
    begin
      minder := false;
      Break;
    end;

    if a[i] = b[j] then Continue;
    if a[i] > b[j] then
    begin
      minder := false;
      Break;
    end;
    if a[i] < b[j] then
    begin
      minder := true;
      Break;
    end;
  end;
  kleiner := minder;
end;

function bitmap(getal : string) : integers;
{
  Convert string to sorted array
}
var
  rest : boolean;
  L,k,m,start,F,t : integer;
  G : string;
  c : byte;
  rij : integers;
begin
  F := Round(ln(10)/ln(2))+1;
  L := Round(Length(getal))*F;
  SetLength(rij,L);
{ Writeln(getal); }
  L := Length(getal);
  G := getal;
  start := 1;
  m := 0; t := 0;
  while (start < L) or (G[L] > '0') do
  begin
    rest := false;
    for k := start to L do
    begin
      c := byte(G[k])-48;
      if rest then c := c + 10;
      rest := ((c mod 2) = 1);
      c := c div 2;
      G[k] := char(c+48);
    end;
    for k := start to L do
    begin
      c := byte(G[k])-48;
      if c = 0 then Continue;
      start := k;
      Break;
    end;
  { if rest then Writeln(' ',m); }
    if rest then
    begin
      rij[t] := m;
      t := t + 1;
    end;
  { Writeln(Copy(G,start,L-start+1)); }
    m := m + 1;
  end;
  SetLength(rij,t);
  bitmap := rij;
end;

function tien(n : integer) : string;
{
  10^n as string
}
var
  t : string;
  k : integer;
begin
  t := '1';
  for k := 1 to n do
  begin
    t := t + '0';
  end;
  tien := t;
end;

procedure wortel_twee(noemer : string);
{
  Walk through the Stern-Brocot tree
  all along until sqrt(2). Modified.
}
var
  m1,m2,n1,n2,vgl,L,R : integers;
  tel : integer;
begin
  Writeln;
{ Square of sqrt(2) }
  SetLength(vgl,1);
  vgl[0] := 1; { = 2 }
{ Initialize tree }
  SetLength(m1,0); { = 0 }
  SetLength(n2,0); { = 0 }
  SetLength(n1,1);
  n1[0] := 0; { = 1 }
  SetLength(m2,1);
  m2[0] := 0; { = 1 }
  SetLength(L,0);
  SetLength(R,0);
{ # iterates }
  tel := 0;
  while true do
  begin
    tel := tel + 1;
  { if sqrt(2) < approximating fraction }
    L := plus(n1,n2);
    L := maal(L,L);
    R := plus(m1,m2); 
    R := maal(R,R);
    if kleiner(maal(L,vgl),R) then
  { Tightening }
    begin
    { Upper Bound }
      m2 := plus(m1,m2);
      n2 := plus(n1,n2);
    end else begin
    { Lower Bound }
      m1 := plus(m1,m2);
      n1 := plus(n1,n2);
    end;
  { if Upper Bound - Lower Bound < absolute tolerance }
    if kleiner(bitmap(noemer),maal(n1,n2)) then Break;
  end;
  Writeln('   ',decimaal(m1)); Writeln(' / ',decimaal(n1));
  Writeln(' < sqrt(2) < ');
  Writeln('   ',decimaal(m2)); Writeln(' / ',decimaal(n2));
  Writeln(' ',tel,' iterations');
end;

function faculteit(N : integer) : string;
{
  Factorial
}
var
  een,bij,fac : integers;
  k : integer;
begin
  SetLength(een,1);
  een[0] := 0;
  SetLength(bij,0);
  SetLength(fac,1);
  fac[0] := 0;
  for k := 1 to N do
  begin
    bij := plus(bij,een);
    fac := maal(bij,fac);
  end; 
  faculteit := decimaal(fac);
end;

END.