OK, here is the encoder and decoder from libzpaq. It is public domain (actually MIT modified to remove the restriction of including the license, so same thing). It requires a bit of explanation. ZPAQ is an archiving format that mixes compressed and uncompressed data, so it needs a way to identify where the compressed data ends. It does this with a marker of 4 zero bytes. To support this, the encoder ensures that 4 zero bytes never appears in the compressed data.

The method Encoder::compress(c) compresses a byte c in the range 0..255, or -1 at EOF. Decoder::decompress() returns a byte 0..255 or -1 at EOF. compress() encodes 9 bits either a 0 followed by 8 bits of c in MSB to LSB order, or a 1 to indicate EOF. Each bit of c is coded with a prediction, p, for that bit in the range 0..65535 estimating the probability that the bit is a 1 as a 16 bit number. For the EOF bit, p = 0. For the bits of c, the predictions come from a Predictor, which has the following 3 methods:

p() returns next bit prediction in the range 0..32767 as a 15 bit number. The encoder scales this up to an odd 16 bit number, 2*p+1.

update(bit) is passed 0 or 1 for the true value of the last predicted bit.

init() initializes. A predictor is initialized with a ZPAQL object which contains code describing how to model the data.

The encoder has initial state low = 1, high = 0xFFFFFFFF. To encode a bit with probability p (that it is a 1), the range is split into two parts, low..mid and mid+1..high, with a 1 encoded in the lower half with size proportional to p. The statement:

Code:

U32 mid=low+((high-low)>>16)*p+((((high-low)&0xffff)*p)>>16);

multiplies in two 16 bit parts to avoid 32 bit overflow. It could be simplified using 64 bit arithmetic to

Code:

uint32_t mid=low+((uint64_t(high-low)*p)>>16);

U32 is uint32_t. The statement

ensures that 4 zero bytes are never output. If you don't care, you could remove it from both the encoder and decoder.

in and out are abstract data types Reader and Writer. A Reader has virtual method get() which returns a byte or -1 at EOF. Writer has put(c) which outputs a byte.

As the range (high - low) gets smaller, the leading bytes of 32 bit values high and low are flushed so that they are always different. During compression, these bytes are written. During decompression, these bytes are discarded and a byte is read into the low end of curr. You can imagine that low and high are infinite precision numbers with implicit trailing 0 bits after low and 1 bits after high.

The EOF bit preceding each byte is always modeled with p = 0, not using the predictor. This has the effect of setting high = low and forcing 4 bytes to be flushed and the new state high = 0xFFFFFFFF and low = one of 1, 0x100, 0x10000, or 0x1000000. During encoding, the application would then write 4 zero bytes. During decoding, these bytes would be read into curr = 0, resulting in curr < low, which would be an error (corrupted input) if it occurred anywhere else. The decoder uses curr = 0 as a signal to initialize the decoder by reading the first 4 bytes of input into curr.

Code:

// Encoder compresses using an arithmetic code
class Encoder {
public:
Encoder(ZPAQL& z):
out(0), low(1), high(0xFFFFFFFF), pr(z) {}
void init();
void compress(int c); // c is 0..255 or EOF
Writer* out; // destination
private:
U32 low, high; // range
Predictor pr; // to get p
void encode(int y, int p); // encode bit y (0..1) with probability p (0..65535)
};
// Initialize for start of block
void Encoder::init() {
low=1;
high=0xFFFFFFFF;
pr.init();
}
// compress bit y having probability p/64K
void Encoder::encode(int y, int p) {
assert(out);
assert(p>=0 && p<65536);
assert(y==0 || y==1);
assert(high>low && low>0);
U32 mid=low+((high-low)>>16)*p+((((high-low)&0xffff)*p)>>16); // split range
assert(high>mid && mid>=low);
if (y) high=mid; else low=mid+1; // pick half
while ((high^low)<0x1000000) { // write identical leading bytes
out->put(high>>24); // same as low>>24
high=high<<8|255;
low=low<<8;
low+=(low==0); // so we don't code 4 0 bytes in a row
}
}
// compress byte c (0..255 or -1=EOS)
void Encoder::compress(int c) {
assert(out);
if (c==-1)
encode(1, 0);
else {
assert(c>=0 && c<=255);
encode(0, 0);
for (int i=7; i>=0; --i) {
int p=pr.predict()*2+1;
assert(p>0 && p<65536);
int y=c>>i&1;
encode(y, p);
pr.update(y);
}
}
}
// Decoder decompresses using an arithmetic code
class Decoder {
public:
Reader* in; // destination
Decoder(ZPAQL& z);
int decompress(); // return a byte or EOF
void init() {pr.init(); low=1; high=0xFFFFFFFF; curr=0;}
private:
U32 low, high; // range
U32 curr; // last 4 bytes of archive
Predictor pr; // to get p
int decode(int p); // return decoded bit (0..1) with prob. p (0..65535)
};
Decoder::Decoder(ZPAQL& z):
in(0), low(1), high(0xFFFFFFFF), curr(0), pr(z) {}
// Return next bit of decoded input, which has 16 bit probability p of being 1
int Decoder::decode(int p) {
assert(p>=0 && p<65536);
assert(high>low && low>0);
if (curr<low || curr>high) error("archive corrupted");
assert(curr>=low && curr<=high);
U32 mid=low+((high-low)>>16)*p+((((high-low)&0xffff)*p)>>16); // split range
assert(high>mid && mid>=low);
int y=curr<=mid;
if (y) high=mid; else low=mid+1; // pick half
while ((high^low)<0x1000000) { // shift out identical leading bytes
high=high<<8|255;
low=low<<8;
low+=(low==0);
int c=in->get();
if (c<0) error("unexpected end of file");
curr=curr<<8|c;
}
return y;
}
// Decompress 1 byte or -1 at end of input
int Decoder::decompress() {
if (curr==0) { // segment initialization
for (int i=0; i<4; ++i)
curr=curr<<8|in->get();
}
if (decode(0)) {
if (curr!=0) error("decoding end of stream");
return -1;
}
else {
int c=1;
while (c<256) { // get 8 bits
int p=pr.predict()*2+1;
c+=c+decode(p);
pr.update(c&1);
}
return c-256;
}
}