dynhuff.p - optimal dynamic huffman

[Graham Toal]

Archive-name: dynhuff.p

All this talk of compression prodded me into typing in a program from
an old faded photocopy that's been kicking around my pending file
for several months.  Surprisingly it worked first time, which is just
as well because if there were any bugs in the logic I'd never have
found them ;-)

---- cut here ----
{
  ALGORITHM 673,

  ACM Transactions on Mathematical Software, Vol 15, No 2, Pages 158-167. 

  Jeffrey Scott Vitter
  Brown University

  This file coded up from the paper above by Graham Toal <gtoal at ed.ac.uk>

  This is a one-pass dynamic Huffman code generator.  I supply a trivial
  interface for testing.  Real-world use would require that you firstly
  translate this into C (please post here when you do) and secondly
  modify the IO to write binary files. (This writes a text file of
  0's and 1's for demonstration purposes)

  Also it needs some logic for what happens on end of file.

}

program huff(input, output);
const
  n = 256;
  b = n*2 + 1;
var
  alpha: array [1..n] of integer;
  rep: array [1..n] of integer;
  block: array [1..b] of integer;
  weight: array [1..b] of integer;
  parent: array [1..b] of integer;
  parity: array [1..b] of integer;
  rtChild: array [1..b] of integer;
  first: array [1..b] of integer;
  last: array [1..b] of integer;
  prevBlock: array [1..b] of integer;
  nextBlock: array [1..b] of integer;
  availBlock: integer;
  stack: array [1..n] of integer;
  a: integer; c: char;
  M, E, R, Z: integer;
  outpos: integer;

procedure Initialize;
var
  i: integer;
begin
  M := 0; E := 0; R := -1; Z := 2*n - 1;
  for i := 1 to n do begin
    M := M+1; R := R+1;
    if R*2 = M then begin E := E+1; R := 0 end;
    alpha[i] := i; rep[i] := i;
  end;
  { Initialize node n as the 0-node }
  block[n] := 1; prevBlock[1] := 1; nextBlock[1] := 1; weight[1] := 0;
  first[1] := n; last[1] := n; parity[1] := 0;
  { Initialize available block list }
  availBlock := 2;
  for i := availBlock to Z-1 do nextBlock[i] := i+1;
  nextBlock[Z] := 0;
end;

procedure Transmit(i: integer);
begin
  outpos := outpos + 1;
  write(i:1);
  if outpos = 64 then begin outpos := 0; writeln end;
end;

procedure Receive: integer;
var
  a: char;
begin
  if eoln then readln;
  if eof then begin
    {writeln('Unexpected end of file');} halt;
  end;
  read(a);
  Receive := ORD(a)-ORD('0');
end;

procedure EncodeAndTransmit(j: integer);
var
  i, ii, q, t, root: integer;
begin
  q := rep[j]; i := 0;
  if q <= M then begin { Encode letter of zero weight }
    q := q - 1;
    if q < R*2 then t := E+1 else begin q := q - R; t := E end;
    for ii := 1 to t do begin
      i := i + 1; stack[i] := q mod 2;
      q := q div 2
    end;
    q := M;
  end;
  if M = n then root := n else root := Z;
  while q <> root do begin { Traverse up the tree }
    i := i + 1; stack[i] := (first[block[q]] - q + parity[block[q]]) mod 2;
    q := parent[block[q]] - (first[block[q]] - q + 1 - parity[block[q]]) div 2
  end;
  for ii := i downto 1 do Transmit(stack[ii])
end;

function FindChild(j, parity: integer): integer;
var
  delta, right, gap: integer;
begin
  delta := 2*(first[block[j]]-j) + 1 - parity;
  right := rtChild[block[j]]; gap := right - last[block[right]];
  if delta <= gap then FindChild := right-delta
  else begin
    delta := delta - gap - 1;
    right := first[prevBlock[block[right]]]; gap := right - last[block[right]];
    if delta <= gap then FindChild := right-delta
    else FindChild := first[prevBlock[block[right]]] - delta + gap + 1
  end;
end;

procedure ReceiveAndDecode: integer;
var
  i, q: integer;
begin
  if M = n then q := n else q := Z; { Set q to the root node }
  while q > n do { Traverse down the tree }
    q := FindChild(q, Receive);
  if q = M then begin { Decode 0-node }
    q := 0;
    for i := 1 to E do q := q*2 + Receive;
    if q < R then q := q*2 + Receive else q := q + R;
    q := q + 1;
  end;
  ReceiveAndDecode := alpha[q];
end;

procedure InterchangeLeaves(e1, e2: integer);
var
  temp: integer;
begin
  rep[alpha[e1]] := e2; rep[alpha[e2]] := e1;
  temp := alpha[e1]; alpha[e1] := alpha[e2]; alpha[e2] := temp;
end;

procedure Update(k: integer);
var
  q, leafToIncrement, bq, b, oldParent, oldParity, nbq, par, bpar: integer;
  slide: boolean;

