2 Copyright (C) 2008, Gonzalo Navarro, all rights reserved.
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.
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.
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
21 // implements canonical Huffman
35 static void sort (Ttree *tree, int lo, int up)
44 { while (tree[j].freq > temp.freq) j--;
46 while (i<j && tree[i].freq <= temp.freq) i++;
50 if (i-lo < up-i) { sort(tree,lo,i-1); lo = i+1; }
51 else { sort(tree,i+1,up); up = i-1; }
55 static void setdepths (Ttree *tree, uint node, int depth)
57 { if (tree[node].ch1 == -1) // leaf
58 { tree[node].h.depth = depth;
61 setdepths (tree,tree[node].ch1,depth+1);
62 setdepths (tree,tree[node].ch2,depth+1);
65 THuff createHuff (uint *freq, uint lim)
71 // remove zero frequencies
73 tree = (Ttree*)malloc((2*(lim+1)-1)*sizeof(Ttree));
75 for (i=0;i<=(int)lim;i++)
77 { tree[j].freq = freq[i];
83 // now run Huffman algorithm
85 for (i=0;i<=(int)lim;i++)
86 { tree[i].h.prev = i+1;
87 tree[i].ch1 = tree[i].ch2 = -1;
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;
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;
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
114 { H.s.spos[tree[i].symb] = i;
115 while ((int)tree[i].h.depth > d)
116 { H.num[d] = i+1; 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));
124 for (i=0;i<=(int)lim;i++)
125 H.total += freq[tree[i].symb] * tree[i].h.depth;
130 void bitzero (register uint *e, register uint p,
144 *e &= ~(((1<<len)-1)<<p);
147 ulong encodeHuff (THuff H, uint symb, uint *stream, ulong ptr)
152 pos = H.s.spos[symb];
155 while (pos >= H.num[d])
156 { code = (code + H.num[d]) >> 1;
160 if (d > W) { bitzero(stream,ptr,d-W); ptr += d-W; d = W; }
162 { if ((code >> d) & 1) bitset(stream,ptr);
163 else bitclean(stream,ptr);
169 ulong decodeHuff (THuff H, uint *symb, uint *stream, ulong ptr)
175 while (pos < H.fst[d])
176 { pos = (pos << 1) | bitget(stream,ptr);
179 *symb = H.s.symb[H.num[d]+pos-H.fst[d]];
183 /* { uint pos; // This "improved" code is actually slower!
188 wrd = *stream >> off;
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; }
196 *symb = H.s.symb[H.num[d]+pos-H.fst[d]];
200 void saveHuff (THuff H, FILE *f)
202 { uint *symb = new uint[H.lim+1];
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);
215 uint sizeHuff (THuff H)
217 { return (4+(H.lim+1)+(H.depth+1))*sizeof(uint);
220 void freeHuff (THuff H)
222 { free (H.s.spos); free (H.num);
225 THuff loadHuff (FILE *f, int enc)
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);
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;
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);
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--)
250 H.fst[d] = (H.fst[d+1]+dact) >> 1;