making codec draft more compliant with IETF submission rules
[opus.git] / libcelt / mfrngdec.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 multiply-free 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   The coder is made multiply-free by replacing the standard multiply/divide
65    used to partition the current interval according to the total frequency
66    count.
67   The new partition function scales the count so that it differs from the size
68    of the interval by no more than a factor of two and then assigns each symbol
69    one or two code words in the interval.
70   For details see \cite{SM98}.
71
72   This coder also handles the end of the stream in a slightly more graceful
73    fashion than most arithmetic or range coders.
74   Once the final symbol has been encoded, the coder selects the code word with
75    the shortest number of bits that still falls within the final interval.
76   This method is not novel.
77   Here, by the length of the code word, we refer to the number of bits until
78    its final 1.
79   Any trailing zeros may be discarded, since the encoder, once it runs out of
80    input, will pad its buffer with zeros.
81
82   But this means that no encoded stream would ever have any zero bytes at the
83    end.
84   Since there are some coded representations we cannot produce, it implies that
85    there is still some redundancy in the stream.
86   In this case, we can pick a special byte value, RSV1, and should the stream
87    end in a sequence of zeros, followed by the RSV1 byte, we can code the
88    zeros, and discard the RSV1 byte.
89   The decoder, knowing that the encoder would never produce a sequence of zeros
90    at the end, would then know to add in the RSV1 byte if it observed it.
91
92   Now, the encoder would never produce a stream that ended in a sequence of
93    zeros followed by a RSV1 byte.
94   So, if the stream ends in a non-empty sequence of zeros, followed by any
95    positive number of RSV1 bytes, the last RSV1 byte is discarded.
96   The decoder, if it encounters a stream that ends in non-empty sequence of
97    zeros followed by any non-negative number of RSV1 bytes, adds an additional
98    RSV1 byte to the stream.
99   With this strategy, every possible sequence of input bytes is transformed to
100    one that could actually be produced by the encoder.
101
102   The only question is what non-zero value to use for RSV1.
103   We select 0x80, since it has the nice property of producing the shortest
104    possible byte streams when using our strategy for selecting a number within
105    the final interval to encode.
106   Clearly if the shortest possible code word that falls within the interval has
107    its last one bit as the most significant bit of the final byte, and the
108    previous bytes were a non-empty sequence of zeros followed by a non-negative
109    number of 0x80 bytes, then the last byte would be discarded.
110   If the shortest code word is not so formed, then no other code word in the
111    interval would result in any more bytes being discarded.
112   Any longer code word would have an additional one bit somewhere, and so would
113    require at a minimum that that byte would be coded.
114   If the shortest code word has a 1 before the final one that is preventing the
115    stream from ending in a non-empty sequence of zeros followed by a
116    non-negative number of 0x80's, then there is no code word of the same length
117    which contains that bit as a zero.
118   If there were, then we could simply leave that bit a 1, and drop all the bits
119    after it without leaving the interval, thus producing a shorter code word.
120
121   In this case, RSV1 can only drop 1 bit off the final stream.
122   Other choices could lead to savings of up to 8 bits for particular streams,
123    but this would produce the odd situation that a stream with more non-zero
124    bits is actually encoded in fewer bytes.
125
126   @PHDTHESIS{Pas76,
127     author="Richard Clark Pasco",
128     title="Source coding algorithms for fast data compression",
129     school="Dept. of Electrical Engineering, Stanford University",
130     address="Stanford, CA",
131     month=May,
132     year=1976
133   }
134   @INPROCEEDINGS{Mar79,
135    author="Martin, G.N.N.",
136    title="Range encoding: an algorithm for removing redundancy from a digitised
137     message",
138    booktitle="Video & Data Recording Conference",
139    year=1979,
140    address="Southampton",
141    month=Jul
142   }
143   @ARTICLE{MNW98,
144    author="Alistair Moffat and Radford Neal and Ian H. Witten",
145    title="Arithmetic Coding Revisited",
146    journal="{ACM} Transactions on Information Systems",
147    year=1998,
148    volume=16,
149    number=3,
150    pages="256--294",
151    month=Jul,
152    URL="http://www.stanford.edu/class/ee398/handouts/papers/Moffat98ArithmCoding.pdf"
153   }
154   @INPROCEEDINGS{SM98,
155    author="Lang Stuiver and Alistair Moffat",
156    title="Piecewise Integer Mapping for Arithmetic Coding",
157    booktitle="Proceedings of the {IEEE} Data Compression Conference",
158    pages="1--10",
159    address="Snowbird, UT",
160    month="Mar./Apr.",
161    year=1998
162   }*/
163
164
165 /*Gets the next byte of input.
166   After all the bytes in the current packet have been consumed, and the extra
167    end code returned if needed, this function will continue to return zero each
168    time it is called.
169   Return: The next byte of input.*/
170 static int ec_dec_in(ec_dec *_this){
171   int ret;
172   ret=ec_byte_read1(_this->buf);
173   if(ret<0){
174     ret=0;
175     /*Needed to keep oc_dec_tell() operating correctly.*/
176     ec_byte_adv1(_this->buf);
177   }
178   return ret;
179 }
180
181 /*Normalizes the contents of dif and rng so that rng lies entirely in the
182    high-order symbol.*/
183 static inline void ec_dec_normalize(ec_dec *_this){
184   /*If the range is too small, rescale it and input some bits.*/
185   while(_this->rng<=EC_CODE_BOT){
186     int sym;
187     _this->rng<<=EC_SYM_BITS;
188     /*Use up the remaining bits from our last symbol.*/
189     sym=_this->rem<<EC_CODE_EXTRA&EC_SYM_MAX;
190     /*Read the next value from the input.*/
191     _this->rem=ec_dec_in(_this);
192     /*Take the rest of the bits we need from this new symbol.*/
193     sym|=_this->rem>>EC_SYM_BITS-EC_CODE_EXTRA;
194     _this->dif=(_this->dif<<EC_SYM_BITS)+sym&EC_CODE_MASK;
195     /*dif can never be larger than EC_CODE_TOP.
196       This is equivalent to the slightly more readable:
197       if(_this->dif>EC_CODE_TOP)_this->dif-=EC_CODE_TOP;*/
198     _this->dif^=_this->dif&_this->dif-1&EC_CODE_TOP;
199   }
200 }
201
202 void ec_dec_init(ec_dec *_this,ec_byte_buffer *_buf){
203   _this->buf=_buf;
204   _this->rem=ec_dec_in(_this);
205   _this->rng=1U<<EC_CODE_EXTRA;
206   _this->dif=_this->rem>>EC_SYM_BITS-EC_CODE_EXTRA;
207   /*Normalize the interval.*/
208   ec_dec_normalize(_this);
209 }
210
211 unsigned ec_decode(ec_dec *_this,unsigned _ft){
212   ec_uint32 ft;
213   ec_uint32 d;
214   unsigned  e;
215   /*Step 1: Compute the normalization factor for the frequency counts.*/
216   _this->nrm=EC_ILOG(_this->rng)-EC_ILOG(_ft);
217   ft=(ec_uint32)_ft<<_this->nrm;
218   e=ft>_this->rng;
219   ft>>=e;
220   _this->nrm-=e;
221   /*Step 2: invert the partition function.*/
222   d=_this->rng-ft;
223   return EC_MAXI((ec_int32)(_this->dif>>1),(ec_int32)(_this->dif-d))>>
224    _this->nrm;
225   /*Step 3: The caller locates the range [fl,fh) containing the return value
226      and calls ec_dec_update().*/
227 }
228
229 unsigned ec_decode_bin(ec_dec *_this,unsigned bits){
230   return ec_decode(_this, 1U<<bits);
231 }
232
233 void ec_dec_update(ec_dec *_this,unsigned _fl,unsigned _fh,unsigned _ft){
234   ec_uint32 fl;
235   ec_uint32 fh;
236   ec_uint32 ft;
237   ec_uint32 r;
238   ec_uint32 s;
239   ec_uint32 d;
240   /*Step 4: Evaluate the two partition function values.*/
241   fl=(ec_uint32)_fl<<_this->nrm;
242   fh=(ec_uint32)_fh<<_this->nrm;
243   ft=(ec_uint32)_ft<<_this->nrm;
244   d=_this->rng-ft;
245   r=fh+EC_MINI(fh,d);
246   s=fl+EC_MINI(fl,d);
247   /*Step 5: Update the interval.*/
248   _this->rng=r-s;
249   _this->dif-=s;
250   /*Step 6: Normalize the interval.*/
251   ec_dec_normalize(_this);
252 }
253
254 long ec_dec_tell(ec_dec *_this,int _b){
255   ec_uint32 r;
256   int       l;
257   long      nbits;
258   nbits=(ec_byte_bytes(_this->buf)-(EC_CODE_BITS+EC_SYM_BITS-1)/EC_SYM_BITS)*
259    EC_SYM_BITS;
260   /*To handle the non-integral number of bits still left in the encoder state,
261      we compute the number of bits of low that must be encoded to ensure that
262      the value is inside the range for any possible subsequent bits.
263     Note that this is subtly different than the actual value we would end the
264      stream with, which tries to make as many of the trailing bits zeros as
265      possible.*/
266   nbits+=EC_CODE_BITS;
267   nbits<<=_b;
268   l=EC_ILOG(_this->rng);
269   r=_this->rng>>l-16;
270   while(_b-->0){
271     int b;
272     r=r*r>>15;
273     b=(int)(r>>16);
274     l=l<<1|b;
275     r>>=b;
276   }
277   return nbits-l;
278 }
279
280 #if 0
281 int ec_dec_done(ec_dec *_this){
282   unsigned low;
283   int      ret;
284   /*Check to make sure we've used all the input bytes.
285     This ensures that no more ones would ever be inserted into the decoder.*/
286   if(_this->buf->ptr-ec_byte_get_buffer(_this->buf)<=
287    ec_byte_bytes(_this->buf)){
288     return 0;
289   }
290   /*We compute the smallest finitely odd fraction that fits inside the current
291      range, and write that to the stream.
292     This is guaranteed to yield the smallest possible encoding.*/
293   /*TODO: Fix this line, as it is wrong.
294     It doesn't seem worth being able to make this check to do an extra
295      subtraction for every symbol decoded.*/
296   low=/*What we want: _this->top-_this->rng; What we have:*/_this->dif
297   if(low){
298     unsigned end;
299     end=EC_CODE_TOP;
300     /*Ensure that the next free end is in the range.*/
301     if(end-low>=_this->rng){
302       unsigned msk;
303       msk=EC_CODE_TOP-1;
304       do{
305         msk>>=1;
306         end=low+msk&~msk|msk+1;
307       }
308       while(end-low>=_this->rng);
309     }
310     /*The remaining input should have been the next free end.*/
311     return end-low!=_this->dif;
312   }
313   return 1;
314 }
315 #endif