procedure FindNode;
begin
  q := rep[k]; leafToIncrement := 0;
  if q <= M then begin { A zero weight becomes positive }
    InterchangeLeaves(q, M);
    if R = 0 then begin R := M div 2; if R > 0 then E := E-1 end;
    M := M-1; R := R-1; q := M+1; bq := block[q];
    if M > 0 then begin
      { Split the 0-node into an internal node with two children.  The new
        0-node is node M; the old 0-node is node M+1; the new parent of
        nodes M and M+1 is node M+n }
      block[M] := bq; last[bq] := M; oldParent := parent[bq];
      parent[bq] := M+n; parity[bq] := 1;
      { Create a new internal block of zero weight for node M+n }
      b := availBlock; availBlock := nextBlock[availBlock];
      prevBlock[b] := bq; nextBlock[b] := nextBlock[bq];
      prevBlock[nextBlock[bq]] := b; nextBlock[bq] := b;
      parent[b] := oldParent; parity[b] := 0; rtChild[b] := q;
      block[M+n] := b; weight[b] := 0;
      first[b] := M+n; last[b] := M+n;
      leafToIncrement := q; q := M+n;
    end;
  end else begin { Interchange q with the first node in q's block }
    InterchangeLeaves(q, first[block[q]]);
    q := first[block[q]];
    if (q = M+1) and (M > 0) then begin
      leafToIncrement := q; q := parent[block[q]]
    end;
  end
end;

procedure SlideAndIncrement;
begin { q is currently the first node in its block }
  bq := block[q]; nbq := nextBlock[bq];
  par := parent[bq]; oldParent := par; oldParity := parity[bq];
  if (
       (q <= n) and (first[nbq] > n) and (weight[nbq] = weight[bq]))
     or
       ((q > n) and (first[nbq] <= n) and (weight[nbq] = weight[bq]+1)
     ) then begin     { Slide q over the next block }
    slide := true;
    oldParent := parent[nbq]; oldParity := parity[nbq];
    { Adjust child pointers for next higher level in tree }
    if par > 0 then begin
      bpar := block[par];
      if rtChild[bpar] = q then rtChild[bpar] := last[nbq]
      else if rtChild[bpar] = first[nbq] then rtChild[bpar] := q
      else rtChild[bpar] := rtChild[bpar]+1;
      if par <> Z then
        if block[par+1] <> bpar then
          if rtChild[block[par+1]] = first[nbq] then
            rtChild[block[par+1]] := q
          else if block[rtChild[block[par+1]]] = nbq then
            rtChild[block[par+1]] := rtChild[block[par+1]]+1
    end;
    { Adjust parent pointers for block nbq }
    parent[nbq] := parent[nbq] - 1 + parity[nbq]; parity[nbq] := 1-parity[nbq];
    nbq := nextBlock[nbq];
  end else slide := false;
  if (
      ((q <= n) and (first[nbq] <= n))
      or
      ((q > n) and (first[nbq] > n))
     )
     and (weight[nbq] = weight[bq]+1) then begin
    { Merge q into the block of weight one higher }
    block[q] := nbq; last[nbq] := q;
    if last[bq] = q then begin { q's old block disappears }
      nextblock[prevBlock[bq]] := nextBlock[bq];
      prevBlock[nextBlock[bq]] := prevBlock[bq];
      nextBlock[bq] := availBlock; availBlock := bq;
    end else begin
      if q > n then rtChild[bq] := FindChild(q-1, 1);
      if parity[bq] = 0 then parent[bq] := parent[bq] - 1;
      parity[bq] := 1-parity[bq];
      first[bq] := q-1;
    end
  end else if last[bq] = q then begin
    if slide then begin { q's block is slid forward in the block list }
      prevBlock[nextBlock[bq]] := prevBlock[bq];
      nextBlock[prevBlock[bq]] := nextBlock[bq];
      prevBlock[bq] := prevBlock[nbq]; nextBlock[bq] := nbq;
      prevBlock[nbq] := bq; nextBlock[prevBlock[bq]] := bq;
      parent[bq] := oldParent; parity[bq] := oldParity;
    end;
    weight[bq] := weight[bq]+1;
  end else begin
    { A new block is created for q }
    b := availBlock; availBlock := nextBlock[availBlock];
    block[q] := b; first[b] := q; last[b] := q;
    if q > n then begin
      rtChild[b] := rtChild[bq];
      rtChild[bq] := FindChild(q-1, 1);
      if rtChild[b] = q-1 then parent[bq] := q
      else if parity[bq] = 0 then parent[bq] := parent[bq]-1
    end else if parity[bq] = 0 then parent[bq] := parent[bq]-1;
    first[bq] := q-1; parity[bq] := 1-parity[bq];
    { Insert q's in block in its proper place in the block list }
    prevBlock[b] := prevBlock[nbq]; nextBlock[b] := nbq;
    prevBlock[nbq] := b; nextBlock[prevBlock[b]] := b;
    weight[b] := weight[bq]+1;
    parent[b] := oldParent; parity[b] := oldParity
  end;
  { Move q one level higher in the tree }
  if q <= n then q := oldParent else q := par
end;

begin
  { Set q to the node whose weight should increase }
  FindNode;
  while q > 0 do
    { At this point, q is the first node in its block.  Increment q's weight
      by 1 and slide q if necessary over the next block to maintain the
      invariant.  Then set q to the node one level higher that needs
      incrementing next }
    SlideAndIncrement;
  { Finish up some special cases involving the 0-node }
  if LeafToIncrement <> 0 then begin
    q := leafToIncrement;
    SlideAndIncrement
  end
end;

begin
  outpos := 0;
  Initialize;
  { Decide to encode or decode depending on data in file! }
  { This is only for pedagogical purposes of course.      }
  if (input^ <> '0') and (input^ <> '1') then begin
    while (not eof) do begin
      while (not eoln) do begin
        read(c);
        EncodeAndTransmit(ORD(c));
        Update(ORD(c));
      end;
      readln;
      EncodeAndTransmit(10);
      Update(10);
    end;
  end else begin
    while (not eof) do begin
      a := ReceiveAndDecode;
      if a = 10 then writeln else write(CHR(a));
      Update(a);
    end;
    writeln;
  end
end.