Adds a fuzzing mode that causes the encoder to make random decisions
[opus.git] / libcelt / cwrs.c
1 /* Copyright (c) 2007-2008 CSIRO
2    Copyright (c) 2007-2009 Xiph.Org Foundation
3    Copyright (c) 2007-2009 Timothy B. Terriberry
4    Written by Timothy B. Terriberry and Jean-Marc Valin */
5 /*
6    Redistribution and use in source and binary forms, with or without
7    modification, are permitted provided that the following conditions
8    are met:
9
10    - Redistributions of source code must retain the above copyright
11    notice, this list of conditions and the following disclaimer.
12
13    - Redistributions in binary form must reproduce the above copyright
14    notice, this list of conditions and the following disclaimer in the
15    documentation and/or other materials provided with the distribution.
16
17    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
18    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
19    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
20    A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
21    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
22    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
23    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
24    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
25    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
26    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
27    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28 */
29
30 #ifdef HAVE_CONFIG_H
31 #include "config.h"
32 #endif
33
34 #include "os_support.h"
35 #include <stdlib.h>
36 #include <string.h>
37 #include "cwrs.h"
38 #include "mathops.h"
39 #include "arch.h"
40
41 #ifdef CUSTOM_MODES
42
43 /*Guaranteed to return a conservatively large estimate of the binary logarithm
44    with frac bits of fractional precision.
45   Tested for all possible 32-bit inputs with frac=4, where the maximum
46    overestimation is 0.06254243 bits.*/
47 int log2_frac(opus_uint32 val, int frac)
48 {
49   int l;
50   l=EC_ILOG(val);
51   if(val&val-1){
52     /*This is (val>>l-16), but guaranteed to round up, even if adding a bias
53        before the shift would cause overflow (e.g., for 0xFFFFxxxx).*/
54     if(l>16)val=(val>>l-16)+((val&(1<<l-16)-1)+(1<<l-16)-1>>l-16);
55     else val<<=16-l;
56     l=l-1<<frac;
57     /*Note that we always need one iteration, since the rounding up above means
58        that we might need to adjust the integer part of the logarithm.*/
59     do{
60       int b;
61       b=(int)(val>>16);
62       l+=b<<frac;
63       val=val+b>>b;
64       val=val*val+0x7FFF>>15;
65     }
66     while(frac-->0);
67     /*If val is not exactly 0x8000, then we have to round up the remainder.*/
68     return l+(val>0x8000);
69   }
70   /*Exact powers of two require no rounding.*/
71   else return l-1<<frac;
72 }
73 #endif
74
75 #ifndef SMALL_FOOTPRINT
76
77 #define MASK32 (0xFFFFFFFF)
78
79 /*INV_TABLE[i] holds the multiplicative inverse of (2*i+1) mod 2**32.*/
80 static const opus_uint32 INV_TABLE[64]={
81   0x00000001,0xAAAAAAAB,0xCCCCCCCD,0xB6DB6DB7,
82   0x38E38E39,0xBA2E8BA3,0xC4EC4EC5,0xEEEEEEEF,
83   0xF0F0F0F1,0x286BCA1B,0x3CF3CF3D,0xE9BD37A7,
84   0xC28F5C29,0x684BDA13,0x4F72C235,0xBDEF7BDF,
85   0x3E0F83E1,0x8AF8AF8B,0x914C1BAD,0x96F96F97,
86   0xC18F9C19,0x2FA0BE83,0xA4FA4FA5,0x677D46CF,
87   0x1A1F58D1,0xFAFAFAFB,0x8C13521D,0x586FB587,
88   0xB823EE09,0xA08AD8F3,0xC10C9715,0xBEFBEFBF,
89   0xC0FC0FC1,0x07A44C6B,0xA33F128D,0xE327A977,
90   0xC7E3F1F9,0x962FC963,0x3F2B3885,0x613716AF,
91   0x781948B1,0x2B2E43DB,0xFCFCFCFD,0x6FD0EB67,
92   0xFA3F47E9,0xD2FD2FD3,0x3F4FD3F5,0xD4E25B9F,
93   0x5F02A3A1,0xBF5A814B,0x7C32B16D,0xD3431B57,
94   0xD8FD8FD9,0x8D28AC43,0xDA6C0965,0xDB195E8F,
95   0x0FDBC091,0x61F2A4BB,0xDCFDCFDD,0x46FDD947,
96   0x56BE69C9,0xEB2FDEB3,0x26E978D5,0xEFDFBF7F,
97   /*
98   0x0FE03F81,0xC9484E2B,0xE133F84D,0xE1A8C537,
99   0x077975B9,0x70586723,0xCD29C245,0xFAA11E6F,
100   0x0FE3C071,0x08B51D9B,0x8CE2CABD,0xBF937F27,
101   0xA8FE53A9,0x592FE593,0x2C0685B5,0x2EB11B5F,
102   0xFCD1E361,0x451AB30B,0x72CFE72D,0xDB35A717,
103   0xFB74A399,0xE80BFA03,0x0D516325,0x1BCB564F,
104   0xE02E4851,0xD962AE7B,0x10F8ED9D,0x95AEDD07,
105   0xE9DC0589,0xA18A4473,0xEA53FA95,0xEE936F3F,
106   0x90948F41,0xEAFEAFEB,0x3D137E0D,0xEF46C0F7,
107   0x028C1979,0x791064E3,0xC04FEC05,0xE115062F,
108   0x32385831,0x6E68575B,0xA10D387D,0x6FECF2E7,
109   0x3FB47F69,0xED4BFB53,0x74FED775,0xDB43BB1F,
110   0x87654321,0x9BA144CB,0x478BBCED,0xBFB912D7,
111   0x1FDCD759,0x14B2A7C3,0xCB125CE5,0x437B2E0F,
112   0x10FEF011,0xD2B3183B,0x386CAB5D,0xEF6AC0C7,
113   0x0E64C149,0x9A020A33,0xE6B41C55,0xFEFEFEFF*/
114 };
115
116 /*Computes (_a*_b-_c)/(2*_d+1) when the quotient is known to be exact.
117   _a, _b, _c, and _d may be arbitrary so long as the arbitrary precision result
118    fits in 32 bits, but currently the table for multiplicative inverses is only
119    valid for _d<128.*/
120 static inline opus_uint32 imusdiv32odd(opus_uint32 _a,opus_uint32 _b,
121  opus_uint32 _c,int _d){
122   return (_a*_b-_c)*INV_TABLE[_d]&MASK32;
123 }
124
125 /*Computes (_a*_b-_c)/_d when the quotient is known to be exact.
126   _d does not actually have to be even, but imusdiv32odd will be faster when
127    it's odd, so you should use that instead.
128   _a and _d are assumed to be small (e.g., _a*_d fits in 32 bits; currently the
129    table for multiplicative inverses is only valid for _d<=256).
130   _b and _c may be arbitrary so long as the arbitrary precision reuslt fits in
131    32 bits.*/
132 static inline opus_uint32 imusdiv32even(opus_uint32 _a,opus_uint32 _b,
133  opus_uint32 _c,int _d){
134   opus_uint32 inv;
135   int           mask;
136   int           shift;
137   int           one;
138   celt_assert(_d>0);
139   shift=EC_ILOG(_d^_d-1);
140   celt_assert(_d<=256);
141   inv=INV_TABLE[_d-1>>shift];
142   shift--;
143   one=1<<shift;
144   mask=one-1;
145   return (_a*(_b>>shift)-(_c>>shift)+
146    (_a*(_b&mask)+one-(_c&mask)>>shift)-1)*inv&MASK32;
147 }
148
149 #endif /* SMALL_FOOTPRINT */
150
151 /*Although derived separately, the pulse vector coding scheme is equivalent to
152    a Pyramid Vector Quantizer \cite{Fis86}.
153   Some additional notes about an early version appear at
154    http://people.xiph.org/~tterribe/notes/cwrs.html, but the codebook ordering
155    and the definitions of some terms have evolved since that was written.
156
157   The conversion from a pulse vector to an integer index (encoding) and back
158    (decoding) is governed by two related functions, V(N,K) and U(N,K).
159
160   V(N,K) = the number of combinations, with replacement, of N items, taken K
161    at a time, when a sign bit is added to each item taken at least once (i.e.,
162    the number of N-dimensional unit pulse vectors with K pulses).
163   One way to compute this is via
164     V(N,K) = K>0 ? sum(k=1...K,2**k*choose(N,k)*choose(K-1,k-1)) : 1,
165    where choose() is the binomial function.
166   A table of values for N<10 and K<10 looks like:
167   V[10][10] = {
168     {1,  0,   0,    0,    0,     0,     0,      0,      0,       0},
169     {1,  2,   2,    2,    2,     2,     2,      2,      2,       2},
170     {1,  4,   8,   12,   16,    20,    24,     28,     32,      36},
171     {1,  6,  18,   38,   66,   102,   146,    198,    258,     326},
172     {1,  8,  32,   88,  192,   360,   608,    952,   1408,    1992},
173     {1, 10,  50,  170,  450,  1002,  1970,   3530,   5890,    9290},
174     {1, 12,  72,  292,  912,  2364,  5336,  10836,  20256,   35436},
175     {1, 14,  98,  462, 1666,  4942, 12642,  28814,  59906,  115598},
176     {1, 16, 128,  688, 2816,  9424, 27008,  68464, 157184,  332688},
177     {1, 18, 162,  978, 4482, 16722, 53154, 148626, 374274,  864146}
178   };
179
180   U(N,K) = the number of such combinations wherein N-1 objects are taken at
181    most K-1 at a time.
182   This is given by
183     U(N,K) = sum(k=0...K-1,V(N-1,k))
184            = K>0 ? (V(N-1,K-1) + V(N,K-1))/2 : 0.
185   The latter expression also makes clear that U(N,K) is half the number of such
186    combinations wherein the first object is taken at least once.
187   Although it may not be clear from either of these definitions, U(N,K) is the
188    natural function to work with when enumerating the pulse vector codebooks,
189    not V(N,K).
190   U(N,K) is not well-defined for N=0, but with the extension
191     U(0,K) = K>0 ? 0 : 1,
192    the function becomes symmetric: U(N,K) = U(K,N), with a similar table:
193   U[10][10] = {
194     {1, 0,  0,   0,    0,    0,     0,     0,      0,      0},
195     {0, 1,  1,   1,    1,    1,     1,     1,      1,      1},
196     {0, 1,  3,   5,    7,    9,    11,    13,     15,     17},
197     {0, 1,  5,  13,   25,   41,    61,    85,    113,    145},
198     {0, 1,  7,  25,   63,  129,   231,   377,    575,    833},
199     {0, 1,  9,  41,  129,  321,   681,  1289,   2241,   3649},
200     {0, 1, 11,  61,  231,  681,  1683,  3653,   7183,  13073},
201     {0, 1, 13,  85,  377, 1289,  3653,  8989,  19825,  40081},
202     {0, 1, 15, 113,  575, 2241,  7183, 19825,  48639, 108545},
203     {0, 1, 17, 145,  833, 3649, 13073, 40081, 108545, 265729}
204   };
205
206   With this extension, V(N,K) may be written in terms of U(N,K):
207     V(N,K) = U(N,K) + U(N,K+1)
208    for all N>=0, K>=0.
209   Thus U(N,K+1) represents the number of combinations where the first element
210    is positive or zero, and U(N,K) represents the number of combinations where
211    it is negative.
212   With a large enough table of U(N,K) values, we could write O(N) encoding
213    and O(min(N*log(K),N+K)) decoding routines, but such a table would be
214    prohibitively large for small embedded devices (K may be as large as 32767
215    for small N, and N may be as large as 200).
216
217   Both functions obey the same recurrence relation:
218     V(N,K) = V(N-1,K) + V(N,K-1) + V(N-1,K-1),
219     U(N,K) = U(N-1,K) + U(N,K-1) + U(N-1,K-1),
220    for all N>0, K>0, with different initial conditions at N=0 or K=0.
221   This allows us to construct a row of one of the tables above given the
222    previous row or the next row.
223   Thus we can derive O(NK) encoding and decoding routines with O(K) memory
224    using only addition and subtraction.
225
226   When encoding, we build up from the U(2,K) row and work our way forwards.
227   When decoding, we need to start at the U(N,K) row and work our way backwards,
228    which requires a means of computing U(N,K).
229   U(N,K) may be computed from two previous values with the same N:
230     U(N,K) = ((2*N-1)*U(N,K-1) - U(N,K-2))/(K-1) + U(N,K-2)
231    for all N>1, and since U(N,K) is symmetric, a similar relation holds for two
232    previous values with the same K:
233     U(N,K>1) = ((2*K-1)*U(N-1,K) - U(N-2,K))/(N-1) + U(N-2,K)
234    for all K>1.
235   This allows us to construct an arbitrary row of the U(N,K) table by starting
236    with the first two values, which are constants.
237   This saves roughly 2/3 the work in our O(NK) decoding routine, but costs O(K)
238    multiplications.
239   Similar relations can be derived for V(N,K), but are not used here.
240
241   For N>0 and K>0, U(N,K) and V(N,K) take on the form of an (N-1)-degree
242    polynomial for fixed N.
243   The first few are
244     U(1,K) = 1,
245     U(2,K) = 2*K-1,
246     U(3,K) = (2*K-2)*K+1,
247     U(4,K) = (((4*K-6)*K+8)*K-3)/3,
248     U(5,K) = ((((2*K-4)*K+10)*K-8)*K+3)/3,
249    and
250     V(1,K) = 2,
251     V(2,K) = 4*K,
252     V(3,K) = 4*K*K+2,
253     V(4,K) = 8*(K*K+2)*K/3,
254     V(5,K) = ((4*K*K+20)*K*K+6)/3,
255    for all K>0.
256   This allows us to derive O(N) encoding and O(N*log(K)) decoding routines for
257    small N (and indeed decoding is also O(N) for N<3).
258
259   @ARTICLE{Fis86,
260     author="Thomas R. Fischer",
261     title="A Pyramid Vector Quantizer",
262     journal="IEEE Transactions on Information Theory",
263     volume="IT-32",
264     number=4,
265     pages="568--583",
266     month=Jul,
267     year=1986
268   }*/
269
270 #ifndef SMALL_FOOTPRINT
271
272 /*Compute V(1,_k).*/
273 static inline unsigned ncwrs1(int _k){
274   return _k?2:1;
275 }
276
277 /*Compute U(2,_k).
278   Note that this may be called with _k=32768 (maxK[2]+1).*/
279 static inline unsigned ucwrs2(unsigned _k){
280   return _k?_k+(_k-1):0;
281 }
282
283 /*Compute V(2,_k).*/
284 static inline opus_uint32 ncwrs2(int _k){
285   return _k?4*(opus_uint32)_k:1;
286 }
287
288 /*Compute U(3,_k).
289   Note that this may be called with _k=32768 (maxK[3]+1).*/
290 static inline opus_uint32 ucwrs3(unsigned _k){
291   return _k?(2*(opus_uint32)_k-2)*_k+1:0;
292 }
293
294 /*Compute V(3,_k).*/
295 static inline opus_uint32 ncwrs3(int _k){
296   return _k?2*(2*(unsigned)_k*(opus_uint32)_k+1):1;
297 }
298
299 /*Compute U(4,_k).*/
300 static inline opus_uint32 ucwrs4(int _k){
301   return _k?imusdiv32odd(2*_k,(2*_k-3)*(opus_uint32)_k+4,3,1):0;
302 }
303
304 /*Compute V(4,_k).*/
305 static inline opus_uint32 ncwrs4(int _k){
306   return _k?((_k*(opus_uint32)_k+2)*_k)/3<<3:1;
307 }
308
309 /*Compute U(5,_k).*/
310 static inline opus_uint32 ucwrs5(int _k){
311   return _k?(((((_k-2)*(unsigned)_k+5)*(opus_uint32)_k-4)*_k)/3<<1)+1:0;
312 }
313
314 /*Compute V(5,_k).*/
315 static inline opus_uint32 ncwrs5(int _k){
316   return _k?(((_k*(unsigned)_k+5)*(opus_uint32)_k*_k)/3<<2)+2:1;
317 }
318
319 #endif /* SMALL_FOOTPRINT */
320
321 /*Computes the next row/column of any recurrence that obeys the relation
322    u[i][j]=u[i-1][j]+u[i][j-1]+u[i-1][j-1].
323   _ui0 is the base case for the new row/column.*/
324 static inline void unext(opus_uint32 *_ui,unsigned _len,opus_uint32 _ui0){
325   opus_uint32 ui1;
326   unsigned      j;
327   /*This do-while will overrun the array if we don't have storage for at least
328      2 values.*/
329   j=1; do {
330     ui1=UADD32(UADD32(_ui[j],_ui[j-1]),_ui0);
331     _ui[j-1]=_ui0;
332     _ui0=ui1;
333   } while (++j<_len);
334   _ui[j-1]=_ui0;
335 }
336
337 /*Computes the previous row/column of any recurrence that obeys the relation
338    u[i-1][j]=u[i][j]-u[i][j-1]-u[i-1][j-1].
339   _ui0 is the base case for the new row/column.*/
340 static inline void uprev(opus_uint32 *_ui,unsigned _n,opus_uint32 _ui0){
341   opus_uint32 ui1;
342   unsigned      j;
343   /*This do-while will overrun the array if we don't have storage for at least
344      2 values.*/
345   j=1; do {
346     ui1=USUB32(USUB32(_ui[j],_ui[j-1]),_ui0);
347     _ui[j-1]=_ui0;
348     _ui0=ui1;
349   } while (++j<_n);
350   _ui[j-1]=_ui0;
351 }
352
353 /*Compute V(_n,_k), as well as U(_n,0..._k+1).
354   _u: On exit, _u[i] contains U(_n,i) for i in [0..._k+1].*/
355 static opus_uint32 ncwrs_urow(unsigned _n,unsigned _k,opus_uint32 *_u){
356   opus_uint32 um2;
357   unsigned      len;
358   unsigned      k;
359   len=_k+2;
360   /*We require storage at least 3 values (e.g., _k>0).*/
361   celt_assert(len>=3);
362   _u[0]=0;
363   _u[1]=um2=1;
364 #ifndef SMALL_FOOTPRINT
365   if(_n<=6 || _k>255)
366 #endif
367  {
368     /*If _n==0, _u[0] should be 1 and the rest should be 0.*/
369     /*If _n==1, _u[i] should be 1 for i>1.*/
370     celt_assert(_n>=2);
371     /*If _k==0, the following do-while loop will overflow the buffer.*/
372     celt_assert(_k>0);
373     k=2;
374     do _u[k]=(k<<1)-1;
375     while(++k<len);
376     for(k=2;k<_n;k++)unext(_u+1,_k+1,1);
377   }
378 #ifndef SMALL_FOOTPRINT
379   else{
380     opus_uint32 um1;
381     opus_uint32 n2m1;
382     _u[2]=n2m1=um1=(_n<<1)-1;
383     for(k=3;k<len;k++){
384       /*U(N,K) = ((2*N-1)*U(N,K-1)-U(N,K-2))/(K-1) + U(N,K-2)*/
385       _u[k]=um2=imusdiv32even(n2m1,um1,um2,k-1)+um2;
386       if(++k>=len)break;
387       _u[k]=um1=imusdiv32odd(n2m1,um2,um1,k-1>>1)+um1;
388     }
389   }
390 #endif /* SMALL_FOOTPRINT */
391   return _u[_k]+_u[_k+1];
392 }
393
394 /*Returns the _i'th combination of _k elements (at most 32767) chosen from a
395    set of size 1 with associated sign bits.
396   _y: Returns the vector of pulses.*/
397 static inline void cwrsi1(int _k,opus_uint32 _i,int *_y){
398   int s;
399   s=-(int)_i;
400   _y[0]=_k+s^s;
401 }
402
403 #ifndef SMALL_FOOTPRINT
404
405 /*Returns the _i'th combination of _k elements (at most 32767) chosen from a
406    set of size 2 with associated sign bits.
407   _y: Returns the vector of pulses.*/
408 static inline void cwrsi2(int _k,opus_uint32 _i,int *_y){
409   opus_uint32 p;
410   int           s;
411   int           yj;
412   p=ucwrs2(_k+1U);
413   s=-(_i>=p);
414   _i-=p&s;
415   yj=_k;
416   _k=_i+1>>1;
417   p=ucwrs2(_k);
418   _i-=p;
419   yj-=_k;
420   _y[0]=yj+s^s;
421   cwrsi1(_k,_i,_y+1);
422 }
423
424 /*Returns the _i'th combination of _k elements (at most 32767) chosen from a
425    set of size 3 with associated sign bits.
426   _y: Returns the vector of pulses.*/
427 static void cwrsi3(int _k,opus_uint32 _i,int *_y){
428   opus_uint32 p;
429   int           s;
430   int           yj;
431   p=ucwrs3(_k+1U);
432   s=-(_i>=p);
433   _i-=p&s;
434   yj=_k;
435   /*Finds the maximum _k such that ucwrs3(_k)<=_i (tested for all
436      _i<2147418113=U(3,32768)).*/
437   _k=_i>0?isqrt32(2*_i-1)+1>>1:0;
438   p=ucwrs3(_k);
439   _i-=p;
440   yj-=_k;
441   _y[0]=yj+s^s;
442   cwrsi2(_k,_i,_y+1);
443 }
444
445 /*Returns the _i'th combination of _k elements (at most 1172) chosen from a set
446    of size 4 with associated sign bits.
447   _y: Returns the vector of pulses.*/
448 static void cwrsi4(int _k,opus_uint32 _i,int *_y){
449   opus_uint32 p;
450   int           s;
451   int           yj;
452   int           kl;
453   int           kr;
454   p=ucwrs4(_k+1);
455   s=-(_i>=p);
456   _i-=p&s;
457   yj=_k;
458   /*We could solve a cubic for k here, but the form of the direct solution does
459      not lend itself well to exact integer arithmetic.
460     Instead we do a binary search on U(4,K).*/
461   kl=0;
462   kr=_k;
463   for(;;){
464     _k=kl+kr>>1;
465     p=ucwrs4(_k);
466     if(p<_i){
467       if(_k>=kr)break;
468       kl=_k+1;
469     }
470     else if(p>_i)kr=_k-1;
471     else break;
472   }
473   _i-=p;
474   yj-=_k;
475   _y[0]=yj+s^s;
476   cwrsi3(_k,_i,_y+1);
477 }
478
479 /*Returns the _i'th combination of _k elements (at most 238) chosen from a set
480    of size 5 with associated sign bits.
481   _y: Returns the vector of pulses.*/
482 static void cwrsi5(int _k,opus_uint32 _i,int *_y){
483   opus_uint32 p;
484   int           s;
485   int           yj;
486   p=ucwrs5(_k+1);
487   s=-(_i>=p);
488   _i-=p&s;
489   yj=_k;
490   /* A binary search on U(5,K) avoids the need for 64-bit arithmetic */
491   {
492     int kl=0;
493     int kr=_k;
494     for(;;){
495       _k=kl+kr>>1;
496       p=ucwrs5(_k);
497       if(p<_i){
498         if(_k>=kr)break;
499         kl=_k+1;
500       }
501       else if(p>_i)kr=_k-1;
502       else break;
503     }
504   }
505   _i-=p;
506   yj-=_k;
507   _y[0]=yj+s^s;
508   cwrsi4(_k,_i,_y+1);
509 }
510 #endif /* SMALL_FOOTPRINT */
511
512 /*Returns the _i'th combination of _k elements chosen from a set of size _n
513    with associated sign bits.
514   _y: Returns the vector of pulses.
515   _u: Must contain entries [0..._k+1] of row _n of U() on input.
516       Its contents will be destructively modified.*/
517 static void cwrsi(int _n,int _k,opus_uint32 _i,int *_y,opus_uint32 *_u){
518   int j;
519   celt_assert(_n>0);
520   j=0;
521   do{
522     opus_uint32 p;
523     int           s;
524     int           yj;
525     p=_u[_k+1];
526     s=-(_i>=p);
527     _i-=p&s;
528     yj=_k;
529     p=_u[_k];
530     while(p>_i)p=_u[--_k];
531     _i-=p;
532     yj-=_k;
533     _y[j]=yj+s^s;
534     uprev(_u,_k+2,0);
535   }
536   while(++j<_n);
537 }
538
539 /*Returns the index of the given combination of K elements chosen from a set
540    of size 1 with associated sign bits.
541   _y: The vector of pulses, whose sum of absolute values is K.
542   _k: Returns K.*/
543 static inline opus_uint32 icwrs1(const int *_y,int *_k){
544   *_k=abs(_y[0]);
545   return _y[0]<0;
546 }
547
548 #ifndef SMALL_FOOTPRINT
549
550 /*Returns the index of the given combination of K elements chosen from a set
551    of size 2 with associated sign bits.
552   _y: The vector of pulses, whose sum of absolute values is K.
553   _k: Returns K.*/
554 static inline opus_uint32 icwrs2(const int *_y,int *_k){
555   opus_uint32 i;
556   int           k;
557   i=icwrs1(_y+1,&k);
558   i+=ucwrs2(k);
559   k+=abs(_y[0]);
560   if(_y[0]<0)i+=ucwrs2(k+1U);
561   *_k=k;
562   return i;
563 }
564
565 /*Returns the index of the given combination of K elements chosen from a set
566    of size 3 with associated sign bits.
567   _y: The vector of pulses, whose sum of absolute values is K.
568   _k: Returns K.*/
569 static inline opus_uint32 icwrs3(const int *_y,int *_k){
570   opus_uint32 i;
571   int           k;
572   i=icwrs2(_y+1,&k);
573   i+=ucwrs3(k);
574   k+=abs(_y[0]);
575   if(_y[0]<0)i+=ucwrs3(k+1U);
576   *_k=k;
577   return i;
578 }
579
580 /*Returns the index of the given combination of K elements chosen from a set
581    of size 4 with associated sign bits.
582   _y: The vector of pulses, whose sum of absolute values is K.
583   _k: Returns K.*/
584 static inline opus_uint32 icwrs4(const int *_y,int *_k){
585   opus_uint32 i;
586   int           k;
587   i=icwrs3(_y+1,&k);
588   i+=ucwrs4(k);
589   k+=abs(_y[0]);
590   if(_y[0]<0)i+=ucwrs4(k+1);
591   *_k=k;
592   return i;
593 }
594
595 /*Returns the index of the given combination of K elements chosen from a set
596    of size 5 with associated sign bits.
597   _y: The vector of pulses, whose sum of absolute values is K.
598   _k: Returns K.*/
599 static inline opus_uint32 icwrs5(const int *_y,int *_k){
600   opus_uint32 i;
601   int           k;
602   i=icwrs4(_y+1,&k);
603   i+=ucwrs5(k);
604   k+=abs(_y[0]);
605   if(_y[0]<0)i+=ucwrs5(k+1);
606   *_k=k;
607   return i;
608 }
609 #endif /* SMALL_FOOTPRINT */
610
611 /*Returns the index of the given combination of K elements chosen from a set
612    of size _n with associated sign bits.
613   _y:  The vector of pulses, whose sum of absolute values must be _k.
614   _nc: Returns V(_n,_k).*/
615 opus_uint32 icwrs(int _n,int _k,opus_uint32 *_nc,const int *_y,
616  opus_uint32 *_u){
617   opus_uint32 i;
618   int           j;
619   int           k;
620   /*We can't unroll the first two iterations of the loop unless _n>=2.*/
621   celt_assert(_n>=2);
622   _u[0]=0;
623   for(k=1;k<=_k+1;k++)_u[k]=(k<<1)-1;
624   i=icwrs1(_y+_n-1,&k);
625   j=_n-2;
626   i+=_u[k];
627   k+=abs(_y[j]);
628   if(_y[j]<0)i+=_u[k+1];
629   while(j-->0){
630     unext(_u,_k+2,0);
631     i+=_u[k];
632     k+=abs(_y[j]);
633     if(_y[j]<0)i+=_u[k+1];
634   }
635   *_nc=_u[k]+_u[k+1];
636   return i;
637 }
638
639 #ifdef CUSTOM_MODES
640 void get_required_bits(opus_int16 *_bits,int _n,int _maxk,int _frac){
641   int k;
642   /*_maxk==0 => there's nothing to do.*/
643   celt_assert(_maxk>0);
644   _bits[0]=0;
645   if (_n==1)
646   {
647     for (k=1;k<=_maxk;k++)
648       _bits[k] = 1<<_frac;
649   }
650   else {
651     VARDECL(opus_uint32,u);
652     SAVE_STACK;
653     ALLOC(u,_maxk+2U,opus_uint32);
654     ncwrs_urow(_n,_maxk,u);
655     for(k=1;k<=_maxk;k++)
656       _bits[k]=log2_frac(u[k]+u[k+1],_frac);
657     RESTORE_STACK;
658   }
659 }
660 #endif /* CUSTOM_MODES */
661
662 void encode_pulses(const int *_y,int _n,int _k,ec_enc *_enc){
663   opus_uint32 i;
664   if (_k==0)
665      return;
666   switch(_n){
667     case 1:{
668       i=icwrs1(_y,&_k);
669       celt_assert(ncwrs1(_k)==2);
670       ec_enc_bits(_enc,i,1);
671     }break;
672 #ifndef SMALL_FOOTPRINT
673     case 2:{
674       i=icwrs2(_y,&_k);
675       ec_enc_uint(_enc,i,ncwrs2(_k));
676     }break;
677     case 3:{
678       i=icwrs3(_y,&_k);
679       ec_enc_uint(_enc,i,ncwrs3(_k));
680     }break;
681     case 4:{
682       i=icwrs4(_y,&_k);
683       ec_enc_uint(_enc,i,ncwrs4(_k));
684     }break;
685     case 5:{
686       i=icwrs5(_y,&_k);
687       ec_enc_uint(_enc,i,ncwrs5(_k));
688     }break;
689 #endif
690      default:
691     {
692       VARDECL(opus_uint32,u);
693       opus_uint32 nc;
694       SAVE_STACK;
695       ALLOC(u,_k+2U,opus_uint32);
696       i=icwrs(_n,_k,&nc,_y,u);
697       ec_enc_uint(_enc,i,nc);
698       RESTORE_STACK;
699     };
700   }
701 }
702
703 void decode_pulses(int *_y,int _n,int _k,ec_dec *_dec)
704 {
705    if (_k==0) {
706       int i;
707       for (i=0;i<_n;i++)
708          _y[i] = 0;
709       return;
710    }
711    switch(_n){
712     case 1:{
713       celt_assert(ncwrs1(_k)==2);
714       cwrsi1(_k,ec_dec_bits(_dec,1),_y);
715     }break;
716 #ifndef SMALL_FOOTPRINT
717     case 2:cwrsi2(_k,ec_dec_uint(_dec,ncwrs2(_k)),_y);break;
718     case 3:cwrsi3(_k,ec_dec_uint(_dec,ncwrs3(_k)),_y);break;
719     case 4:cwrsi4(_k,ec_dec_uint(_dec,ncwrs4(_k)),_y);break;
720     case 5:cwrsi5(_k,ec_dec_uint(_dec,ncwrs5(_k)),_y);break;
721 #endif
722       default:
723     {
724       VARDECL(opus_uint32,u);
725       SAVE_STACK;
726       ALLOC(u,_k+2U,opus_uint32);
727       cwrsi(_n,_k,ec_dec_uint(_dec,ncwrs_urow(_n,_k,u)),_y,u);
728       RESTORE_STACK;
729     }
730   }
731 }
732