tools/bwt/dcs-bwt/src/bbb_compress.cc

   1 // (C) 2006, Matt Mahoney, mmahoney (at) cs.fit.edu
   2 // Copyright 2007 Google Inc.
   3 //
   4 // This program is free software; you can redistribute it and/or
   5 // modify it under the terms of the GNU General Public License
   6 // as published by the Free Software Foundation; either version 2
   7 // of the License, or (at your option) any later version.
   8 //
   9 // This program is distributed in the hope that it will be useful,
  10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12 // GNU General Public License for more details.
  13 //
  14 // You should have received a copy of the GNU General Public License
  15 // along with this program; if not, write to the Free Software
  16 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  17
  18 #include "bbb_compress.h"
  19 #include "inttypes.h"
  20
  21 #include <cassert>
  22 #include <cstring>
  23 #include <cstdlib>
  24 #include <cstdio>
  25 #include <algorithm>
  26
  27 namespace dcsbwt {
  28 namespace bbb {
  29
  30 namespace {
  31
  32 ///////////////////////// state table ////////////////////////
  33
  34 // State table:
  35 //   nex(state, 0) = next state if bit y is 0, 0 <= state < 256
  36 //   nex(state, 1) = next state if bit y is 1
  37 //   nex(state, 2) = number of zeros in bit history represented by state
  38 //   nex(state, 3) = number of ones represented
  39 //
  40 // States represent a bit history within some context.
  41 // State 0 is the starting state (no bits seen).
  42 // States 1-30 represent all possible sequences of 1-4 bits.
  43 // States 31-252 represent a pair of counts, (n0,n1), the number
  44 //   of 0 and 1 bits respectively.  If n0+n1 < 16 then there are
  45 //   two states for each pair, depending on if a 0 or 1 was the last
  46 //   bit seen.
  47 // If n0 and n1 are too large, then there is no state to represent this
  48 // pair, so another state with about the same ratio of n0/n1 is substituted.
  49 // Also, when a bit is observed and the count of the opposite bit is large,
  50 // then part of this count is discarded to favor newer data over old.
  51
  52 static const uint8 State_table[256][4]={
  53   {  1,  2, 0, 0},{  3,  5, 1, 0},{  4,  6, 0, 1},{  7, 10, 2, 0}, // 0-3
  54   {  8, 12, 1, 1},{  9, 13, 1, 1},{ 11, 14, 0, 2},{ 15, 19, 3, 0}, // 4-7
  55   { 16, 23, 2, 1},{ 17, 24, 2, 1},{ 18, 25, 2, 1},{ 20, 27, 1, 2}, // 8-11
  56   { 21, 28, 1, 2},{ 22, 29, 1, 2},{ 26, 30, 0, 3},{ 31, 33, 4, 0}, // 12-15
  57   { 32, 35, 3, 1},{ 32, 35, 3, 1},{ 32, 35, 3, 1},{ 32, 35, 3, 1}, // 16-19
  58   { 34, 37, 2, 2},{ 34, 37, 2, 2},{ 34, 37, 2, 2},{ 34, 37, 2, 2}, // 20-23
  59   { 34, 37, 2, 2},{ 34, 37, 2, 2},{ 36, 39, 1, 3},{ 36, 39, 1, 3}, // 24-27
  60   { 36, 39, 1, 3},{ 36, 39, 1, 3},{ 38, 40, 0, 4},{ 41, 43, 5, 0}, // 28-31
  61   { 42, 45, 4, 1},{ 42, 45, 4, 1},{ 44, 47, 3, 2},{ 44, 47, 3, 2}, // 32-35
  62   { 46, 49, 2, 3},{ 46, 49, 2, 3},{ 48, 51, 1, 4},{ 48, 51, 1, 4}, // 36-39
  63   { 50, 52, 0, 5},{ 53, 43, 6, 0},{ 54, 57, 5, 1},{ 54, 57, 5, 1}, // 40-43
  64   { 56, 59, 4, 2},{ 56, 59, 4, 2},{ 58, 61, 3, 3},{ 58, 61, 3, 3}, // 44-47
  65   { 60, 63, 2, 4},{ 60, 63, 2, 4},{ 62, 65, 1, 5},{ 62, 65, 1, 5}, // 48-51
  66   { 50, 66, 0, 6},{ 67, 55, 7, 0},{ 68, 57, 6, 1},{ 68, 57, 6, 1}, // 52-55
  67   { 70, 73, 5, 2},{ 70, 73, 5, 2},{ 72, 75, 4, 3},{ 72, 75, 4, 3}, // 56-59
  68   { 74, 77, 3, 4},{ 74, 77, 3, 4},{ 76, 79, 2, 5},{ 76, 79, 2, 5}, // 60-63
  69   { 62, 81, 1, 6},{ 62, 81, 1, 6},{ 64, 82, 0, 7},{ 83, 69, 8, 0}, // 64-67
  70   { 84, 71, 7, 1},{ 84, 71, 7, 1},{ 86, 73, 6, 2},{ 86, 73, 6, 2}, // 68-71
  71   { 44, 59, 5, 3},{ 44, 59, 5, 3},{ 58, 61, 4, 4},{ 58, 61, 4, 4}, // 72-75
  72   { 60, 49, 3, 5},{ 60, 49, 3, 5},{ 76, 89, 2, 6},{ 76, 89, 2, 6}, // 76-79
  73   { 78, 91, 1, 7},{ 78, 91, 1, 7},{ 80, 92, 0, 8},{ 93, 69, 9, 0}, // 80-83
  74   { 94, 87, 8, 1},{ 94, 87, 8, 1},{ 96, 45, 7, 2},{ 96, 45, 7, 2}, // 84-87
  75   { 48, 99, 2, 7},{ 48, 99, 2, 7},{ 88,101, 1, 8},{ 88,101, 1, 8}, // 88-91
  76   { 80,102, 0, 9},{103, 69,10, 0},{104, 87, 9, 1},{104, 87, 9, 1}, // 92-95
  77   {106, 57, 8, 2},{106, 57, 8, 2},{ 62,109, 2, 8},{ 62,109, 2, 8}, // 96-99
  78   { 88,111, 1, 9},{ 88,111, 1, 9},{ 80,112, 0,10},{113, 85,11, 0}, // 100-103
  79   {114, 87,10, 1},{114, 87,10, 1},{116, 57, 9, 2},{116, 57, 9, 2}, // 104-107
  80   { 62,119, 2, 9},{ 62,119, 2, 9},{ 88,121, 1,10},{ 88,121, 1,10}, // 108-111
  81   { 90,122, 0,11},{123, 85,12, 0},{124, 97,11, 1},{124, 97,11, 1}, // 112-115
  82   {126, 57,10, 2},{126, 57,10, 2},{ 62,129, 2,10},{ 62,129, 2,10}, // 116-119
  83   { 98,131, 1,11},{ 98,131, 1,11},{ 90,132, 0,12},{133, 85,13, 0}, // 120-123
  84   {134, 97,12, 1},{134, 97,12, 1},{136, 57,11, 2},{136, 57,11, 2}, // 124-127
  85   { 62,139, 2,11},{ 62,139, 2,11},{ 98,141, 1,12},{ 98,141, 1,12}, // 128-131
  86   { 90,142, 0,13},{143, 95,14, 0},{144, 97,13, 1},{144, 97,13, 1}, // 132-135
  87   { 68, 57,12, 2},{ 68, 57,12, 2},{ 62, 81, 2,12},{ 62, 81, 2,12}, // 136-139
  88   { 98,147, 1,13},{ 98,147, 1,13},{100,148, 0,14},{149, 95,15, 0}, // 140-143
  89   {150,107,14, 1},{150,107,14, 1},{108,151, 1,14},{108,151, 1,14}, // 144-147
  90   {100,152, 0,15},{153, 95,16, 0},{154,107,15, 1},{108,155, 1,15}, // 148-151
  91   {100,156, 0,16},{157, 95,17, 0},{158,107,16, 1},{108,159, 1,16}, // 152-155
  92   {100,160, 0,17},{161,105,18, 0},{162,107,17, 1},{108,163, 1,17}, // 156-159
  93   {110,164, 0,18},{165,105,19, 0},{166,117,18, 1},{118,167, 1,18}, // 160-163
  94   {110,168, 0,19},{169,105,20, 0},{170,117,19, 1},{118,171, 1,19}, // 164-167
  95   {110,172, 0,20},{173,105,21, 0},{174,117,20, 1},{118,175, 1,20}, // 168-171
  96   {110,176, 0,21},{177,105,22, 0},{178,117,21, 1},{118,179, 1,21}, // 172-175
  97   {110,180, 0,22},{181,115,23, 0},{182,117,22, 1},{118,183, 1,22}, // 176-179
  98   {120,184, 0,23},{185,115,24, 0},{186,127,23, 1},{128,187, 1,23}, // 180-183
  99   {120,188, 0,24},{189,115,25, 0},{190,127,24, 1},{128,191, 1,24}, // 184-187
 100   {120,192, 0,25},{193,115,26, 0},{194,127,25, 1},{128,195, 1,25}, // 188-191
 101   {120,196, 0,26},{197,115,27, 0},{198,127,26, 1},{128,199, 1,26}, // 192-195
 102   {120,200, 0,27},{201,115,28, 0},{202,127,27, 1},{128,203, 1,27}, // 196-199
 103   {120,204, 0,28},{205,115,29, 0},{206,127,28, 1},{128,207, 1,28}, // 200-203
 104   {120,208, 0,29},{209,125,30, 0},{210,127,29, 1},{128,211, 1,29}, // 204-207
 105   {130,212, 0,30},{213,125,31, 0},{214,137,30, 1},{138,215, 1,30}, // 208-211
 106   {130,216, 0,31},{217,125,32, 0},{218,137,31, 1},{138,219, 1,31}, // 212-215
 107   {130,220, 0,32},{221,125,33, 0},{222,137,32, 1},{138,223, 1,32}, // 216-219
 108   {130,224, 0,33},{225,125,34, 0},{226,137,33, 1},{138,227, 1,33}, // 220-223
 109   {130,228, 0,34},{229,125,35, 0},{230,137,34, 1},{138,231, 1,34}, // 224-227
 110   {130,232, 0,35},{233,125,36, 0},{234,137,35, 1},{138,235, 1,35}, // 228-231
 111   {130,236, 0,36},{237,125,37, 0},{238,137,36, 1},{138,239, 1,36}, // 232-235
 112   {130,240, 0,37},{241,125,38, 0},{242,137,37, 1},{138,243, 1,37}, // 236-239
 113   {130,244, 0,38},{245,135,39, 0},{246,137,38, 1},{138,247, 1,38}, // 240-243
 114   {140,248, 0,39},{249,135,40, 0},{250, 69,39, 1},{ 80,251, 1,39}, // 244-247
 115   {140,252, 0,40},{249,135,41, 0},{250, 69,40, 1},{ 80,251, 1,40}, // 248-251
 116   {140,252, 0,41}};  // 253-255 are reserved
 117
 118 #define nex(state,sel) State_table[state][sel]
 119
 120 ///////////////////////////// Squash //////////////////////////////
 121
 122 // return p = 1/(1 + exp(-d)), d scaled by 8 bits, p scaled by 12 bits
 123 int squash(int d) {
 124   static const int t[33]={
 125     1,2,3,6,10,16,27,45,73,120,194,310,488,747,1101,
 126     1546,2047,2549,2994,3348,3607,3785,3901,3975,4022,
 127     4050,4068,4079,4085,4089,4092,4093,4094};
 128   if (d>2047) return 4095;
 129   if (d<-2047) return 0;
 130   int w=d&127;
 131   d=(d>>7)+16;
 132   return (t[d]*(128-w)+t[(d+1)]*w+64) >> 7;
 133 }
 134
 135 //////////////////////////// Stretch ///////////////////////////////
 136
 137 // Inverse of squash. d = ln(p/(1-p)), d scaled by 8 bits, p by 12 bits.
 138 // d has range -2047 to 2047 representing -8 to 8.  p has range 0 to 4095.
 139
 140 class Stretch {
 141   Array<short> t;
 142 public:
 143   Stretch();
 144   int operator()(int p) const {
 145     assert(p>=0 && p<4096);
 146     return t[p];
 147   }
 148 } stretch;
 149
 150 Stretch::Stretch(): t(4096) {
 151   int pi=0;
 152   for (int x=-2047; x<=2047; ++x) {  // invert squash()
 153     int i=squash(x);
 154     for (int j=pi; j<=i; ++j)
 155       t[j]=x;
 156     pi=i+1;
 157   }
 158   t[4095]=2047;
 159 }
 160
 161 }  // namespace
 162
 163 //////////////////////////// StateMap //////////////////////////
 164
 165 StateMap::StateMap(): cxt(0), t(256) {
 166   for (int i=0; i<256; ++i) {
 167     int n0=nex(i,2);
 168     int n1=nex(i,3);
 169     if (n0==0) n1*=128;
 170     if (n1==0) n0*=128;
 171     t[i] = 65536*(n1+1)/(n0+n1+2);
 172   }
 173 }
 174
 175 //////////////////////////// APM //////////////////////////////
 176
 177 // maps p, cxt -> p initially
 178 APM::APM(int n): index(0), N(n), t(n*33) {
 179   for (int i=0; i<N; ++i)
 180     for (int j=0; j<33; ++j)
 181       t[i*33+j] = i==0 ? squash((j-16)*128)*16 : t[j];
 182 }
 183
 184 int APM::p(int y, int pr, int cxt, int rate) {
 185   assert(pr>=0 && pr<4096 && cxt>=0 && cxt<N && rate>0 && rate<32);
 186   pr = stretch(pr);
 187   int g = (y<<16)+(y<<rate)-y-y;
 188   t[index] += (g-t[index]) >> rate;
 189   t[index+1] += (g-t[index+1]) >> rate;
 190   const int w = pr&127;  // interpolation weight (33 points)
 191   index = ((pr+2048)>>7)+cxt*33;
 192   return (t[index]*(128-w)+t[index+1]*w) >> 11;
 193 }
 194
 195 //////////////////////////// Predictor //////////////////////////
 196
 197 void Predictor::update(int y) {
 198   // update model
 199   *cp=nex(*cp, y);
 200
 201   // update context
 202   c0+=c0+y;
 203   if (++bpos==8) {
 204     bpos=0;
 205     c4=(c4<<8)|(c0-256);
 206     c0=1;
 207     bpos=0;
 208     if (((c4^c4>>8)&255)==0) {
 209       if (run<65535)
 210         ++run;
 211       if (run==1 || run==2 || run==4) runcxt += 1;
 212     }
 213     else run=0, runcxt=0;
 214   }
 215
 216   int c1c = ((c4 & 255) + 256) >> (8-bpos);
 217   int c1d = (c4 >> (7-bpos)) & 1;
 218   if (c1c == c0) c1d += 2;
 219   //int c1 = c4 & (255 ^ (127>>bpos));
 220
 221   // predict
 222   cp = &t1[c0];
 223   pr = sm.p(y, *cp);
 224   pr = a2.p(y, pr, c0|(c1d<<8), 7);
 225   pr = a3.p(y, pr, runcxt, 8);
 226   pr = (a5.p(y, pr, c0^((c4*123456791)>>21), 7)+pr+1)>>1;
 227   pr = (a4.p(y, pr, c0, 7)+pr+1)>>1;
 228 }
 229
 230 }  // namespace bbb
 231
 232 ///////////////////// BbbCompressor /////////////////////
 233
 234 void BbbCompressor::CompressBit(int bit) {
 235   int p = predictor_.p();
 236   assert(p >= 0);
 237   assert(p < 4096);
 238   p += (p < 2048);
 239   uint32 xmid = x1_ + ((x2_-x1_) >> 12) * p
 240       + ((((x2_-x1_) & 0xfff) * p) >> 12);
 241   assert(xmid >= x1_);
 242   assert(xmid < x2_);
 243   bit ? (x2_ = xmid) : (x1_ = xmid + 1);
 244   predictor_.update(bit);
 245   while (((x1_ ^ x2_) & 0xff000000) == 0) {
 246     EmitByte(x2_>>24);
 247     x1_ <<= 8;
 248     x2_ = (x2_<<8) + 255;
 249   }
 250 }
 251
 252 void BbbCompressor::WriteBegin() {
 253   x1_ = 0;
 254   x2_ = 0xffffffff;
 255 }
 256
 257 void BbbCompressor::WriteEndPrivate() {
 258   EmitByte(x1_ >> 24);
 259   EmitByte(255);
 260   EmitByte(255);
 261   EmitByte(255);
 262 }
 263
 264 ///////////////////// BbbDecompressor /////////////////////
 265
 266 void BbbDecompressor::ReadBegin() {
 267   x1_ = 0;
 268   x2_ = 0xffffffff;
 269   x_ = GetCompressedByte();
 270   for (int i = 3; i; --i) x_ = (x_ << 8) + GetCompressedByte();
 271 }
 272
 273 int BbbDecompressor::GetUncompressedBit() {
 274   int p = predictor_.p();
 275   assert(p >= 0);
 276   assert(p < 4096);
 277   p += (p < 2048);
 278   uint32 xmid = x1_ + ((x2_-x1_) >> 12) * p
 279       + ((((x2_-x1_) & 0xfff) * p) >> 12);
 280   assert(xmid >= x1_);
 281   assert(xmid < x2_);
 282   int bit = (x_ <= xmid);
 283   bit ? (x2_ = xmid) : (x1_ = xmid + 1);
 284   predictor_.update(bit);
 285   while (((x1_ ^ x2_) & 0xff000000) == 0) {
 286     x1_ <<= 8;
 287     x2_ = (x2_<<8) + 255;
 288     x_ =(x_<<8) + GetCompressedByte();
 289   }
 290   return bit;
 291 }
 292
 293 }  // namespace dcsbwt