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