|
|
1.1 root 1: /* @(#)huff.c 1.3 */
2: #include <stdio.h>
3: #define BYTE 8
4: #define QW 1 /* width of bas-q digit in bits */
5:
6: /* this stuff should be local and hidden; it was made
7: * accessible outside for dirty reasons: 20% faster spell
8: */
9: #include "huff.h"
10: struct huff huffcode;
11:
12: /* Infinite Huffman code
13: *
14: * Let the messages be exponentially distributed with ratio r:
15: * P{message k} = r^k*(1-r), k=0,1,...
16: * Let the messages be coded in base q, and suppose
17: * r^n = 1/q
18: * If each decade (base q) contains n codes, then
19: * the messages assigned to each decade will be q times
20: * as probable as the next. Moreover the code for the tail of
21: * the distribution after truncating one decade should look
22: * just like the original, but longer by one leading digit q-1.
23: * q(z+n) = z + (q-1)q^w
24: * where z is first code of decade, w is width of code (in shortest
25: * full decade). Examples, base 2:
26: * r^1 = 1/2 r^5 = 1/2
27: * 0 0110
28: * 10 0111
29: * 110 1000
30: * 1110 1001
31: * ... 1010
32: * 10110
33: * w=1,z=0 w=4,z=0110
34: * Rewriting slightly
35: * (q-1)z + q*n = (q-1)q^w
36: * whence z is a multiple of q and n is a multiple of q-1. Let
37: * z = cq, n = d(q-1)
38: * We pick w to be the least integer such that
39: * d = n/(q-1) <= q^(w-1)
40: * Then solve for c
41: * c = q^(w-1) - d
42: * If c is not zero, the first decade may be preceded by
43: * even shorter (w-1)-digit codes 0,1,...,c-1. Thus
44: * the example code with r^5 = 1/2 becomes
45: * 000
46: * 001
47: * 010
48: * 0110
49: * 0111
50: * 1000
51: * 1001
52: * 1010
53: * 10110
54: * ...
55: * 110110
56: * ...
57: * The expected number of base-q digits in a codeword is then
58: * w - 1 + r^c/(1-r^n)
59: * The present routines require q to be a power of 2
60: */
61: /* There is a lot of hanky-panky with left justification against
62: * sign instead of simple left justification because
63: * unsigned long is not available
64: */
65: #define L (BYTE*(sizeof(long))-1) /* length of signless long */
66: #define MASK (~(1L<<L)) /* mask out sign */
67:
68: /* decode the prefix of word y (which is left justified against sign)
69: * place mesage number into place pointed to by kp
70: * return length (in bits) of decoded prefix or 0 if code is out of
71: * range
72: */
73: decode(y,pk)
74: long y;
75: long *pk;
76: {
77: register l;
78: long v;
79: if(y < cs) {
80: *pk = y >> (L+QW-w);
81: return(w-QW);
82: }
83: for(l=w,v=v0; y>=qcs; y=(y<<QW)&MASK,v+=n)
84: if((l+=QW) > L)
85: return(0);
86: *pk = v + (y>>(L-w));
87: return(l);
88: }
89:
90: /* encode message k and put result (right justified) into
91: * place pointed to by py.
92: * return length (in bits) of result,
93: * or 0 if code is too long
94: */
95: encode(k,py)
96: long k;
97: long *py;
98: {
99: register l;
100: long y;
101: if(k < c) {
102: *py = k;
103: return(w-QW);
104: }
105: for(k-=c,y=1,l=w; k>=n; k-=n,y<<=QW)
106: if((l+=QW) > L)
107: return(0);
108: *py = ((y-1)<<w) + cq + k;
109: return(l);
110: }
111: /* Initialization code, given expected value of k
112: * E(k) = r/(1-r) = a
113: * and given base width b
114: * return expected length of coded messages
115: */
116: struct qlog {
117: long p;
118: double u;
119: };
120: double
121: huff(a)
122: float a;
123: {
124: struct qlog qlog();
125: struct qlog z;
126: register i, q;
127: long d, j;
128: double r = a/(1.0 + a);
129: double rc, rq;
130: for(i=0,q=1,rq=r; i<QW; i++,q*=2,rq*=rq)
131: continue;
132: rq /= r; /* rq = r^(q-1) */
133: z = qlog(rq, 1./q, 1L, rq);
134: d = z.p;
135: n = d*(q-1);
136: if(n!=d*(q-1))
137: abort(); /*time to make n long*/
138: for(w=QW,j=1; j<d; w+=QW,j*=q)
139: continue;
140: c = j - d;
141: cq = c*q;
142: cs = cq<<(L-w);
143: qcs = (((long)(q-1)<<w) + cq) << (L-QW-w);
144: v0 = c - cq;
145: for(i=0,rc=1; i<c; i++,rc*=r) /* rc = r^c */
146: continue;
147: return(w + QW*(rc/(1-z.u) - 1));
148: }
149:
150: struct qlog
151: qlog(x, y, p, u) /*find smallest p so x^p<=y */
152: double x,y,u;
153: long p;
154: {
155: struct qlog z;
156: if(u/x <= y) {
157: z.p = 0;
158: z.u = 1;
159: } else {
160: z = qlog(x, y, p+p, u*u);
161: if(u*z.u/x > y) {
162: z.p += p;
163: z.u *= u;
164: }
165: }
166: return(z);
167: }
168:
169: whuff()
170: {
171: fwrite((char*)&huffcode,sizeof(huffcode),1,stdout);
172: }
173:
174: rhuff(f)
175: FILE *f;
176: {
177: return(read(fileno(f), (char*)&huffcode, sizeof(huffcode)) == sizeof(huffcode));
178: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.