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