Making the MDCT produce interleaved data
[opus.git] / libcelt / entdec.c
index e5617ba..47d710b 100644 (file)
+/* Copyright (c) 2001-2011 Timothy B. Terriberry
+   Copyright (c) 2008-2009 Xiph.Org Foundation */
+/*
+   Redistribution and use in source and binary forms, with or without
+   modification, are permitted provided that the following conditions
+   are met:
+
+   - Redistributions of source code must retain the above copyright
+   notice, this list of conditions and the following disclaimer.
+
+   - Redistributions in binary form must reproduce the above copyright
+   notice, this list of conditions and the following disclaimer in the
+   documentation and/or other materials provided with the distribution.
+
+   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+   ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
+   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
 #include <stddef.h>
+#include "os_support.h"
+#include "arch.h"
 #include "entdec.h"
+#include "mfrngcod.h"
 
+/*A range decoder.
+  This is an entropy decoder based upon \cite{Mar79}, which is itself a
+   rediscovery of the FIFO arithmetic code introduced by \cite{Pas76}.
+  It is very similar to arithmetic encoding, except that encoding is done with
+   digits in any base, instead of with bits, and so it is faster when using
+   larger bases (i.e.: a byte).
+  The author claims an average waste of $\frac{1}{2}\log_b(2b)$ bits, where $b$
+   is the base, longer than the theoretical optimum, but to my knowledge there
+   is no published justification for this claim.
+  This only seems true when using near-infinite precision arithmetic so that
+   the process is carried out with no rounding errors.
 
+  IBM (the author's employer) never sought to patent the idea, and to my
+   knowledge the algorithm is unencumbered by any patents, though its
+   performance is very competitive with proprietary arithmetic coding.
+  The two are based on very similar ideas, however.
+  An excellent description of implementation details is available at
+   http://www.arturocampos.com/ac_range.html
+  A recent work \cite{MNW98} which proposes several changes to arithmetic
+   encoding for efficiency actually re-discovers many of the principles
+   behind range encoding, and presents a good theoretical analysis of them.
 
-void ec_byte_readinit(ec_byte_buffer *_b,unsigned char *_buf,long _bytes){
-  _b->buf=_b->ptr=_buf;
-  _b->storage=_bytes;
+  End of stream is handled by writing out the smallest number of bits that
+   ensures that the stream will be correctly decoded regardless of the value of
+   any subsequent bits.
+  ec_tell() can be used to determine how many bits were needed to decode
+   all the symbols thus far; other data can be packed in the remaining bits of
+   the input buffer.
+  @PHDTHESIS{Pas76,
+    author="Richard Clark Pasco",
+    title="Source coding algorithms for fast data compression",
+    school="Dept. of Electrical Engineering, Stanford University",
+    address="Stanford, CA",
+    month=May,
+    year=1976
+  }
+  @INPROCEEDINGS{Mar79,
+   author="Martin, G.N.N.",
+   title="Range encoding: an algorithm for removing redundancy from a digitised
+    message",
+   booktitle="Video & Data Recording Conference",
+   year=1979,
+   address="Southampton",
+   month=Jul
+  }
+  @ARTICLE{MNW98,
+   author="Alistair Moffat and Radford Neal and Ian H. Witten",
+   title="Arithmetic Coding Revisited",
+   journal="{ACM} Transactions on Information Systems",
+   year=1998,
+   volume=16,
+   number=3,
+   pages="256--294",
+   month=Jul,
+   URL="http://www.stanford.edu/class/ee398/handouts/papers/Moffat98ArithmCoding.pdf"
+  }*/
+
+static int ec_read_byte(ec_dec *_this){
+  return _this->offs<_this->storage?_this->buf[_this->offs++]:0;
 }
 
