libcds/src/coders/huff.cpp

   1 /* huff.cpp
   2    Copyright (C) 2008, Gonzalo Navarro, all rights reserved.
   3
   4    Canonical Huffman
   5
   6    This library is free software; you can redistribute it and/or
   7    modify it under the terms of the GNU Lesser General Public
   8    License as published by the Free Software Foundation; either
   9    version 2.1 of the License, or (at your option) any later version.
  10
  11    This library is distributed in the hope that it will be useful,
  12    but WITHOUT ANY WARRANTY; without even the implied warranty of
  13    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14    Lesser General Public License for more details.
  15
  16    You should have received a copy of the GNU Lesser General Public
  17    License along with this library; if not, write to the Free Software
  18    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
  19
  20 */
  21 // implements canonical Huffman
  22
  23 #include <huff.h>
  24
  25 typedef struct
  26    { uint freq;
  27      uint symb;
  28      union
  29        { int prev;
  30          uint depth;
  31        } h;
  32      int ch1,ch2;
  33    } Ttree;
  34
  35 static void sort (Ttree *tree, int lo, int up)
  36
  37    { uint i, j;
  38      Ttree temp;
  39      while (up>lo)
  40         { i = lo;
  41           j = up;
  42           temp = tree[lo];
  43           while (i<j)
  44              { while (tree[j].freq > temp.freq) j--;
  45                tree[i] = tree[j];
  46                while (i<j && tree[i].freq <= temp.freq) i++;
  47                tree[j] = tree[i];
  48              }
  49           tree[i] = temp;
  50           if (i-lo < up-i) { sort(tree,lo,i-1); lo = i+1; }
  51           else { sort(tree,i+1,up); up = i-1; }
  52         }
  53    }
  54
  55 static void setdepths (Ttree *tree, uint node, int depth)
  56
  57    { if (tree[node].ch1 == -1) // leaf
  58         { tree[node].h.depth = depth;
  59           return;
  60         }
  61      setdepths (tree,tree[node].ch1,depth+1);
  62      setdepths (tree,tree[node].ch2,depth+1);
  63    }
  64
  65 THuff createHuff (uint *freq, uint lim)
  66
  67    { THuff H;
  68      int i,j,d;
  69      Ttree *tree;
  70      uint ptr,last,fre;
  71          // remove zero frequencies
  72      H.max = lim;
  73      tree = (Ttree*)malloc((2*(lim+1)-1)*sizeof(Ttree));
  74      j = 0;
  75      for (i=0;i<=(int)lim;i++)
  76         { if (freq[i]>0)
  77              { tree[j].freq = freq[i];
  78                tree[j].symb = i;
  79                j++;
  80              }
  81         }
  82      H.lim = lim = j-1;
  83         // now run Huffman algorithm
  84      sort (tree,0,lim);
  85      for (i=0;i<=(int)lim;i++)
  86          { tree[i].h.prev = i+1;
  87            tree[i].ch1 = tree[i].ch2 = -1;
  88          }
  89      tree[lim].h.prev = -1;
  90         // last = next node to process, ptr = search point, fre = next free cell
  91         // leaves are in 0..lim in decreasing freq order
  92         // internal nodes are in lim+1.. 2*lim, created in incr. fre order
  93      last=0; ptr = 0; fre = lim+1;
  94      for (i=0;i<(int)lim;i++)
  95         { tree[fre].ch1 = last;
  96           last = tree[last].h.prev;
  97           tree[fre].ch2 = last;
  98           tree[fre].freq = tree[tree[fre].ch1].freq+tree[tree[fre].ch2].freq;
  99           while ((tree[ptr].h.prev != -1) &&
 100                  (tree[tree[ptr].h.prev].freq <= tree[fre].freq))
 101               ptr = tree[ptr].h.prev;
 102           tree[fre].h.prev = tree[ptr].h.prev;
 103           tree[ptr].h.prev = fre;
 104           last = tree[last].h.prev;
 105           fre++;
 106         }
 107         // now assign depths recursively
 108      setdepths (tree,2*lim,0);
 109      H.s.spos = (uint*)malloc ((H.max+1)*sizeof(uint));
 110      for (i=0;i<=(int)H.max;i++) H.s.spos[i] = ~0;
 111      H.num = (uint*)malloc ((lim+1)*sizeof(uint)); // max possible depth
 112      d=0;
 113      for (i=lim;i>=0;i--)
 114         { H.s.spos[tree[i].symb] = i;
 115           while ((int)tree[i].h.depth > d)
 116               { H.num[d] = i+1; d++; }
 117         }
 118      H.num[d] = 0;
 119      H.depth = d;
 120      for (d=H.depth;d>0;d--) H.num[d] = H.num[d-1] - H.num[d];
 121      H.num[0] = (lim == 0);
 122      H.num = (uint*)realloc(H.num,(H.depth+1)*sizeof(uint));
 123      H.total = 0;
 124      for (i=0;i<=(int)lim;i++)
 125        H.total += freq[tree[i].symb] * tree[i].h.depth;
 126      free (tree);
 127      return H;
 128    }
 129
 130 void bitzero (register uint *e, register uint p,
 131          register uint len)
 132
 133    { e += p/W; p %= W;
 134      if (p+len >= W)
 135   { *e &= ~((1<<p)-1);
 136     len -= p;
 137     e++; p = 0;
 138   }
 139      while (len >= W)
 140   { *e++ = 0;
 141     len -= W;
 142   }
 143      if (len > 0)
 144   *e &= ~(((1<<len)-1)<<p);
 145    }
 146
 147 ulong encodeHuff (THuff H, uint symb, uint *stream, ulong ptr)
 148
 149    { uint pos;
 150      uint code;
 151      uint d;
 152      pos = H.s.spos[symb];
 153      code = 0;
 154      d = H.depth;
 155      while (pos >= H.num[d])
 156        { code = (code + H.num[d]) >> 1;
 157          pos -= H.num[d--];
 158        }
 159      code += pos;
 160      if (d > W) { bitzero(stream,ptr,d-W); ptr += d-W; d = W; }
 161      while (d--)
 162         { if ((code >> d) & 1) bitset(stream,ptr);
 163           else bitclean(stream,ptr);
 164           ptr++;
 165         }
 166      return ptr;
 167    }
 168
 169 ulong decodeHuff (THuff H, uint *symb, uint *stream, ulong ptr)
 170
 171    { uint pos;
 172      uint d;
 173      pos = 0;
 174      d = 0;
 175      while (pos < H.fst[d])
 176         { pos = (pos << 1) | bitget(stream,ptr);
 177           ptr++; d++;
 178         }
 179      *symb = H.s.symb[H.num[d]+pos-H.fst[d]];
 180      return ptr;
 181    }
 182
 183 /*   { uint pos;  // This "improved" code is actually slower!
 184      int d;
 185      uint wrd,off;
 186      stream += ptr/W;
 187      off = ptr & (W-1);
 188      wrd = *stream >> off;
 189      pos = 0;
 190      d = 0;
 191      while (pos < H.fst[d])
 192        { pos = (pos << 1) | (wrd & 1);
 193          d++; wrd >>= 1; off++;
 194          if (off == W) { wrd = *++stream; off = 0; }
 195        }
 196      *symb = H.s.symb[H.num[d]+pos-H.fst[d]];
 197      return ptr+d;
 198    }
 199 */
 200 void saveHuff (THuff H, FILE *f)
 201
 202    { uint *symb = new uint[H.lim+1];
 203      uint i;
 204                  for(i=0;i<(H.lim+1);i++) symb[i] = 0;
 205      for (i=0;i<=H.max;i++)
 206                          if (H.s.spos[i] != (uint)~0) symb[H.s.spos[i]] = i;
 207      uint l=fwrite (&H.max,sizeof(uint),1,f);
 208      l += fwrite (&H.lim,sizeof(uint),1,f);
 209      l += fwrite (&H.depth,sizeof(uint),1,f);
 210      l += fwrite (symb,sizeof(uint),H.lim+1,f);
 211      l += fwrite (H.num,sizeof(uint),H.depth+1,f);
 212      delete [] (symb);
 213    }
 214
 215 uint sizeHuff (THuff H)
 216
 217    { return (4+(H.lim+1)+(H.depth+1))*sizeof(uint);
 218    }
 219
 220 void freeHuff (THuff H)
 221
 222   { free (H.s.spos); free (H.num);
 223   }
 224
 225 THuff loadHuff (FILE *f, int enc)
 226
 227    { THuff H;
 228      uint *symb;
 229      //uint *num;
 230      uint i,d,dold,dact;
 231      uint l = fread (&H.max,sizeof(uint),1,f);
 232      l += fread (&H.lim,sizeof(uint),1,f);
 233      l += fread (&H.depth,sizeof(uint),1,f);
 234      symb = (uint*)malloc ((H.lim+1)*sizeof(uint));
 235      l += fread (symb,sizeof(uint),H.lim+1,f);
 236      if (enc)
 237           { H.s.spos = (uint*)malloc ((H.max+1)*sizeof(uint));
 238             for (i=0;i<=H.max;i++) H.s.spos[i] = (uint)~0;
 239             for (i=0;i<=H.lim;i++) H.s.spos[symb[i]] = i;
 240             free (symb);
 241           }
 242      else H.s.symb = symb;
 243      H.num = (uint*)malloc ((H.depth+1)*sizeof(uint));
 244      l += fread (H.num,sizeof(uint),H.depth+1,f);
 245      if (!enc)
 246           { H.fst = (uint*)malloc ((H.depth+1)*sizeof(uint));
 247             H.fst[H.depth] = 0; dold = 0;
 248             for (d=H.depth-1;d>=0;d--)
 249                 { dact = H.num[d+1];
 250                   H.fst[d] = (H.fst[d+1]+dact) >> 1;
 251                   H.num[d+1] = dold;
 252                   dold += dact;
 253                 }
 254             H.num[0] = dold;
 255           }
 256      return H;
 257    }