=p && q0 && n<=N);
assert(c_>=0 && c_0);
while (sum>2000000000/PSCALE) sum/=4, n1/=4;
return int((PSCALE-1)*n1/sum);
}
// Adjust the weights by gradient descent to reduce cost of bit y
void Mixer::update(int y) {
U32 s0=0, s1=0;
for (int i=0; i0 && s1>0) {
const U32 s=s0+s1;
const U32 sy=y?s1:s0;
const U32 sy1=0xffffffff/sy+(rnd()&1023) >> 10;
const U32 s1 =0xffffffff/s +(rnd()&1023) >> 10;
for (int i=0; i> 8;
wt[c][i]=min(65535, max(1, int(wt[c][i]+dw)));
}
}
n=0;
}
Mixer::Mixer(int C_): C(C_), bc0(new U32[N]), bc1(new U32[N]),
wt(new U32[C_][N]), n(0), c(0) {
for (int i=0; i> 6)+4*(ch(2) >> 6));
return (p1+p2)/2;
}
void update(int y) {
m1.update(y);
m2.update(y);
}
U32 getC() const {return 256;}
U32 getN() const {return m1.getN();}
};
MultiMixer mixer;
//////////////////////////// CounterMap ////////////////////////////
/* CounterMap maintains a model and one context
Countermap cm(N); -- Create, size 2^N bytes
cm.update(h); -- Update model, then set next context hash to h
cm.write(); -- Predict next bit and write counts to mixer
cm.add(); -- Predict and add to previously written counts
There should be 8 calls to either write() or add() between each update(h).
h is a 32-bit hash of the context which should be set after a whole number
of bytes are read. */
// Stores only the most recent byte and its count per context
// in a hash table without collision detection
class CounterMap1 {
const int N;
struct S {
U8 c; // char
U8 n; // count
S(): c(0), n(0) {}
};
S* t; // cxt -> c repeated last n times
U32 cxt;
public:
CounterMap1(int n): N(n-1), t(new S[1<> 32-N;
}
void add() {
if ((U32)((t[cxt].c+256) >> 8-ch.bpos())==ch()) {
if ((t[cxt].c >> 7-ch.bpos()) & 1)
mixer.add(0, t[cxt].n);
else
mixer.add(t[cxt].n, 0);
}
}
void write() {
mixer.write(0, 0);
add();
}
};
// Uses a nibble-oriented hash table of contexts
class CounterMap2 {
private:
const U32 N2; // Size of ht2 in elements
U32 cxt; // Major context
Hashtable ht2; // Secondary hash table
Counter* cp[8]; // Pointers into ht2 or 0 if not used
public:
CounterMap2(int n); // Use 2^n bytes memory
void add();
void update(U32 h);
void write() {
mixer.write(0, 0);
add();
}
};
CounterMap2::CounterMap2(int n): N2(n), cxt(0), ht2(N2) {
for (int i=0; i<8; ++i)
cp[i]=0;
}
// Predict the next bit given the bits so far in ch()
void CounterMap2::add() {
const U32 bcount = ch.bpos();
if (bcount==4) {
cxt^=hash(ch.hi(), cxt);
ht2.set(cxt);
}
cp[bcount]=&ht2(ch.lo());
mixer.add(cp[bcount]->get0(), cp[bcount]->get1());
}
// After 8 predictions, update the models with the last input char, ch(1),
// then set the new context hash to h
void CounterMap2::update(U32 h) {
const U32 c=ch(1);
// Update the secondary context
for (int i=0; i<8; ++i) {
if (cp[i]) {
cp[i]->add((c>>(7-i))&1);
cp[i]=0;
}
}
cxt=h;
ht2.set(cxt);
}
// Combines 1 and 2 above.
class CounterMap3 {
CounterMap1 cm1;
CounterMap2 cm2;
public:
CounterMap3(int n): cm1(n-2), cm2(n) {}
void update(U32 h) {
cm1.update(h);
cm2.update(h);
}
void write() {
cm2.write();
cm1.add();
}
void add() {
cm2.add();
cm1.add();
}
};
#define CounterMap CounterMap3
//////////////////////////// Model ////////////////////////////
// All models have a function model() which updates the model with the
// last bit of input (in ch) then writes probabilities for the following
// bit into mixer.
class Model {
public:
virtual void model() = 0;
virtual ~Model() {}
};
//////////////////////////// defaultModel ////////////////////////////
// DefaultModel predicts P(1) = 0.5
class DefaultModel: public Model {
public:
void model() {mixer.write(1, 1);}
} defaultModel;
//////////////////////////// charModel ////////////////////////////
// A CharModel contains n-gram models from 0 to 7
class CharModel: public Model {
enum {N=8}; // Number of models
Counter *t0, *t1; // Model orders 0, 1 [256], [65536]
CounterMap t2, t3, t4, t5, t6, t7;
U32 *cxt; // Context hashes [N]
Counter *cp0, *cp1; // Pointers to counters in t0, t1
public:
CharModel(): t0(new Counter[256]), t1(new Counter[65536]),
t2(12+MEM), t3(14+MEM), t4(16+MEM), t5(16+MEM),
t6(16+MEM), t7(16+MEM), cxt(new U32[N]) {
cp0=&t0[0];
cp1=&t1[0];
memset(cxt, 0, N*sizeof(U32));
memset(t0, 0, 256*sizeof(Counter));
memset(t1, 0, 65536*sizeof(Counter));
}
void model(); // Update and predict
} charModel;
// Update with bit y, put array of 0 counts in n0 and 1 counts in n1
inline void CharModel::model() {
// Update models
int y = ch(ch.bpos()==0)&1; // last input bit
cp0->add(y);
cp1->add(y);
// Update context
if (ch.bpos()==0) { // Start new byte
for (int i=N-1; i>0; --i)
cxt[i]=cxt[i-1]^hash(ch(1), i);
t2.update(cxt[2]);
t3.update(cxt[3]);
t4.update(cxt[4]);
t5.update(cxt[5]);
t6.update(cxt[6]);
t7.update(cxt[7]);
}
cp0=&t0[ch()];
cp1=&t1[ch()+256*ch(1)];
// Write predictions to the mixer
mixer.write(cp0->get0(), cp0->get1());
mixer.write(cp1->get0(), cp1->get1());
t2.write();
t3.add();
t4.write();
t5.add();
t6.write();
t7.add();
}
//////////////////////////// matchModel ////////////////////////////
/* A MatchModel looks for a match of length n >= 8 bytes between
the current context and previous input, and predicts the next bit
in the previous context with weight n. If the next bit is 1, then
the mixer is assigned (0, n), else (n, 0). Matchies are found using
an index (a hash table of pointers into ch). */
class MatchModel: public Model {
enum {M=4}; // Number of strings to match
U32 hash[2]; // Hashes of current context up to pos-1
U32 begin[M]; // Points to first matching byte
U32 end[M]; // Points to last matching byte + 1, 0 if no match
U32 *ptr; // Hash table of pointers [2^(MEM+14)]
public:
MatchModel(): ptr(new U32[1 << (MEM+14)]) {
memset(ptr, 0, (1 << (MEM+14))*sizeof(U32));
hash[0]=hash[1]=0;
for (int i=0; i> (18-MEM);
if ((hash[0]>>28)==0)
h=hash[1] >> (18-MEM); // 1/16 of 8-contexts are hashed to 32 bytes
for (int i=0; i0) {
begin[i]=end[i];
U32 p=ch.pos();
while (begin[i]>0 && p>0 && begin[i]!=p+1
&& ch[begin[i]-1]==ch[p-1]) {
--begin[i];
--p;
}
}
if (end[i]==begin[i]) // No match found
begin[i]=end[i]=0;
break;
}
}
ptr[h]=ch.pos();
}
// Test whether the current context is valid in the last 0-7 bits
for (int i=0; i> (8-ch.bpos())) != ch())
begin[i]=end[i]=0;
}
// Predict the bit found in the matching contexts
int n0=0, n1=0;
for (int i=0; i511)
wt=511;
int y=(ch[end[i]]>>(7-ch.bpos()))&1;
if (y)
n1+=wt;
else
n0+=wt;
}
}
mixer.write(n0, n1);
}
//////////////////////////// recordModel ////////////////////////////
/* A RecordModel finds fixed length records and models bits in the context
of the two bytes above (the same position in the two previous records)
and in the context of the byte above and to the left (the previous byte).
The record length is assumed to be the interval in the most recent
occurrence of a byte occuring 4 times in a row equally spaced, e.g.
"x..x..x..x" would imply a record size of 3 (the mimimum). */
class RecordModel: public Model {
enum {N=2}; // Number of models
CounterMap t0, t1; // Model
int repeat; // Cycle length
public:
RecordModel(): t0(14+MEM), t1(14+MEM), repeat(1) {}
void model();
} recordModel;
// Update the model with bit y, then put predictions of the next update
// as 0 counts in n0[0..N-1] and 1 counts in n1[0..N-1]
inline void RecordModel::model() {
if (ch.bpos()==0) {
// Check for a repeating pattern of interval 3 or more
const int c=ch(1);
const int d1=ch.pos(c,0)-ch.pos(c,1);
const int d2=ch.pos(c,1)-ch.pos(c,2);
const int d3=ch.pos(c,2)-ch.pos(c,3);
if (d1>2 && d1==d2 && d2==d3)
repeat=d1;
// Compute context hashes
t0.update(hash(ch(repeat), ch(repeat*2), repeat&255)); // 2 above
t1.update(hash(ch(1), ch(repeat), repeat&255)); // above and left
}
t0.write();
t1.add();
}
//////////////////////////// sparseModel ////////////////////////////
// A SparseModel models several order-2 contexts with gaps
class SparseModel: public Model {
enum {N=6}; // Number of models
CounterMap t0, t1, t2, t3, t4, t5; // Sparse models
public:
SparseModel(): t0(12+MEM), t1(12+MEM), t2(12+MEM),
t3(12+MEM), t4(12+MEM), t5(12+MEM) {}
void model(); // Update and predict
} sparseModel;
// Update with bit y, put array of 0 counts in n0 and 1 counts in n1
inline void SparseModel::model() {
// Update context
if (ch.bpos()==0) { // Start new byte
t0.update(hash(ch(1), ch(3)));
t1.update(hash(ch(1), ch(4)));
t2.update(hash(ch(2), ch(4)));
t3.update(hash(ch(4), ch(8)));
t4.update(hash(ch(2), ch(3)));
t5.update(hash(ch(3), ch(4)));
}
// Predict
t0.write();
t1.add();
t2.write();
t3.add();
t4.write();
t5.add();
}
//////////////////////////// analogModel ////////////////////////////
// An AnalogModel is intended for 16-bit mono or stereo (WAV files)
// 24-bit images (BMP files), and 8 bit analog data (such as grayscale
// images).
class AnalogModel: public Model {
enum {N=6};
CounterMap t0, t1, t2, t3, t4, t5;
int pos3; // pos % 3
public:
AnalogModel(): t0(10+MEM), t1(10+MEM), t2(10+MEM), t3(10+MEM),
t4(10+MEM), t5 (10+MEM), pos3(0) {}
void model() {
if (ch.bpos()==0) {
if (++pos3==3) pos3=0;
t0.update(hash(ch(2)/4, ch(4)/4, ch.pos()%2)); // 16 bit mono model
t1.update(hash(ch(2)/16, ch(4)/16, ch.pos()%2));
t2.update(hash(ch(2)/4, ch(4)/4, ch(8)/4, ch.pos()%4)); // Stereo
t3.update(hash(ch(3), ch(6)/4, pos3)); // 24 bit image models
t4.update(hash(ch(1)/16, ch(2)/16, ch(3)/4, pos3));
t5.update(hash(ch(1)/2, ch(2)/8, ch(3)/32)); // 8-bit data model
}
t0.write();
t1.add();
t2.add();
t3.write();
t4.add();
t5.write();
}
} analogModel;
//////////////////////////// wordModel ////////////////////////////
// A WordModel models words, which are any characters > 32 separated
// by whitespace ( <= 32). There is a unigram, bigram and sparse
// bigram model (skipping 1 word).
class WordModel: public Model {
enum {N=3};
CounterMap t0, t1, t2;
U32 cxt[N]; // Hashes of last N words
public:
WordModel(): t0(16+MEM), t1(15+MEM), t2(15+MEM) {
for (int i=0; i32) {
cxt[0]^=hash(cxt[0], c);
}
else if (cxt[0]) {
for (int i=N-1; i>0; --i)
cxt[i]=cxt[i-1];
cxt[0]=0;
}
t0.update(cxt[0]);
t1.update(cxt[1]+cxt[0]);
t2.update(cxt[2]+cxt[0]);
}
t0.write();
t1.write();
t2.add();
}
} wordModel;
//////////////////////////// Predictor ////////////////////////////
/* A Predictor adjusts the model probability using SSE and passes it
to the encoder. An SSE model is a table of counters, sse[SSE1][SSE2]
which maps a context and a probability into a new, more accurate
probability. The context, SSE1, consists of the 0-7 bits of the current
byte and the 2 leading bits of the previous byte. The probability
to be mapped, SSE2 is first stretched near 0 and 1 using SSEMap, then
quantized into SSE2=32 intervals. Each SSE element is a pair of 0
and 1 counters of the bits seen so far in the current context and
probability range. Both the bin below and above the current probability
is updated by adding 1 to the appropriate count (n0 or n1). The
output probability for an SSE element is n1/(n0+n1) interpolated between
the bins below and above the input probability. This is averaged
with the original probability with 25% weight to give the final
probability passed to the encoder. */
class Predictor {
enum {SSE1=256*4, SSE2=32, // SSE dimensions (contexts, probability bins)
SSESCALE=1024/SSE2}; // Number of mapped probabilities between bins
// Scale probability p into a context in the range 0 to 1K-1 by
// stretching the ends of the range.
class SSEMap {
U16 table[PSCALE];
public:
int operator()(int p) const {return table[p];}
SSEMap();
} ssemap; // functoid
// Secondary source encoder element
struct SSEContext {
U8 c1, n; // Count of 1's, count of bits
int p() const {return PSCALE*(c1*64+1)/(n*64+2);}
void update(int y) {
if (y)
++c1;
if (++n>254) { // Roll over count overflows
c1/=2;
n/=2;
}
}
SSEContext(): c1(0), n(0) {}
};
SSEContext (*sse)[SSE2+1]; // [SSE1][SSE2+1] context, mapped probability
U32 nextp; // p()
U32 ssep; // Output of sse
U32 context; // SSE context
public:
Predictor();
int p() const {return nextp;} // Returns pr(y = 1) * PSCALE
void update(int y); // Update model with bit y = 0 or 1
};
Predictor::SSEMap::SSEMap() {
for (int i=0; i1023) p=1023;
if (p<0) p=0;
table[i]=p;
}
}
Predictor::Predictor():
sse((SSEContext(*)[SSE2+1]) new SSEContext[SSE1][SSE2+1]),
nextp(PSCALE/2), ssep(512), context(0) {
// Initialize to sse[context][ssemap(p)] = p
int N=PSCALE;
int oldp=SSE2+1;
for (int i=N-1; i>=0; --i) {
int p=(ssemap(i*PSCALE/N)+SSESCALE/2)/SSESCALE;
int n=1+N*N/((i+1)*(N-i));
if (n>254) n=254;
int c1=(i*n+N/2)/N;
for (int j=oldp-1; j>=p; --j) {
for (int k=0; k=0x4000000) xmid+=(xdiff>>12)*p;
else if (xdiff>=0x100000) xmid+=((xdiff>>6)*p)>>6;
else xmid+=(xdiff*p)>>12;
// Update the range
if (y)
x2=xmid;
else
x1=xmid+1;
predictor.update(y);
// Shift equal MSB's out
while (((x1^x2)&0xff000000)==0) {
putc(x2>>24, archive);
x1<<=8;
x2=(x2<<8)+255;
}
}
/* Decode one bit from the archive, splitting [x1, x2] as in the encoder
and returning 1 or 0 depending on which subrange the archive point x is in.
*/
inline int Encoder::decode() {
// Split the range
const U32 p=predictor.p()*(4096/PSCALE)+2048/PSCALE; // P(1) * 4K
assert(p<4096);
const U32 xdiff=x2-x1;
U32 xmid=x1; // = x1+p*(x2-x1) multiply without overflow, round down
if (xdiff>=0x4000000) xmid+=(xdiff>>12)*p;
else if (xdiff>=0x100000) xmid+=((xdiff>>6)*p)>>6;
else xmid+=(xdiff*p)>>12;
// Update the range
int y=0;
if (x<=xmid) {
y=1;
x2=xmid;
}
else
x1=xmid+1;
predictor.update(y);
// Shift equal MSB's out
while (((x1^x2)&0xff000000)==0) {
x1<<=8;
x2=(x2<<8)+255;
int c=getc(archive);
if (c==EOF) c=0;
x=(x<<8)+c;
}
return y;
}
// Should be called when there is no more to compress
void Encoder::flush() {
// In COMPRESS mode, write out the remaining bytes of x, x1 < x < x2
if (mode==COMPRESS) {
while (((x1^x2)&0xff000000)==0) {
putc(x2>>24, archive);
x1<<=8;
x2=(x2<<8)+255;
}
putc(x2>>24, archive); // First unequal byte
}
}
//////////////////////////// main ////////////////////////////
// Read one byte from encoder e
int decompress(Encoder& e) { // Decompress 8 bits, MSB first
int c=0;
for (int i=0; i<8; ++i)
c=c+c+e.decode();
return c;
}
// Write one byte c to encoder e
void compress(Encoder& e, int c) {
for (int i=7; i>=0; --i)
e.encode((c>>i)&1);
}
// Read and return a line of input from FILE f (default stdin) up to
// first control character except tab. Skips CR in CR LF.
string getline(FILE* f=stdin) {
int c;
string result="";
while ((c=getc(f))!=EOF && (c>=32 || c=='\t'))
result+=char(c);
if (c=='\r')
(void) getc(f);
return result;
}
// User interface
int main(int argc, char** argv) {
// Check arguments
if (argc<2) {
printf(
PROGNAME " file compressor/archiver, (C) 2003, Matt Mahoney, mmahoney@cs.fit.edu\n"
"This program is free software distributed without warranty under the terms\n"
"of the GNU General Public License, see http://www.gnu.org/licenses/gpl.txt\n"
"\n"
"To compress: " PROGNAME " archive filenames... (archive will be created)\n"
" or (MSDOS): dir/b | " PROGNAME " archive (reads file names from input)\n"
"To extract/compare: " PROGNAME " archive (does not clobber existing files)\n"
"To view contents: more < archive\n");
return 1;
}
// File names and sizes from input or archive
vector filename; // List of names
vector filesize; // Size or -1 if error
int uncompressed_bytes=0, compressed_bytes=0; // Input, output sizes
// Extract files
FILE* archive=fopen(argv[1], "rb");
if (archive) {
if (argc>2) {
printf("File %s already exists\n", argv[1]);
return 1;
}
printf("Extracting archive %s ...\n", argv[1]);
// Read PROGNAME "\r\n" at start of archive
if (getline(archive) != PROGNAME) {
printf("Archive file %s not in " PROGNAME " format\n", argv[1]);
return 1;
}
// Read "size filename" in "%d\t%s\r\n" format
while (true) {
string s=getline(archive);
if (s.size()>1) {
filesize.push_back(atol(s.c_str()));
string::iterator tab=find(s.begin(), s.end(), '\t');
if (tab!=s.end())
filename.push_back(string(tab+1, s.end()));
else
filename.push_back("");
}
else
break;
}
// Test end of header for "\f\0"
{
int c1=0, c2=0;
if ((c1=getc(archive))!='\f' || (c2=getc(archive))!=0) {
printf("%s: Bad " PROGNAME " header format %d %d\n", argv[1],
c1, c2);
return 1;
}
}
// Extract files from archive data
Encoder e(DECOMPRESS, archive);
for (int i=0; i2)
for (int i=2; i=0)
fprintf(archive, "%ld\t%s\r\n", filesize[i], filename[i].c_str());
}
putc(032, archive); // MSDOS EOF
putc('\f', archive);
putc(0, archive);
// Write data
Encoder e(COMPRESS, archive);
long file_start=ftell(archive);
for (int i=0; i=0) {
uncompressed_bytes+=size;
printf("%-23s %10ld -> ", filename[i].c_str(), size);
FILE* f=fopen(filename[i].c_str(), "rb");
int c;
for (long j=0; j0 && elapsed_time>0) {
printf(" (%1.4f bpc, %1.2f%% at %1.0f KB/s)",
compressed_bytes*8.0/uncompressed_bytes,
compressed_bytes*100.0/uncompressed_bytes,
uncompressed_bytes/(elapsed_time*1000.0));
}
printf("\n");
return 0;
}