Annotation of researchv8dc/cmd/spell/huff.c, revision 1.1.1.1

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: }

unix.superglobalmegacorp.com

This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.