libcds/src/coders/huff.cpp

   1
   2 // implements canonical Huffman
   3
   4 #include <huff.h>
   5
   6 typedef struct
   7    { uint freq;
   8      uint symb;
   9      union
  10        { int prev;
  11          uint depth;
  12        } h;
  13      int ch1,ch2;
  14    } Ttree;
  15
  16 static void sort (Ttree *tree, int lo, int up)
  17
  18    { uint i, j;
  19      Ttree temp;
  20      while (up>lo)
  21         { i = lo;
  22           j = up;
  23           temp = tree[lo];
  24           while (i<j)
  25              { while (tree[j].freq > temp.freq) j--;
  26                tree[i] = tree[j];
  27                while (i<j && tree[i].freq <= temp.freq) i++;
  28                tree[j] = tree[i];
  29              }
  30           tree[i] = temp;
  31           if (i-lo < up-i) { sort(tree,lo,i-1); lo = i+1; }
  32           else { sort(tree,i+1,up); up = i-1; }
  33         }
  34    }
  35
  36 static void setdepths (Ttree *tree, uint node, int depth)
  37
  38    { if (tree[node].ch1 == -1) // leaf
  39         { tree[node].h.depth = depth;
  40           return;
  41         }
  42      setdepths (tree,tree[node].ch1,depth+1);
  43      setdepths (tree,tree[node].ch2,depth+1);
  44    }
  45
  46 THuff createHuff (uint *freq, uint lim)
  47
  48    { THuff H;
  49      int i,j,d;
  50      Ttree *tree;
  51      uint ptr,last,fre;
  52          // remove zero frequencies
  53      H.max = lim;
  54      tree = (Ttree*)malloc((2*(lim+1)-1)*sizeof(Ttree));
  55      j = 0;
  56      for (i=0;i<=(int)lim;i++)
  57         { if (freq[i]>0)
  58              { tree[j].freq = freq[i];
  59                tree[j].symb = i;
  60                j++;
  61              }
  62         }
  63      H.lim = lim = j-1;
  64         // now run Huffman algorithm
  65      sort (tree,0,lim);
  66      for (i=0;i<=(int)lim;i++)
  67          { tree[i].h.prev = i+1;
  68            tree[i].ch1 = tree[i].ch2 = -1;
  69          }
  70      tree[lim].h.prev = -1;
  71         // last = next node to process, ptr = search point, fre = next free cell
  72         // leaves are in 0..lim in decreasing freq order
  73         // internal nodes are in lim+1.. 2*lim, created in incr. fre order
  74      last=0; ptr = 0; fre = lim+1;
  75      for (i=0;i<(int)lim;i++)
  76         { tree[fre].ch1 = last;
  77           last = tree[last].h.prev;
  78           tree[fre].ch2 = last;
  79           tree[fre].freq = tree[tree[fre].ch1].freq+tree[tree[fre].ch2].freq;
  80           while ((tree[ptr].h.prev != -1) &&
  81                  (tree[tree[ptr].h.prev].freq <= tree[fre].freq))
  82               ptr = tree[ptr].h.prev;
  83           tree[fre].h.prev = tree[ptr].h.prev;
  84           tree[ptr].h.prev = fre;
  85           last = tree[last].h.prev;
  86           fre++;
  87         }
  88         // now assign depths recursively
  89      setdepths (tree,2*lim,0);
  90      H.s.spos = (uint*)malloc ((H.max+1)*sizeof(uint));
  91      for (i=0;i<=(int)H.max;i++) H.s.spos[i] = ~0;
  92      H.num = (uint*)malloc ((lim+1)*sizeof(uint)); // max possible depth
  93      d=0;
  94      for (i=lim;i>=0;i--)
  95         { H.s.spos[tree[i].symb] = i;
  96           while ((int)tree[i].h.depth > d)
  97               { H.num[d] = i+1; d++; }
  98         }
  99      H.num[d] = 0;
 100      H.depth = d;
 101      for (d=H.depth;d>0;d--) H.num[d] = H.num[d-1] - H.num[d];
 102      H.num[0] = (lim == 0);
 103      H.num = (uint*)realloc(H.num,(H.depth+1)*sizeof(uint));
 104      H.total = 0;
 105      for (i=0;i<=(int)lim;i++)
 106        H.total += freq[tree[i].symb] * tree[i].h.depth;
 107      free (tree);
 108      return H;
 109    }
 110
 111 void bitzero (register uint *e, register uint p,
 112          register uint len)
 113
 114    { e += p/W; p %= W;
 115      if (p+len >= W)
 116   { *e &= ~((1<<p)-1);
 117     len -= p;
 118     e++; p = 0;
 119   }
 120      while (len >= W)
 121   { *e++ = 0;
 122     len -= W;
 123   }
 124      if (len > 0)
 125   *e &= ~(((1<<len)-1)<<p);
 126    }
 127
 128 ulong encodeHuff (THuff H, uint symb, uint *stream, ulong ptr)
 129
 130    { uint pos;
 131      uint code;
 132      uint d;
 133      pos = H.s.spos[symb];
 134      code = 0;
 135      d = H.depth;
 136      while (pos >= H.num[d])
 137        { code = (code + H.num[d]) >> 1;
 138          pos -= H.num[d--];
 139        }
 140      code += pos;
 141      if (d > W) { bitzero(stream,ptr,d-W); ptr += d-W; d = W; }
 142      while (d--)
 143         { if ((code >> d) & 1) bitset(stream,ptr);
 144           else bitclean(stream,ptr);
 145           ptr++;
 146         }
 147      return ptr;
 148    }
 149
 150 ulong decodeHuff (THuff H, uint *symb, uint *stream, ulong ptr)
 151
 152    { uint pos;
 153      uint d;
 154      pos = 0;
 155      d = 0;
 156      while (pos < H.fst[d])
 157         { pos = (pos << 1) | bitget(stream,ptr);
 158           ptr++; d++;
 159         }
 160      *symb = H.s.symb[H.num[d]+pos-H.fst[d]];
 161      return ptr;
 162    }
 163
 164 /*   { uint pos;  // This "improved" code is actually slower!
 165      int d;
 166      uint wrd,off;
 167      stream += ptr/W;
 168      off = ptr & (W-1);
 169      wrd = *stream >> off;
 170      pos = 0;
 171      d = 0;
 172      while (pos < H.fst[d])
 173        { pos = (pos << 1) | (wrd & 1);
 174          d++; wrd >>= 1; off++;
 175          if (off == W) { wrd = *++stream; off = 0; }
 176        }
 177      *symb = H.s.symb[H.num[d]+pos-H.fst[d]];
 178      return ptr+d;
 179    }
 180 */
 181 void saveHuff (THuff H, FILE *f)
 182
 183    { uint *symb = (uint*)malloc((H.lim+1)*sizeof(uint));
 184      uint i;
 185      for (i=0;i<=H.max;i++)
 186          if (H.s.spos[i] != (uint)~0) symb[H.s.spos[i]] = i;
 187      uint l=fwrite (&H.max,sizeof(uint),1,f);
 188      l += fwrite (&H.lim,sizeof(uint),1,f);
 189      l += fwrite (&H.depth,sizeof(uint),1,f);
 190      l += fwrite (symb,sizeof(uint),H.lim+1,f);
 191      l += fwrite (H.num,sizeof(uint),H.depth+1,f);
 192      free (symb);
 193    }
 194
 195 uint sizeHuff (THuff H)
 196
 197    { return (4+(H.lim+1)+(H.depth+1))*sizeof(uint);
 198    }
 199
 200 void freeHuff (THuff H)
 201
 202   { free (H.s.spos); free (H.num);
 203   }
 204
 205 THuff loadHuff (FILE *f, int enc)
 206
 207    { THuff H;
 208      uint *symb;
 209      //uint *num;
 210      uint i,d,dold,dact;
 211      uint l = fread (&H.max,sizeof(uint),1,f);
 212      l += fread (&H.lim,sizeof(uint),1,f);
 213      l += fread (&H.depth,sizeof(uint),1,f);
 214      symb = (uint*)malloc ((H.lim+1)*sizeof(uint));
 215      l += fread (symb,sizeof(uint),H.lim+1,f);
 216      if (enc)
 217           { H.s.spos = (uint*)malloc ((H.max+1)*sizeof(uint));
 218             for (i=0;i<=H.max;i++) H.s.spos[i] = (uint)~0;
 219             for (i=0;i<=H.lim;i++) H.s.spos[symb[i]] = i;
 220             free (symb);
 221           }
 222      else H.s.symb = symb;
 223      H.num = (uint*)malloc ((H.depth+1)*sizeof(uint));
 224      l += fread (H.num,sizeof(uint),H.depth+1,f);
 225      if (!enc)
 226           { H.fst = (uint*)malloc ((H.depth+1)*sizeof(uint));
 227             H.fst[H.depth] = 0; dold = 0;
 228             for (d=H.depth-1;d>=0;d--)
 229                 { dact = H.num[d+1];
 230                   H.fst[d] = (H.fst[d+1]+dact) >> 1;
 231                   H.num[d+1] = dold;
 232                   dold += dact;
 233                 }
 234             H.num[0] = dold;
 235           }
 236      return H;
 237    }
 238
 239