Add a generic CDF decoding routine.

This decodes a value encoded with ec_encode_bin() without using any
 divisions.
It is only meant for small alphabets.
If a symbol can take on a large number of possible values, a binary
 search would be better.

This patch also converts spread_decision to use it, since it is
 faster and introduces less rounding error to encode a single
 decision for the entire value than to encode it a bit at a time.

libcelt/celt.c

@@ -54,6 +54,7 @@
 #include "plc.h"
 
 static const int trim_cdf[12] = {0, 2, 4, 9, 19, 41, 87, 109, 119, 124, 126, 128};
+static const int spread_cdf[5] = {0, 7, 9, 30, 32};
 
 #define COMBFILTER_MAXPERIOD 1024
 #define COMBFILTER_MINPERIOD 16
@@ -957,9 +958,8 @@ int celt_encode_with_ec_float(CELTEncoder * restrict st, const celt_sig * pcm, i
       st->spread_decision = spreading_decision(st->mode, X, &st->tonal_average, st->spread_decision, effEnd, C, M);
    }
    /* Probs: NONE: 21.875%, LIGHT: 6.25%, NORMAL: 65.625%, AGGRESSIVE: 6.25% */
-   ec_enc_bit_prob(enc, st->spread_decision>>1, 47104);
+   ec_encode_bin(enc, spread_cdf[st->spread_decision],
-   ec_enc_bit_prob(enc, st->spread_decision&1,
+         spread_cdf[st->spread_decision+1], 5);
-         (st->spread_decision>>1) ? 5699 : 14564);
 
    ALLOC(offsets, st->mode->nbEBands, int);
 
@@ -1823,8 +1823,7 @@ int celt_decode_with_ec_float(CELTDecoder * restrict st, const unsigned char *da
    ALLOC(tf_res, st->mode->nbEBands, int);
    tf_decode(st->start, st->end, C, isTransient, tf_res, LM, dec);
 
-   spread_decision = ec_dec_bit_prob(dec, 47104)<<1;
+   spread_decision = ec_dec_cdf(dec, spread_cdf, 5);
-   spread_decision |= ec_dec_bit_prob(dec, (spread_decision>>1) ? 5699 : 14564);
 
    ALLOC(pulses, st->mode->nbEBands, int);
    ALLOC(offsets, st->mode->nbEBands, int);
@@ -1844,14 +1843,7 @@ int celt_decode_with_ec_float(CELTDecoder * restrict st, const unsigned char *da
    }
 
    ALLOC(fine_quant, st->mode->nbEBands, int);
-   {
+   alloc_trim = ec_dec_cdf(dec, trim_cdf, 7);
-      int fl;
-      alloc_trim = 0;
-      fl = ec_decode_bin(dec, 7);
-      while (trim_cdf[alloc_trim+1] <= fl)
-         alloc_trim++;
-      ec_dec_update(dec, trim_cdf[alloc_trim], trim_cdf[alloc_trim+1], 128);
-   }
 
    if (C==2)
    {

libcelt/entdec.h

@@ -110,6 +110,15 @@ ec_uint32 ec_dec_bits(ec_dec *_this,unsigned _ftb);
        This must be at least one, and no more than 2**32-1.
   Return: The decoded bits.*/
 ec_uint32 ec_dec_uint(ec_dec *_this,ec_uint32 _ft);
+/*Decodes a symbol given its CDF.
+  No call to ec_dec_update() is necessary after this call.
+  _cdf: The CDF, such that symbol s falls in the range [_cdf[s],_cdf[s+1]).
+        The first value must be 0, the last value must be (1<<_ftb), and the
+         values must be monotonicly non-decreasing.
+  _ftb: The number of bits of precision in the cumulative distribution.
+  Return: The decoded symbol s, which must have been encoded with
+   ec_encode_bin(enc,_cdf[s],_cdf[s+1],_ftb).*/
+int ec_dec_cdf(ec_dec *_this,const int *_cdf,unsigned _ftb);
 
 /* Decode a bit that has a _prob/65536 probability of being a one */
 int ec_dec_bit_prob(ec_dec *_this,unsigned _prob);

libcelt/rangedec.c

@@ -203,6 +203,28 @@ int ec_dec_bit_prob(ec_dec *_this,unsigned _prob){
   return val;
 }
 
+int ec_dec_cdf(ec_dec *_this,const int *_cdf,unsigned _ftb){
+  ec_uint32 r;
+  ec_uint32 d;
+  ec_uint32 s;
+  ec_uint32 t;
+  int       val;
+  r=_this->rng>>_ftb;
+  d=_this->dif;
+  _cdf++;
+  val=0;
+  t=0;
+  s=IMUL32(r,(1<<_ftb)-_cdf[0]);
+  while(d<=s){
+    t=s;
+    s=IMUL32(r,(1<<_ftb)-_cdf[++val]);
+  }
+  _this->dif=d-s;
+  _this->rng=(val?t:_this->rng)-s;
+  ec_dec_normalize(_this);
+  return val;
+}
+
 ec_uint32 ec_dec_tell(ec_dec *_this,int _b){
   ec_uint32 r;
   int       l;