License for the kiss-fft headers
[opus.git] / libcelt / rangedec.c
1 /* (C) 2001-2008 Timothy B. Terriberry
2    (C) 2008 Jean-Marc Valin */
3 /*
4    Redistribution and use in source and binary forms, with or without
5    modification, are permitted provided that the following conditions
6    are met:
7
8    - Redistributions of source code must retain the above copyright
9    notice, this list of conditions and the following disclaimer.
10
11    - Redistributions in binary form must reproduce the above copyright
12    notice, this list of conditions and the following disclaimer in the
13    documentation and/or other materials provided with the distribution.
14
15    - Neither the name of the Xiph.org Foundation nor the names of its
16    contributors may be used to endorse or promote products derived from
17    this software without specific prior written permission.
18
19    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
20    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
21    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
22    A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
23    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
24    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
25    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
26    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
27    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
28    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
29    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
30 */
31
32 #ifdef HAVE_CONFIG_H
33 #include "config.h"
34 #endif
35
36 #include "arch.h"
37 #include "entdec.h"
38 #include "mfrngcod.h"
39
40
41
42 /*A range decoder.
43   This is an entropy decoder based upon \cite{Mar79}, which is itself a
44    rediscovery of the FIFO arithmetic code introduced by \cite{Pas76}.
45   It is very similar to arithmetic encoding, except that encoding is done with
46    digits in any base, instead of with bits, and so it is faster when using
47    larger bases (i.e.: a byte).
48   The author claims an average waste of $\frac{1}{2}\log_b(2b)$ bits, where $b$
49    is the base, longer than the theoretical optimum, but to my knowledge there
50    is no published justification for this claim.
51   This only seems true when using near-infinite precision arithmetic so that
52    the process is carried out with no rounding errors.
53
54   IBM (the author's employer) never sought to patent the idea, and to my
55    knowledge the algorithm is unencumbered by any patents, though its
56    performance is very competitive with proprietary arithmetic coding.
57   The two are based on very similar ideas, however.
58   An excellent description of implementation details is available at
59    http://www.arturocampos.com/ac_range.html
60   A recent work \cite{MNW98} which proposes several changes to arithmetic
61    encoding for efficiency actually re-discovers many of the principles
62    behind range encoding, and presents a good theoretical analysis of them.
63
64   This coder handles the end of the stream in a slightly more graceful fashion
65    than most arithmetic or range coders.
66   Once the final symbol has been encoded, the coder selects the code word with
67    the shortest number of bits that still falls within the final interval.
68   This method is not novel.
69   Here, by the length of the code word, we refer to the number of bits until
70    its final 1.
71   Any trailing zeros may be discarded, since the encoder, once it runs out of
72    input, will pad its buffer with zeros.
73
74   But this means that no encoded stream would ever have any zero bytes at the
75    end.
76   Since there are some coded representations we cannot produce, it implies that
77    there is still some redundancy in the stream.
78   In this case, we can pick a special byte value, RSV1, and should the stream
79    end in a sequence of zeros, followed by the RSV1 byte, we can code the
80    zeros, and discard the RSV1 byte.
81   The decoder, knowing that the encoder would never produce a sequence of zeros
82    at the end, would then know to add in the RSV1 byte if it observed it.
83
84   Now, the encoder would never produce a stream that ended in a sequence of
85    zeros followed by a RSV1 byte.
86   So, if the stream ends in a non-empty sequence of zeros, followed by any
87    positive number of RSV1 bytes, the last RSV1 byte is discarded.
88   The decoder, if it encounters a stream that ends in non-empty sequence of
89    zeros followed by any non-negative number of RSV1 bytes, adds an additional
90    RSV1 byte to the stream.
91   With this strategy, every possible sequence of input bytes is transformed to
92    one that could actually be produced by the encoder.
93
94   The only question is what non-zero value to use for RSV1.
95   We select 0x80, since it has the nice property of producing the shortest
96    possible byte streams when using our strategy for selecting a number within
97    the final interval to encode.
98   Clearly if the shortest possible code word that falls within the interval has
99    its last one bit as the most significant bit of the final byte, and the
100    previous bytes were a non-empty sequence of zeros followed by a non-negative
101    number of 0x80 bytes, then the last byte would be discarded.
102   If the shortest code word is not so formed, then no other code word in the
103    interval would result in any more bytes being discarded.
104   Any longer code word would have an additional one bit somewhere, and so would
105    require at a minimum that that byte would be coded.
106   If the shortest code word has a 1 before the final one that is preventing the
107    stream from ending in a non-empty sequence of zeros followed by a
108    non-negative number of 0x80's, then there is no code word of the same length
109    which contains that bit as a zero.
110   If there were, then we could simply leave that bit a 1, and drop all the bits
111    after it without leaving the interval, thus producing a shorter code word.
112
113   In this case, RSV1 can only drop 1 bit off the final stream.
114   Other choices could lead to savings of up to 8 bits for particular streams,
115    but this would produce the odd situation that a stream with more non-zero
116    bits is actually encoded in fewer bytes.
117
118   @PHDTHESIS{Pas76,
119     author="Richard Clark Pasco",
120     title="Source coding algorithms for fast data compression",
121     school="Dept. of Electrical Engineering, Stanford University",
122     address="Stanford, CA",
123     month=May,
124     year=1976
125   }
126   @INPROCEEDINGS{Mar79,
127    author="Martin, G.N.N.",
128    title="Range encoding: an algorithm for removing redundancy from a digitised
129     message",
130    booktitle="Video & Data Recording Conference",
131    year=1979,
132    address="Southampton",
133    month=Jul
134   }
135   @ARTICLE{MNW98,
136    author="Alistair Moffat and Radford Neal and Ian H. Witten",
137    title="Arithmetic Coding Revisited",
138    journal="{ACM} Transactions on Information Systems",
139    year=1998,
140    volume=16,
141    number=3,
142    pages="256--294",
143    month=Jul,
144    URL="http://www.stanford.edu/class/ee398/handouts/papers/Moffat98ArithmCoding.pdf"
145   }*/
146
147
148 /*Gets the next byte of input.
149   After all the bytes in the current packet have been consumed, and the extra
150    end code returned if needed, this function will continue to return zero each
151    time it is called.
152   Return: The next byte of input.*/
153 static int ec_dec_in(ec_dec *_this){
154   int ret;
155   ret=ec_byte_read1(_this->buf);
156   if(ret<0){
157     ret=0;
158     /*Needed to keep oc_dec_tell() operating correctly.*/
159     ec_byte_adv1(_this->buf);
160   }
161   return ret;
162 }
163
164 /*Normalizes the contents of dif and rng so that rng lies entirely in the
165    high-order symbol.*/
166 static inline void ec_dec_normalize(ec_dec *_this){
167   /*If the range is too small, rescale it and input some bits.*/
168   while(_this->rng<=EC_CODE_BOT){
169     int sym;
170     _this->rng<<=EC_SYM_BITS;
171     /*Use up the remaining bits from our last symbol.*/
172     sym=_this->rem<<EC_CODE_EXTRA&EC_SYM_MAX;
173     /*Read the next value from the input.*/
174     _this->rem=ec_dec_in(_this);
175     /*Take the rest of the bits we need from this new symbol.*/
176     sym|=_this->rem>>EC_SYM_BITS-EC_CODE_EXTRA;
177     _this->dif=(_this->dif<<EC_SYM_BITS)-sym&EC_CODE_MASK;
178     /*dif can never be larger than EC_CODE_TOP.
179       This is equivalent to the slightly more readable:
180       if(_this->dif>EC_CODE_TOP)_this->dif-=EC_CODE_TOP;*/
181     _this->dif^=_this->dif&_this->dif-1&EC_CODE_TOP;
182   }
183 }
184
185 void ec_dec_init(ec_dec *_this,ec_byte_buffer *_buf){
186   _this->buf=_buf;
187   _this->rem=ec_dec_in(_this);
188   _this->rng=1U<<EC_CODE_EXTRA;
189   _this->dif=_this->rng-(_this->rem>>EC_SYM_BITS-EC_CODE_EXTRA);
190   /*Normalize the interval.*/
191   ec_dec_normalize(_this);
192 }
193
194
195 unsigned ec_decode(ec_dec *_this,unsigned _ft){
196   unsigned s;
197   _this->nrm=_this->rng/_ft;
198   s=(unsigned)((_this->dif-1)/_this->nrm);
199   return _ft-EC_MINI(s+1,_ft);
200 }
201
202 unsigned ec_decode_bin(ec_dec *_this,unsigned bits){
203    unsigned s;
204    ec_uint32 ft;
205    ft = (ec_uint32)1<<bits;
206    _this->nrm=_this->rng>>bits;
207    s=(unsigned)((_this->dif-1)/_this->nrm);
208    return ft-EC_MINI(s+1,ft);
209 }
210
211 void ec_dec_update(ec_dec *_this,unsigned _fl,unsigned _fh,unsigned _ft){
212   ec_uint32 s;
213   s=IMUL32(_this->nrm,(_ft-_fh));
214   _this->dif-=s;
215   _this->rng=_fl>0?IMUL32(_this->nrm,(_fh-_fl)):_this->rng-s;
216   ec_dec_normalize(_this);
217 }
218
219 long ec_dec_tell(ec_dec *_this,int _b){
220   ec_uint32 r;
221   int       l;
222   long      nbits;
223   nbits=(ec_byte_bytes(_this->buf)-(EC_CODE_BITS+EC_SYM_BITS-1)/EC_SYM_BITS)*
224    EC_SYM_BITS;
225   /*To handle the non-integral number of bits still left in the encoder state,
226      we compute the number of bits of low that must be encoded to ensure that
227      the value is inside the range for any possible subsequent bits.
228     Note that this is subtly different than the actual value we would end the
229      stream with, which tries to make as many of the trailing bits zeros as
230      possible.*/
231   nbits+=EC_CODE_BITS;
232   nbits<<=_b;
233   l=EC_ILOG(_this->rng);
234   r=_this->rng>>l-16;
235   while(_b-->0){
236     int b;
237     r=r*r>>15;
238     b=(int)(r>>16);
239     l=l<<1|b;
240     r>>=b;
241   }
242   return nbits-l;
243 }
244
245 #if 0
246 int ec_dec_done(ec_dec *_this){
247   unsigned low;
248   int      ret;
249   /*Check to make sure we've used all the input bytes.
250     This ensures that no more ones would ever be inserted into the decoder.*/
251   if(_this->buf->ptr-ec_byte_get_buffer(_this->buf)<=
252    ec_byte_bytes(_this->buf)){
253     return 0;
254   }
255   /*We compute the smallest finitely odd fraction that fits inside the current
256      range, and write that to the stream.
257     This is guaranteed to yield the smallest possible encoding.*/
258   /*TODO: Fix this line, as it is wrong.
259     It doesn't seem worth being able to make this check to do an extra
260      subtraction for every symbol decoded.*/
261   low=/*What we want: _this->top-_this->rng; What we have:*/_this->dif
262   if(low){
263     unsigned end;
264     end=EC_CODE_TOP;
265     /*Ensure that the next free end is in the range.*/
266     if(end-low>=_this->rng){
267       unsigned msk;
268       msk=EC_CODE_TOP-1;
269       do{
270         msk>>=1;
271         end=(low+msk)&~msk|msk+1;
272       }
273       while(end-low>=_this->rng);
274     }
275     /*The remaining input should have been the next free end.*/
276     return end-low!=_this->dif;
277   }
278   return 1;
279 }
280 #endif