Import libcds.
[SXSI/libcds.git] / src / coders / huff.cpp
diff --git a/src/coders/huff.cpp b/src/coders/huff.cpp
new file mode 100644 (file)
index 0000000..52b95f9
--- /dev/null
@@ -0,0 +1,257 @@
+/* huff.cpp
+   Copyright (C) 2008, Gonzalo Navarro, all rights reserved.
+
+   Canonical Huffman
+
+   This library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   This library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with this library; if not, write to the Free Software
+   Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+
+*/
+// implements canonical Huffman 
+
+#include <huff.h>
+
+typedef struct
+   { uint freq;
+     uint symb;
+     union
+       { int prev;
+        uint depth;
+       } h;
+     int ch1,ch2;
+   } Ttree;
+  
+static void sort (Ttree *tree, int lo, int up)
+
+   { uint i, j;
+     Ttree temp;
+     while (up>lo) 
+        { i = lo;
+         j = up;
+         temp = tree[lo]; 
+         while (i<j) 
+            { while (tree[j].freq > temp.freq) j--;
+              tree[i] = tree[j]; 
+              while (i<j && tree[i].freq <= temp.freq) i++;
+              tree[j] = tree[i];
+            }
+         tree[i] = temp;
+         if (i-lo < up-i) { sort(tree,lo,i-1); lo = i+1; }
+         else { sort(tree,i+1,up); up = i-1; }
+       }
+   }
+
+static void setdepths (Ttree *tree, uint node, int depth)
+
+   { if (tree[node].ch1 == -1) // leaf
+       { tree[node].h.depth = depth; 
+         return;
+       }
+     setdepths (tree,tree[node].ch1,depth+1);
+     setdepths (tree,tree[node].ch2,depth+1);
+   }
+
+THuff createHuff (uint *freq, uint lim)
+
+   { THuff H;
+     int i,j,d;
+     Ttree *tree;
+     uint ptr,last,fre;
+         // remove zero frequencies
+     H.max = lim;
+     tree = (Ttree*)malloc((2*(lim+1)-1)*sizeof(Ttree));
+     j = 0;
+     for (i=0;i<=(int)lim;i++) 
+       { if (freq[i]>0) 
+            { tree[j].freq = freq[i];   
+               tree[j].symb = i;   
+              j++;
+            }
+       }
+     H.lim = lim = j-1;
+       // now run Huffman algorithm
+     sort (tree,0,lim);
+     for (i=0;i<=(int)lim;i++) 
+         { tree[i].h.prev = i+1;
+          tree[i].ch1 = tree[i].ch2 = -1;
+        }
+     tree[lim].h.prev = -1;
+        // last = next node to process, ptr = search point, fre = next free cell
+       // leaves are in 0..lim in decreasing freq order
+       // internal nodes are in lim+1.. 2*lim, created in incr. fre order
+     last=0; ptr = 0; fre = lim+1;
+     for (i=0;i<(int)lim;i++)
+       { tree[fre].ch1 = last;
+         last = tree[last].h.prev;
+         tree[fre].ch2 = last;
+         tree[fre].freq = tree[tree[fre].ch1].freq+tree[tree[fre].ch2].freq;
+         while ((tree[ptr].h.prev != -1) &&
+                (tree[tree[ptr].h.prev].freq <= tree[fre].freq))
+             ptr = tree[ptr].h.prev;
+         tree[fre].h.prev = tree[ptr].h.prev;
+         tree[ptr].h.prev = fre;
+         last = tree[last].h.prev;
+         fre++;
+       }
+        // now assign depths recursively
+     setdepths (tree,2*lim,0);
+     H.s.spos = (uint*)malloc ((H.max+1)*sizeof(uint));
+     for (i=0;i<=(int)H.max;i++) H.s.spos[i] = ~0;
+     H.num = (uint*)malloc ((lim+1)*sizeof(uint)); // max possible depth
+     d=0;
+     for (i=lim;i>=0;i--)
+       { H.s.spos[tree[i].symb] = i;
+         while ((int)tree[i].h.depth > d) 
+             { H.num[d] = i+1; d++; }
+       }
+     H.num[d] = 0;
+     H.depth = d;
+     for (d=H.depth;d>0;d--) H.num[d] = H.num[d-1] - H.num[d];
+     H.num[0] = (lim == 0);
+     H.num = (uint*)realloc(H.num,(H.depth+1)*sizeof(uint));
+     H.total = 0;
+     for (i=0;i<=(int)lim;i++) 
+       H.total += freq[tree[i].symb] * tree[i].h.depth;
+     free (tree);
+     return H;
+   }
+
+void bitzero (register uint *e, register uint p, 
+         register uint len)
+
+   { e += p/W; p %= W;
+     if (p+len >= W)
+  { *e &= ~((1<<p)-1);
+    len -= p;
+    e++; p = 0;
+  }
+     while (len >= W)
+  { *e++ = 0;
+    len -= W;
+  }
+     if (len > 0)
+  *e &= ~(((1<<len)-1)<<p);
+   }
+   
+ulong encodeHuff (THuff H, uint symb, uint *stream, ulong ptr)
+
+   { uint pos;
+     uint code;
+     uint d;
+     pos = H.s.spos[symb];
+     code = 0;
+     d = H.depth; 
+     while (pos >= H.num[d])
+       { code = (code + H.num[d]) >> 1;
+         pos -= H.num[d--];
+       }
+     code += pos;
+     if (d > W) { bitzero(stream,ptr,d-W); ptr += d-W; d = W; }
+     while (d--)
+        { if ((code >> d) & 1) bitset(stream,ptr);
+         else bitclean(stream,ptr);
+         ptr++;
+       }
+     return ptr;
+   } 
+
+ulong decodeHuff (THuff H, uint *symb, uint *stream, ulong ptr)
+
+   { uint pos;
+     uint d;
+     pos = 0;
+     d = 0;
+     while (pos < H.fst[d])
+        { pos = (pos << 1) | bitget(stream,ptr);
+          ptr++; d++;
+        }
+     *symb = H.s.symb[H.num[d]+pos-H.fst[d]];
+     return ptr;
+   }
+
+/*   { uint pos;  // This "improved" code is actually slower!
+     int d;
+     uint wrd,off;
+     stream += ptr/W;
+     off = ptr & (W-1);
+     wrd = *stream >> off;
+     pos = 0;
+     d = 0;
+     while (pos < H.fst[d])
+       { pos = (pos << 1) | (wrd & 1);
+        d++; wrd >>= 1; off++;
+        if (off == W) { wrd = *++stream; off = 0; }
+       }
+     *symb = H.s.symb[H.num[d]+pos-H.fst[d]];
+     return ptr+d;
+   } 
+*/
+void saveHuff (THuff H, FILE *f)
+
+   { uint *symb = new uint[H.lim+1];
+     uint i;
+                for(i=0;i<(H.lim+1);i++) symb[i] = 0;
+     for (i=0;i<=H.max;i++) 
+                        if (H.s.spos[i] != (uint)~0) symb[H.s.spos[i]] = i;
+     uint l=fwrite (&H.max,sizeof(uint),1,f); 
+     l += fwrite (&H.lim,sizeof(uint),1,f); 
+     l += fwrite (&H.depth,sizeof(uint),1,f); 
+     l += fwrite (symb,sizeof(uint),H.lim+1,f); 
+     l += fwrite (H.num,sizeof(uint),H.depth+1,f); 
+     delete [] (symb);
+   }
+
+uint sizeHuff (THuff H)
+
+   { return (4+(H.lim+1)+(H.depth+1))*sizeof(uint);
+   }
+
+void freeHuff (THuff H)
+
+  { free (H.s.spos); free (H.num);
+  }
+
+THuff loadHuff (FILE *f, int enc)
+
+   { THuff H;
+     uint *symb;
+     //uint *num;
+     uint i,d,dold,dact;
+     uint l = fread (&H.max,sizeof(uint),1,f); 
+     l += fread (&H.lim,sizeof(uint),1,f); 
+     l += fread (&H.depth,sizeof(uint),1,f); 
+     symb = (uint*)malloc ((H.lim+1)*sizeof(uint));
+     l += fread (symb,sizeof(uint),H.lim+1,f); 
+     if (enc) 
+          { H.s.spos = (uint*)malloc ((H.max+1)*sizeof(uint));
+            for (i=0;i<=H.max;i++) H.s.spos[i] = (uint)~0;
+            for (i=0;i<=H.lim;i++) H.s.spos[symb[i]] = i;
+            free (symb);
+         }
+     else H.s.symb = symb;
+     H.num = (uint*)malloc ((H.depth+1)*sizeof(uint));
+     l += fread (H.num,sizeof(uint),H.depth+1,f); 
+     if (!enc) 
+          { H.fst = (uint*)malloc ((H.depth+1)*sizeof(uint));
+            H.fst[H.depth] = 0; dold = 0;
+            for (d=H.depth-1;d>=0;d--)
+               { dact = H.num[d+1];
+                 H.fst[d] = (H.fst[d+1]+dact) >> 1;
+                 H.num[d+1] = dold;
+                 dold += dact;
+               }
+           H.num[0] = dold;
+         }
+     return H;
+   }