1 /* b_i_t.c                                                                   */
2 /* this implements bit indexed trees for culmulative frequency storeage      */
3 /* described in                                                              */
4 /* Peter Fenwick: A New Data Structure for Cumulative Probability Tables     */
5 /* Technical Report 88, Dep. of Computer Science, University of Auckland, NZ */
6 /* I used it in Pascal in 1993? for random number generation                 */
7 
8 #include "sz_bit.h"
9 #include <stdio.h>
10 #include <stdlib.h>
11 
12 /* returns the culmulative frequency <= the given symbol */
getcf(symb s,cumtbl * tbl)13 Inline freq getcf(symb s, cumtbl *tbl)
14 { freq sum, *cf;
15   sum = (cf=tbl->cf)[s+1] - tbl->f[s];
16   s++;
17   while (s &= s-1)
18     sum += cf[s];
19   return sum;
20 }
21 
22 /* updates the given culmulative frequency */
updatecumonly(symb s,cumtbl * tbl,int delta)23 Inline void updatecumonly(symb s, cumtbl *tbl, int delta)
24 { freq *cf;
25   symb size;
26   tbl->totfreq += delta;
27   s++;
28   cf = tbl->cf;
29   size = tbl->size;
30   while (s<=size) {
31      cf[s] += delta;
32      s = (s | (s-1)) + 1;
33   }
34 }
35 
36 /* get symbol for this culmulative frequency */
getsym(freq f,cumtbl * tbl)37 Inline symb getsym(freq f, cumtbl *tbl)
38 { symb s, m, x, size;
39   freq *cf;
40   m = tbl->mask;
41   size = tbl->size;
42   cf = tbl->cf;
43   s = 0;
44   do {
45     if ((x=s|m) <= size && f >= cf[x]) {
46       f -= cf[x];
47       s = x;
48     }
49   } while (m >>= 1);
50   return s;
51 }
52 
53 #define RECALC(cond) {                   \
54   symb i;                                \
55   freq *cf;                              \
56   cf = tbl->cf;                          \
57   tbl->totfreq = 0;                      \
58   for (i=1; i<=tbl->size; i<<=1) {       \
59     symb j;                              \
60     for (j=i; j<=tbl->size; j+= i<<1) {  \
61       symb k;                            \
62       if (cond)                          \
63         cf[j] = 0;                       \
64       else                               \
65         tbl->totfreq += cf[j] = tbl->f[j-1];\
66       for (k=i>>1; k; k>>=1)             \
67         cf[j] += cf[j-k];                \
68     } /* end for j */                    \
69   } /* end for i */                      \
70 }
71 
72 /* scales the culmulative frequency tables by 0.5 and keeps nonzero values */
scalefreq(cumtbl * tbl)73 void scalefreq(cumtbl *tbl)
74 { freq *f, *endf;
75   for (f=tbl->f, endf = f+tbl->size; f<endf; f++)
76      *f = (*f + 1)>>1;
77   RECALC(0);
78 }
79 
80 /* scales the culmulative frequency tables by 0.5 and keeps nonzero values */
scalefreqcond(cumtbl * tbl,uint * donotuse)81 void scalefreqcond(cumtbl *tbl, uint *donotuse)
82 { freq *f, *endf;
83   for (f=tbl->f, endf = f+tbl->size; f<endf; f++)
84      *f = (*f + 1)>>1;
85   RECALC(donotuse[j-1]!=3);
86 }
87 
88 /* allocates memory for the frequency table and initializes it */
initfreq(cumtbl * tbl,symb tblsize,freq initvalue)89 int initfreq(cumtbl *tbl, symb tblsize, freq initvalue)
90 { symb i;
91   if ((tbl->f = (freq*)malloc(2*tblsize*sizeof(freq))) == NULL) return 1;
92   tbl->cf = tbl->f + tblsize - 1;
93   tbl->size = tblsize;
94   tbl->mask = 1;
95   while (tblsize>>=1)
96      tbl->mask <<=1;
97   for (i=0; i<tbl->size; i++)
98      tbl->f[i] = initvalue;
99   RECALC(0);
100   return 0;
101 }
102 
103 /* does the obvious thing */
freefreq(cumtbl * tbl)104 void freefreq(cumtbl *tbl)
105 { free(tbl->f);
106 }
107