1*a58d3d2aSXin Li /* Copyright (c) 2001-2011 Timothy B. Terriberry
2*a58d3d2aSXin Li Copyright (c) 2008-2009 Xiph.Org Foundation */
3*a58d3d2aSXin Li /*
4*a58d3d2aSXin Li Redistribution and use in source and binary forms, with or without
5*a58d3d2aSXin Li modification, are permitted provided that the following conditions
6*a58d3d2aSXin Li are met:
7*a58d3d2aSXin Li
8*a58d3d2aSXin Li - Redistributions of source code must retain the above copyright
9*a58d3d2aSXin Li notice, this list of conditions and the following disclaimer.
10*a58d3d2aSXin Li
11*a58d3d2aSXin Li - Redistributions in binary form must reproduce the above copyright
12*a58d3d2aSXin Li notice, this list of conditions and the following disclaimer in the
13*a58d3d2aSXin Li documentation and/or other materials provided with the distribution.
14*a58d3d2aSXin Li
15*a58d3d2aSXin Li THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
16*a58d3d2aSXin Li ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
17*a58d3d2aSXin Li LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
18*a58d3d2aSXin Li A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
19*a58d3d2aSXin Li OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
20*a58d3d2aSXin Li EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
21*a58d3d2aSXin Li PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
22*a58d3d2aSXin Li PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
23*a58d3d2aSXin Li LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
24*a58d3d2aSXin Li NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
25*a58d3d2aSXin Li SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26*a58d3d2aSXin Li */
27*a58d3d2aSXin Li
28*a58d3d2aSXin Li #ifdef HAVE_CONFIG_H
29*a58d3d2aSXin Li #include "config.h"
30*a58d3d2aSXin Li #endif
31*a58d3d2aSXin Li
32*a58d3d2aSXin Li #include <stddef.h>
33*a58d3d2aSXin Li #include "os_support.h"
34*a58d3d2aSXin Li #include "arch.h"
35*a58d3d2aSXin Li #include "entdec.h"
36*a58d3d2aSXin Li #include "mfrngcod.h"
37*a58d3d2aSXin Li
38*a58d3d2aSXin Li /*A range decoder.
39*a58d3d2aSXin Li This is an entropy decoder based upon \cite{Mar79}, which is itself a
40*a58d3d2aSXin Li rediscovery of the FIFO arithmetic code introduced by \cite{Pas76}.
41*a58d3d2aSXin Li It is very similar to arithmetic encoding, except that encoding is done with
42*a58d3d2aSXin Li digits in any base, instead of with bits, and so it is faster when using
43*a58d3d2aSXin Li larger bases (i.e.: a byte).
44*a58d3d2aSXin Li The author claims an average waste of $\frac{1}{2}\log_b(2b)$ bits, where $b$
45*a58d3d2aSXin Li is the base, longer than the theoretical optimum, but to my knowledge there
46*a58d3d2aSXin Li is no published justification for this claim.
47*a58d3d2aSXin Li This only seems true when using near-infinite precision arithmetic so that
48*a58d3d2aSXin Li the process is carried out with no rounding errors.
49*a58d3d2aSXin Li
50*a58d3d2aSXin Li An excellent description of implementation details is available at
51*a58d3d2aSXin Li http://www.arturocampos.com/ac_range.html
52*a58d3d2aSXin Li A recent work \cite{MNW98} which proposes several changes to arithmetic
53*a58d3d2aSXin Li encoding for efficiency actually re-discovers many of the principles
54*a58d3d2aSXin Li behind range encoding, and presents a good theoretical analysis of them.
55*a58d3d2aSXin Li
56*a58d3d2aSXin Li End of stream is handled by writing out the smallest number of bits that
57*a58d3d2aSXin Li ensures that the stream will be correctly decoded regardless of the value of
58*a58d3d2aSXin Li any subsequent bits.
59*a58d3d2aSXin Li ec_tell() can be used to determine how many bits were needed to decode
60*a58d3d2aSXin Li all the symbols thus far; other data can be packed in the remaining bits of
61*a58d3d2aSXin Li the input buffer.
62*a58d3d2aSXin Li @PHDTHESIS{Pas76,
63*a58d3d2aSXin Li author="Richard Clark Pasco",
64*a58d3d2aSXin Li title="Source coding algorithms for fast data compression",
65*a58d3d2aSXin Li school="Dept. of Electrical Engineering, Stanford University",
66*a58d3d2aSXin Li address="Stanford, CA",
67*a58d3d2aSXin Li month=May,
68*a58d3d2aSXin Li year=1976
69*a58d3d2aSXin Li }
70*a58d3d2aSXin Li @INPROCEEDINGS{Mar79,
71*a58d3d2aSXin Li author="Martin, G.N.N.",
72*a58d3d2aSXin Li title="Range encoding: an algorithm for removing redundancy from a digitised
73*a58d3d2aSXin Li message",
74*a58d3d2aSXin Li booktitle="Video & Data Recording Conference",
75*a58d3d2aSXin Li year=1979,
76*a58d3d2aSXin Li address="Southampton",
77*a58d3d2aSXin Li month=Jul
78*a58d3d2aSXin Li }
79*a58d3d2aSXin Li @ARTICLE{MNW98,
80*a58d3d2aSXin Li author="Alistair Moffat and Radford Neal and Ian H. Witten",
81*a58d3d2aSXin Li title="Arithmetic Coding Revisited",
82*a58d3d2aSXin Li journal="{ACM} Transactions on Information Systems",
83*a58d3d2aSXin Li year=1998,
84*a58d3d2aSXin Li volume=16,
85*a58d3d2aSXin Li number=3,
86*a58d3d2aSXin Li pages="256--294",
87*a58d3d2aSXin Li month=Jul,
88*a58d3d2aSXin Li URL="http://www.stanford.edu/class/ee398a/handouts/papers/Moffat98ArithmCoding.pdf"
89*a58d3d2aSXin Li }*/
90*a58d3d2aSXin Li
ec_read_byte(ec_dec * _this)91*a58d3d2aSXin Li static int ec_read_byte(ec_dec *_this){
92*a58d3d2aSXin Li return _this->offs<_this->storage?_this->buf[_this->offs++]:0;
93*a58d3d2aSXin Li }
94*a58d3d2aSXin Li
ec_read_byte_from_end(ec_dec * _this)95*a58d3d2aSXin Li static int ec_read_byte_from_end(ec_dec *_this){
96*a58d3d2aSXin Li return _this->end_offs<_this->storage?
97*a58d3d2aSXin Li _this->buf[_this->storage-++(_this->end_offs)]:0;
98*a58d3d2aSXin Li }
99*a58d3d2aSXin Li
100*a58d3d2aSXin Li /*Normalizes the contents of val and rng so that rng lies entirely in the
101*a58d3d2aSXin Li high-order symbol.*/
ec_dec_normalize(ec_dec * _this)102*a58d3d2aSXin Li static void ec_dec_normalize(ec_dec *_this){
103*a58d3d2aSXin Li /*If the range is too small, rescale it and input some bits.*/
104*a58d3d2aSXin Li while(_this->rng<=EC_CODE_BOT){
105*a58d3d2aSXin Li int sym;
106*a58d3d2aSXin Li _this->nbits_total+=EC_SYM_BITS;
107*a58d3d2aSXin Li _this->rng<<=EC_SYM_BITS;
108*a58d3d2aSXin Li /*Use up the remaining bits from our last symbol.*/
109*a58d3d2aSXin Li sym=_this->rem;
110*a58d3d2aSXin Li /*Read the next value from the input.*/
111*a58d3d2aSXin Li _this->rem=ec_read_byte(_this);
112*a58d3d2aSXin Li /*Take the rest of the bits we need from this new symbol.*/
113*a58d3d2aSXin Li sym=(sym<<EC_SYM_BITS|_this->rem)>>(EC_SYM_BITS-EC_CODE_EXTRA);
114*a58d3d2aSXin Li /*And subtract them from val, capped to be less than EC_CODE_TOP.*/
115*a58d3d2aSXin Li _this->val=((_this->val<<EC_SYM_BITS)+(EC_SYM_MAX&~sym))&(EC_CODE_TOP-1);
116*a58d3d2aSXin Li }
117*a58d3d2aSXin Li }
118*a58d3d2aSXin Li
ec_dec_init(ec_dec * _this,unsigned char * _buf,opus_uint32 _storage)119*a58d3d2aSXin Li void ec_dec_init(ec_dec *_this,unsigned char *_buf,opus_uint32 _storage){
120*a58d3d2aSXin Li _this->buf=_buf;
121*a58d3d2aSXin Li _this->storage=_storage;
122*a58d3d2aSXin Li _this->end_offs=0;
123*a58d3d2aSXin Li _this->end_window=0;
124*a58d3d2aSXin Li _this->nend_bits=0;
125*a58d3d2aSXin Li /*This is the offset from which ec_tell() will subtract partial bits.
126*a58d3d2aSXin Li The final value after the ec_dec_normalize() call will be the same as in
127*a58d3d2aSXin Li the encoder, but we have to compensate for the bits that are added there.*/
128*a58d3d2aSXin Li _this->nbits_total=EC_CODE_BITS+1
129*a58d3d2aSXin Li -((EC_CODE_BITS-EC_CODE_EXTRA)/EC_SYM_BITS)*EC_SYM_BITS;
130*a58d3d2aSXin Li _this->offs=0;
131*a58d3d2aSXin Li _this->rng=1U<<EC_CODE_EXTRA;
132*a58d3d2aSXin Li _this->rem=ec_read_byte(_this);
133*a58d3d2aSXin Li _this->val=_this->rng-1-(_this->rem>>(EC_SYM_BITS-EC_CODE_EXTRA));
134*a58d3d2aSXin Li _this->error=0;
135*a58d3d2aSXin Li /*Normalize the interval.*/
136*a58d3d2aSXin Li ec_dec_normalize(_this);
137*a58d3d2aSXin Li }
138*a58d3d2aSXin Li
ec_decode(ec_dec * _this,unsigned _ft)139*a58d3d2aSXin Li unsigned ec_decode(ec_dec *_this,unsigned _ft){
140*a58d3d2aSXin Li unsigned s;
141*a58d3d2aSXin Li _this->ext=celt_udiv(_this->rng,_ft);
142*a58d3d2aSXin Li s=(unsigned)(_this->val/_this->ext);
143*a58d3d2aSXin Li return _ft-EC_MINI(s+1,_ft);
144*a58d3d2aSXin Li }
145*a58d3d2aSXin Li
ec_decode_bin(ec_dec * _this,unsigned _bits)146*a58d3d2aSXin Li unsigned ec_decode_bin(ec_dec *_this,unsigned _bits){
147*a58d3d2aSXin Li unsigned s;
148*a58d3d2aSXin Li _this->ext=_this->rng>>_bits;
149*a58d3d2aSXin Li s=(unsigned)(_this->val/_this->ext);
150*a58d3d2aSXin Li return (1U<<_bits)-EC_MINI(s+1U,1U<<_bits);
151*a58d3d2aSXin Li }
152*a58d3d2aSXin Li
ec_dec_update(ec_dec * _this,unsigned _fl,unsigned _fh,unsigned _ft)153*a58d3d2aSXin Li void ec_dec_update(ec_dec *_this,unsigned _fl,unsigned _fh,unsigned _ft){
154*a58d3d2aSXin Li opus_uint32 s;
155*a58d3d2aSXin Li s=IMUL32(_this->ext,_ft-_fh);
156*a58d3d2aSXin Li _this->val-=s;
157*a58d3d2aSXin Li _this->rng=_fl>0?IMUL32(_this->ext,_fh-_fl):_this->rng-s;
158*a58d3d2aSXin Li ec_dec_normalize(_this);
159*a58d3d2aSXin Li }
160*a58d3d2aSXin Li
161*a58d3d2aSXin Li /*The probability of having a "one" is 1/(1<<_logp).*/
ec_dec_bit_logp(ec_dec * _this,unsigned _logp)162*a58d3d2aSXin Li int ec_dec_bit_logp(ec_dec *_this,unsigned _logp){
163*a58d3d2aSXin Li opus_uint32 r;
164*a58d3d2aSXin Li opus_uint32 d;
165*a58d3d2aSXin Li opus_uint32 s;
166*a58d3d2aSXin Li int ret;
167*a58d3d2aSXin Li r=_this->rng;
168*a58d3d2aSXin Li d=_this->val;
169*a58d3d2aSXin Li s=r>>_logp;
170*a58d3d2aSXin Li ret=d<s;
171*a58d3d2aSXin Li if(!ret)_this->val=d-s;
172*a58d3d2aSXin Li _this->rng=ret?s:r-s;
173*a58d3d2aSXin Li ec_dec_normalize(_this);
174*a58d3d2aSXin Li return ret;
175*a58d3d2aSXin Li }
176*a58d3d2aSXin Li
ec_dec_icdf(ec_dec * _this,const unsigned char * _icdf,unsigned _ftb)177*a58d3d2aSXin Li int ec_dec_icdf(ec_dec *_this,const unsigned char *_icdf,unsigned _ftb){
178*a58d3d2aSXin Li opus_uint32 r;
179*a58d3d2aSXin Li opus_uint32 d;
180*a58d3d2aSXin Li opus_uint32 s;
181*a58d3d2aSXin Li opus_uint32 t;
182*a58d3d2aSXin Li int ret;
183*a58d3d2aSXin Li s=_this->rng;
184*a58d3d2aSXin Li d=_this->val;
185*a58d3d2aSXin Li r=s>>_ftb;
186*a58d3d2aSXin Li ret=-1;
187*a58d3d2aSXin Li do{
188*a58d3d2aSXin Li t=s;
189*a58d3d2aSXin Li s=IMUL32(r,_icdf[++ret]);
190*a58d3d2aSXin Li }
191*a58d3d2aSXin Li while(d<s);
192*a58d3d2aSXin Li _this->val=d-s;
193*a58d3d2aSXin Li _this->rng=t-s;
194*a58d3d2aSXin Li ec_dec_normalize(_this);
195*a58d3d2aSXin Li return ret;
196*a58d3d2aSXin Li }
197*a58d3d2aSXin Li
ec_dec_icdf16(ec_dec * _this,const opus_uint16 * _icdf,unsigned _ftb)198*a58d3d2aSXin Li int ec_dec_icdf16(ec_dec *_this,const opus_uint16 *_icdf,unsigned _ftb){
199*a58d3d2aSXin Li opus_uint32 r;
200*a58d3d2aSXin Li opus_uint32 d;
201*a58d3d2aSXin Li opus_uint32 s;
202*a58d3d2aSXin Li opus_uint32 t;
203*a58d3d2aSXin Li int ret;
204*a58d3d2aSXin Li s=_this->rng;
205*a58d3d2aSXin Li d=_this->val;
206*a58d3d2aSXin Li r=s>>_ftb;
207*a58d3d2aSXin Li ret=-1;
208*a58d3d2aSXin Li do{
209*a58d3d2aSXin Li t=s;
210*a58d3d2aSXin Li s=IMUL32(r,_icdf[++ret]);
211*a58d3d2aSXin Li }
212*a58d3d2aSXin Li while(d<s);
213*a58d3d2aSXin Li _this->val=d-s;
214*a58d3d2aSXin Li _this->rng=t-s;
215*a58d3d2aSXin Li ec_dec_normalize(_this);
216*a58d3d2aSXin Li return ret;
217*a58d3d2aSXin Li }
218*a58d3d2aSXin Li
ec_dec_uint(ec_dec * _this,opus_uint32 _ft)219*a58d3d2aSXin Li opus_uint32 ec_dec_uint(ec_dec *_this,opus_uint32 _ft){
220*a58d3d2aSXin Li unsigned ft;
221*a58d3d2aSXin Li unsigned s;
222*a58d3d2aSXin Li int ftb;
223*a58d3d2aSXin Li /*In order to optimize EC_ILOG(), it is undefined for the value 0.*/
224*a58d3d2aSXin Li celt_assert(_ft>1);
225*a58d3d2aSXin Li _ft--;
226*a58d3d2aSXin Li ftb=EC_ILOG(_ft);
227*a58d3d2aSXin Li if(ftb>EC_UINT_BITS){
228*a58d3d2aSXin Li opus_uint32 t;
229*a58d3d2aSXin Li ftb-=EC_UINT_BITS;
230*a58d3d2aSXin Li ft=(unsigned)(_ft>>ftb)+1;
231*a58d3d2aSXin Li s=ec_decode(_this,ft);
232*a58d3d2aSXin Li ec_dec_update(_this,s,s+1,ft);
233*a58d3d2aSXin Li t=(opus_uint32)s<<ftb|ec_dec_bits(_this,ftb);
234*a58d3d2aSXin Li if(t<=_ft)return t;
235*a58d3d2aSXin Li _this->error=1;
236*a58d3d2aSXin Li return _ft;
237*a58d3d2aSXin Li }
238*a58d3d2aSXin Li else{
239*a58d3d2aSXin Li _ft++;
240*a58d3d2aSXin Li s=ec_decode(_this,(unsigned)_ft);
241*a58d3d2aSXin Li ec_dec_update(_this,s,s+1,(unsigned)_ft);
242*a58d3d2aSXin Li return s;
243*a58d3d2aSXin Li }
244*a58d3d2aSXin Li }
245*a58d3d2aSXin Li
ec_dec_bits(ec_dec * _this,unsigned _bits)246*a58d3d2aSXin Li opus_uint32 ec_dec_bits(ec_dec *_this,unsigned _bits){
247*a58d3d2aSXin Li ec_window window;
248*a58d3d2aSXin Li int available;
249*a58d3d2aSXin Li opus_uint32 ret;
250*a58d3d2aSXin Li window=_this->end_window;
251*a58d3d2aSXin Li available=_this->nend_bits;
252*a58d3d2aSXin Li if((unsigned)available<_bits){
253*a58d3d2aSXin Li do{
254*a58d3d2aSXin Li window|=(ec_window)ec_read_byte_from_end(_this)<<available;
255*a58d3d2aSXin Li available+=EC_SYM_BITS;
256*a58d3d2aSXin Li }
257*a58d3d2aSXin Li while(available<=EC_WINDOW_SIZE-EC_SYM_BITS);
258*a58d3d2aSXin Li }
259*a58d3d2aSXin Li ret=(opus_uint32)window&(((opus_uint32)1<<_bits)-1U);
260*a58d3d2aSXin Li window>>=_bits;
261*a58d3d2aSXin Li available-=_bits;
262*a58d3d2aSXin Li _this->end_window=window;
263*a58d3d2aSXin Li _this->nend_bits=available;
264*a58d3d2aSXin Li _this->nbits_total+=_bits;
265*a58d3d2aSXin Li return ret;
266*a58d3d2aSXin Li }
267