This is the start of the documentation on the maths behind using the
Laplacian distribution for Huffman codes.  Obviously it is isn't finished
yet.

Define the single sided Laplacian density as:

p(x) = { \sqrt 2 \over \sigma } e^{{- \sqrt 2 \over \sigma} x}

Then:

\int_0^\infinity p(x)  dx = 1

\int_0^\infinity x p(x) dx = {\sigma \over \sqrt 2}

\int_0^\infinity x^2 p(x) dx = \sigma^2

With E(|x|) being the expectation of the absolute value of x, and being
calculated as sumx / blocksize.

We need resn such that half the samples lie below 2^resn.  This ensures that
the Huffman code is resn + 1 bits long with probability 0.5, resn + 2 bits
long with probability 0.25, resn + n bits long with probability 2^{-n},
which is optimal.

Hence with y = 2^resn

0.5 = \int_0^y p(x) dx = - e^{{- \sqrt 2 \over \sigma} y} - 1

y = {ln(2) \over \sqrt 2} \sigma

  = ln(2) E(|x|)

Putting it all together:

resn = log2(ln(2) * sumx / blocksize)

and the code does this, with a bit of fiddling for rounding errors and
implementation of log2().

Similarly, the lpc routine returns the summed squared prediction residual,
e = E(x^2) * blocksize

E(x^2) = \int_0^\infinity x^2 p(x) dx

       = \sigma^2

y = {ln(2) \over \sqrt 2} \sigma

  = \sqrt { ln(2)^2 \sigma^2 \over 2}

  = \sqrt { ln(2)^2 E(x^2) \over 2}

So:

resn = log2(\sqrt { ln(2)^2 E(x^2) \over 2})

     = 0.5 * log2(ln(2)^2 E(x^2) \over 2})

     = 0.5 * log2(ln(2)^2 e \over {2 * blocksize}})