-int ec_byte_look1(ec_byte_buffer *_b){
-  ptrdiff_t endbyte;
-  endbyte=_b->ptr-_b->buf;
-  if(endbyte>=_b->storage)return -1;
-  else return _b->ptr[0];
+static int ec_read_byte_from_end(ec_dec *_this){
+  return _this->end_offs<_this->storage?
+   _this->buf[_this->storage-++(_this->end_offs)]:0;
 }
 
-int ec_byte_look4(ec_byte_buffer *_b,ec_uint32 *_val){
-  ptrdiff_t endbyte;
-  endbyte=_b->ptr-_b->buf;
-  if(endbyte+4>_b->storage){
-    if(endbyte<_b->storage){
-      *_val=_b->ptr[0];
-      endbyte++;
-      if(endbyte<_b->storage){
-        *_val|=(ec_uint32)_b->ptr[1]<<8;
-        endbyte++;
-        if(endbyte<_b->storage)*_val|=(ec_uint32)_b->ptr[2]<<16;
-      }
-    }
-    return -1;
-  }
-  else{
-    *_val=_b->ptr[0];
-    *_val|=(ec_uint32)_b->ptr[1]<<8;
-    *_val|=(ec_uint32)_b->ptr[2]<<16;
-    *_val|=(ec_uint32)_b->ptr[3]<<24;
+/*Normalizes the contents of val and rng so that rng lies entirely in the
+   high-order symbol.*/
+static void ec_dec_normalize(ec_dec *_this){
+  /*If the range is too small, rescale it and input some bits.*/
+  while(_this->rng<=EC_CODE_BOT){
+    int sym;
+    _this->nbits_total+=EC_SYM_BITS;
+    _this->rng<<=EC_SYM_BITS;
+    /*Use up the remaining bits from our last symbol.*/
+    sym=_this->rem;
+    /*Read the next value from the input.*/
+    _this->rem=ec_read_byte(_this);
+    /*Take the rest of the bits we need from this new symbol.*/
+    sym=(sym<<EC_SYM_BITS|_this->rem)>>EC_SYM_BITS-EC_CODE_EXTRA;
+    /*And subtract them from val, capped to be less than EC_CODE_TOP.*/
+    _this->val=(_this->val<<EC_SYM_BITS)+(EC_SYM_MAX&~sym)&EC_CODE_TOP-1;
   }
-  return 0;
 }
 
-void ec_byte_adv1(ec_byte_buffer *_b){
-  _b->ptr++;
+void ec_dec_init(ec_dec *_this,unsigned char *_buf,opus_uint32 _storage){
+  _this->buf=_buf;
+  _this->storage=_storage;
+  _this->end_offs=0;
+  _this->end_window=0;
+  _this->nend_bits=0;
+  _this->offs=0;
+  _this->rng=1U<<EC_CODE_EXTRA;
+  _this->rem=ec_read_byte(_this);
+  _this->val=_this->rng-1-(_this->rem>>EC_SYM_BITS-EC_CODE_EXTRA);
+  _this->error=0;
+  /*Normalize the interval.*/
+  ec_dec_normalize(_this);
+  /*This is the offset from which ec_tell() will subtract partial bits.
+    This must be after the initial ec_dec_normalize(), or you will have to
+     compensate for the bits that are read there.*/
+  _this->nbits_total=EC_CODE_BITS+1;
 }
 
-void ec_byte_adv4(ec_byte_buffer *_b){
-  _b->ptr+=4;
+unsigned ec_decode(ec_dec *_this,unsigned _ft){
+  unsigned s;
+  _this->ext=_this->rng/_ft;
+  s=(unsigned)(_this->val/_this->ext);
+  return _ft-EC_MINI(s+1,_ft);
 }
 
-int ec_byte_read1(ec_byte_buffer *_b){
-  ptrdiff_t endbyte;
-  endbyte=_b->ptr-_b->buf;
-  if(endbyte>=_b->storage)return -1;
-  else return *(_b->ptr++);
+unsigned ec_decode_bin(ec_dec *_this,unsigned _bits){
+   unsigned s;
+   _this->ext=_this->rng>>_bits;
+   s=(unsigned)(_this->val/_this->ext);
+   return (1<<_bits)-EC_MINI(s+1,1<<_bits);
 }
 
-int ec_byte_read4(ec_byte_buffer *_b,ec_uint32 *_val){
-  unsigned char *end;
-  end=_b->buf+_b->storage;
-  if(_b->ptr+4>end){
-    if(_b->ptr<end){
-      *_val=*(_b->ptr++);
-      if(_b->ptr<end){
-        *_val|=(ec_uint32)*(_b->ptr++)<<8;
-        if(_b->ptr<end)*_val|=(ec_uint32)*(_b->ptr++)<<16;
-      }
-    }
-    return -1;
-  }
-  else{
-    *_val=(*_b->ptr++);
-    *_val|=(ec_uint32)*(_b->ptr++)<<8;
-    *_val|=(ec_uint32)*(_b->ptr++)<<16;
-    *_val|=(ec_uint32)*(_b->ptr++)<<24;
-  }
-  return 0;
+void ec_dec_update(ec_dec *_this,unsigned _fl,unsigned _fh,unsigned _ft){
+  opus_uint32 s;
+  s=IMUL32(_this->ext,_ft-_fh);
+  _this->val-=s;
+  _this->rng=_fl>0?IMUL32(_this->ext,_fh-_fl):_this->rng-s;
+  ec_dec_normalize(_this);
 }
 
-
-
-ec_uint32 ec_dec_bits(ec_dec *_this,int _ftb){
-  ec_uint32 t;
-  unsigned  s;
-  unsigned  ft;
-  t=0;
-  while(_ftb>EC_UNIT_BITS){
-    s=ec_decode(_this,EC_UNIT_MASK+1);
-    ec_dec_update(_this,s,s+1,EC_UNIT_MASK+1);
-    t=t<<EC_UNIT_BITS|s;
-    _ftb-=EC_UNIT_BITS;
-  }
-  ft=1U<<_ftb;
-  s=ec_decode(_this,ft);
-  ec_dec_update(_this,s,s+1,ft);
-  t=t<<_ftb|s;
-  return t;
+/*The probability of having a "one" is 1/(1<<_logp).*/
+int ec_dec_bit_logp(ec_dec *_this,unsigned _logp){
+  opus_uint32 r;
+  opus_uint32 d;
+  opus_uint32 s;
+  int         ret;
+  r=_this->rng;
+  d=_this->val;
+  s=r>>_logp;
+  ret=d<s;
+  if(!ret)_this->val=d-s;
+  _this->rng=ret?s:r-s;
+  ec_dec_normalize(_this);
+  return ret;
 }
 
-ec_uint64 ec_dec_bits64(ec_dec *_this,int _ftb){
-  ec_uint32 t;
-  if(_ftb>32){
-    t=ec_dec_bits(_this,_ftb-32);
-    _ftb=32;
+int ec_dec_icdf(ec_dec *_this,const unsigned char *_icdf,unsigned _ftb){
+  opus_uint32 r;
+  opus_uint32 d;
+  opus_uint32 s;
+  opus_uint32 t;
+  int         ret;
+  s=_this->rng;
+  d=_this->val;
+  r=s>>_ftb;
+  ret=-1;
+  do{
+    t=s;
+    s=IMUL32(r,_icdf[++ret]);
   }
-  else t=0;
-  return (ec_uint64)t<<32|ec_dec_bits(_this,_ftb);
+  while(d<s);
+  _this->val=d-s;
+  _this->rng=t-s;
+  ec_dec_normalize(_this);
+  return ret;
 }
 
-ec_uint32 ec_dec_uint(ec_dec *_this,ec_uint32 _ft){
-  ec_uint32 mask;
-  ec_uint32 t;
-  unsigned  ft;
-  unsigned  s;
-  int       ftb;
-  t=0;
+opus_uint32 ec_dec_uint(ec_dec *_this,opus_uint32 _ft){
+  unsigned ft;
+  unsigned s;
+  int      ftb;
+  /*In order to optimize EC_ILOG(), it is undefined for the value 0.*/
+  celt_assert(_ft>1);
   _ft--;
   ftb=EC_ILOG(_ft);
-  while(ftb>EC_UNIT_BITS){
-    ftb-=EC_UNIT_BITS;
+  if(ftb>EC_UINT_BITS){
+    opus_uint32 t;
+    ftb-=EC_UINT_BITS;
     ft=(unsigned)(_ft>>ftb)+1;
     s=ec_decode(_this,ft);
     ec_dec_update(_this,s,s+1,ft);
-    t=t<<EC_UNIT_BITS|s;
-    if(s<ft-1)return t<<ftb|ec_dec_bits(_this,ftb);
-    mask=((ec_uint32)1<<ftb)-1;
-    _ft=_ft&mask;
+    t=(opus_uint32)s<<ftb|ec_dec_bits(_this,ftb);
+    if(t<=_ft)return t;
+    _this->error=1;
+    return _ft;
+  }
+  else{
+    _ft++;
+    s=ec_decode(_this,(unsigned)_ft);
+    ec_dec_update(_this,s,s+1,(unsigned)_ft);
+    return s;
   }
-  _ft++;
-  s=ec_decode(_this,(unsigned)_ft);
-  ec_dec_update(_this,s,s+1,(unsigned)_ft);
-  t=t<<ftb|s;
-  return t;
 }
 
-ec_uint64 ec_dec_uint64(ec_dec *_this,ec_uint64 _ft){
-  ec_uint64 mask;
-  ec_uint64 t;
-  unsigned  ft;
-  unsigned  s;
-  int       ftb;
-  t=0;
-  _ft--;
-  ftb=EC_ILOG64(_ft);
-  while(ftb>EC_UNIT_BITS){
-    ftb-=EC_UNIT_BITS;
-    ft=(unsigned)(_ft>>ftb)+1;
-    s=ec_decode(_this,ft);
-    ec_dec_update(_this,s,s+1,ft);
-    t=t<<EC_UNIT_BITS|s;
-    if(s<ft-1)return t<<ftb|ec_dec_bits64(_this,ftb);
-    mask=((ec_uint64)1<<ftb)-1;
-    _ft=_ft&mask;
+opus_uint32 ec_dec_bits(ec_dec *_this,unsigned _bits){
+  ec_window   window;
+  int         available;
+  opus_uint32 ret;
+  window=_this->end_window;
+  available=_this->nend_bits;
+  if(available<_bits){
+    do{
+      window|=(ec_window)ec_read_byte_from_end(_this)<<available;
+      available+=EC_SYM_BITS;
+    }
+    while(available<=EC_WINDOW_SIZE-EC_SYM_BITS);
   }
-  _ft++;
-  s=ec_decode(_this,(unsigned)_ft);
-  ec_dec_update(_this,s,s+1,(unsigned)_ft);
-  t=t<<ftb|s;
-  return t;
+  ret=(opus_uint32)window&((opus_uint32)1<<_bits)-1;
+  window>>=_bits;
+  available-=_bits;
+  _this->end_window=window;
+  _this->nend_bits=available;
+  _this->nbits_total+=_bits;
+  return ret;
 }