More work on variable frame size (getting rid of FRAMESIZE() )
[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 /*Compute floor(sqrt(_val)) with exact arithmetic.
152   This has been tested on all possible 32-bit inputs.*/
153 static unsigned isqrt32(celt_uint32 _val){
154   unsigned b;
155   unsigned g;
156   int      bshift;
157   /*Uses the second method from
158      http://www.azillionmonkeys.com/qed/sqroot.html
159     The main idea is to search for the largest binary digit b such that
160      (g+b)*(g+b) <= _val, and add it to the solution g.*/
161   g=0;
162   bshift=EC_ILOG(_val)-1>>1;
163   b=1U<<bshift;
164   do{
165     celt_uint32 t;
166     t=((celt_uint32)g<<1)+b<<bshift;
167     if(t<=_val){
168       g+=b;
169       _val-=t;
170     }
171     b>>=1;
172     bshift--;
173   }
174   while(bshift>=0);
175   return g;
176 }
177
178 #endif /* SMALL_FOOTPRINT */
179
180 /*Although derived separately, the pulse vector coding scheme is equivalent to
181    a Pyramid Vector Quantizer \cite{Fis86}.
182   Some additional notes about an early version appear at
183    http://people.xiph.org/~tterribe/notes/cwrs.html, but the codebook ordering
184    and the definitions of some terms have evolved since that was written.
185
186   The conversion from a pulse vector to an integer index (encoding) and back
187    (decoding) is governed by two related functions, V(N,K) and U(N,K).
188
189   V(N,K) = the number of combinations, with replacement, of N items, taken K
190    at a time, when a sign bit is added to each item taken at least once (i.e.,
191    the number of N-dimensional unit pulse vectors with K pulses).
192   One way to compute this is via
193     V(N,K) = K>0 ? sum(k=1...K,2**k*choose(N,k)*choose(K-1,k-1)) : 1,
194    where choose() is the binomial function.
195   A table of values for N<10 and K<10 looks like:
196   V[10][10] = {
197     {1,  0,   0,    0,    0,     0,     0,      0,      0,       0},
198     {1,  2,   2,    2,    2,     2,     2,      2,      2,       2},
199     {1,  4,   8,   12,   16,    20,    24,     28,     32,      36},
200     {1,  6,  18,   38,   66,   102,   146,    198,    258,     326},
201     {1,  8,  32,   88,  192,   360,   608,    952,   1408,    1992},
202     {1, 10,  50,  170,  450,  1002,  1970,   3530,   5890,    9290},
203     {1, 12,  72,  292,  912,  2364,  5336,  10836,  20256,   35436},
204     {1, 14,  98,  462, 1666,  4942, 12642,  28814,  59906,  115598},
205     {1, 16, 128,  688, 2816,  9424, 27008,  68464, 157184,  332688},
206     {1, 18, 162,  978, 4482, 16722, 53154, 148626, 374274,  864146}
207   };
208
209   U(N,K) = the number of such combinations wherein N-1 objects are taken at
210    most K-1 at a time.
211   This is given by
212     U(N,K) = sum(k=0...K-1,V(N-1,k))
213            = K>0 ? (V(N-1,K-1) + V(N,K-1))/2 : 0.
214   The latter expression also makes clear that U(N,K) is half the number of such
215    combinations wherein the first object is taken at least once.
216   Although it may not be clear from either of these definitions, U(N,K) is the
217    natural function to work with when enumerating the pulse vector codebooks,
218    not V(N,K).
219   U(N,K) is not well-defined for N=0, but with the extension
220     U(0,K) = K>0 ? 0 : 1,
221    the function becomes symmetric: U(N,K) = U(K,N), with a similar table:
222   U[10][10] = {
223     {1, 0,  0,   0,    0,    0,     0,     0,      0,      0},
224     {0, 1,  1,   1,    1,    1,     1,     1,      1,      1},
225     {0, 1,  3,   5,    7,    9,    11,    13,     15,     17},
226     {0, 1,  5,  13,   25,   41,    61,    85,    113,    145},
227     {0, 1,  7,  25,   63,  129,   231,   377,    575,    833},
228     {0, 1,  9,  41,  129,  321,   681,  1289,   2241,   3649},
229     {0, 1, 11,  61,  231,  681,  1683,  3653,   7183,  13073},
230     {0, 1, 13,  85,  377, 1289,  3653,  8989,  19825,  40081},
231     {0, 1, 15, 113,  575, 2241,  7183, 19825,  48639, 108545},
232     {0, 1, 17, 145,  833, 3649, 13073, 40081, 108545, 265729}
233   };
234
235   With this extension, V(N,K) may be written in terms of U(N,K):
236     V(N,K) = U(N,K) + U(N,K+1)
237    for all N>=0, K>=0.
238   Thus U(N,K+1) represents the number of combinations where the first element
239    is positive or zero, and U(N,K) represents the number of combinations where
240    it is negative.
241   With a large enough table of U(N,K) values, we could write O(N) encoding
242    and O(min(N*log(K),N+K)) decoding routines, but such a table would be
243    prohibitively large for small embedded devices (K may be as large as 32767
244    for small N, and N may be as large as 200).
245
246   Both functions obey the same recurrence relation:
247     V(N,K) = V(N-1,K) + V(N,K-1) + V(N-1,K-1),
248     U(N,K) = U(N-1,K) + U(N,K-1) + U(N-1,K-1),
249    for all N>0, K>0, with different initial conditions at N=0 or K=0.
250   This allows us to construct a row of one of the tables above given the
251    previous row or the next row.
252   Thus we can derive O(NK) encoding and decoding routines with O(K) memory
253    using only addition and subtraction.
254
255   When encoding, we build up from the U(2,K) row and work our way forwards.
256   When decoding, we need to start at the U(N,K) row and work our way backwards,
257    which requires a means of computing U(N,K).
258   U(N,K) may be computed from two previous values with the same N:
259     U(N,K) = ((2*N-1)*U(N,K-1) - U(N,K-2))/(K-1) + U(N,K-2)
260    for all N>1, and since U(N,K) is symmetric, a similar relation holds for two
261    previous values with the same K:
262     U(N,K>1) = ((2*K-1)*U(N-1,K) - U(N-2,K))/(N-1) + U(N-2,K)
263    for all K>1.
264   This allows us to construct an arbitrary row of the U(N,K) table by starting
265    with the first two values, which are constants.
266   This saves roughly 2/3 the work in our O(NK) decoding routine, but costs O(K)
267    multiplications.
268   Similar relations can be derived for V(N,K), but are not used here.
269
270   For N>0 and K>0, U(N,K) and V(N,K) take on the form of an (N-1)-degree
271    polynomial for fixed N.
272   The first few are
273     U(1,K) = 1,
274     U(2,K) = 2*K-1,
275     U(3,K) = (2*K-2)*K+1,
276     U(4,K) = (((4*K-6)*K+8)*K-3)/3,
277     U(5,K) = ((((2*K-4)*K+10)*K-8)*K+3)/3,
278    and
279     V(1,K) = 2,
280     V(2,K) = 4*K,
281     V(3,K) = 4*K*K+2,
282     V(4,K) = 8*(K*K+2)*K/3,
283     V(5,K) = ((4*K*K+20)*K*K+6)/3,
284    for all K>0.
285   This allows us to derive O(N) encoding and O(N*log(K)) decoding routines for
286    small N (and indeed decoding is also O(N) for N<3).
287
288   @ARTICLE{Fis86,
289     author="Thomas R. Fischer",
290     title="A Pyramid Vector Quantizer",
291     journal="IEEE Transactions on Information Theory",
292     volume="IT-32",
293     number=4,
294     pages="568--583",
295     month=Jul,
296     year=1986
297   }*/
298
299 /*Determines if V(N,K) fits in a 32-bit unsigned integer.
300   N and K are themselves limited to 15 bits.*/
301 int fits_in32(int _n, int _k)
302 {
303    static const celt_int16 maxN[15] = {
304       32767, 32767, 32767, 1476, 283, 109,  60,  40,
305        29,  24,  20,  18,  16,  14,  13};
306    static const celt_int16 maxK[15] = {
307       32767, 32767, 32767, 32767, 1172, 238,  95,  53,
308        36,  27,  22,  18,  16,  15,  13};
309    if (_n>=14)
310    {
311       if (_k>=14)
312          return 0;
313       else
314          return _n <= maxN[_k];
315    } else {
316       return _k <= maxK[_n];
317    }
318 }
319
320 #ifndef SMALL_FOOTPRINT
321
322 /*Compute U(1,_k).*/
323 static inline unsigned ucwrs1(int _k){
324   return _k?1:0;
325 }
326
327 /*Compute V(1,_k).*/
328 static inline unsigned ncwrs1(int _k){
329   return _k?2:1;
330 }
331
332 /*Compute U(2,_k).
333   Note that this may be called with _k=32768 (maxK[2]+1).*/
334 static inline unsigned ucwrs2(unsigned _k){
335   return _k?_k+(_k-1):0;
336 }
337
338 /*Compute V(2,_k).*/
339 static inline celt_uint32 ncwrs2(int _k){
340   return _k?4*(celt_uint32)_k:1;
341 }
342
343 /*Compute U(3,_k).
344   Note that this may be called with _k=32768 (maxK[3]+1).*/
345 static inline celt_uint32 ucwrs3(unsigned _k){
346   return _k?(2*(celt_uint32)_k-2)*_k+1:0;
347 }
348
349 /*Compute V(3,_k).*/
350 static inline celt_uint32 ncwrs3(int _k){
351   return _k?2*(2*(unsigned)_k*(celt_uint32)_k+1):1;
352 }
353
354 /*Compute U(4,_k).*/
355 static inline celt_uint32 ucwrs4(int _k){
356   return _k?imusdiv32odd(2*_k,(2*_k-3)*(celt_uint32)_k+4,3,1):0;
357 }
358
359 /*Compute V(4,_k).*/
360 static inline celt_uint32 ncwrs4(int _k){
361   return _k?((_k*(celt_uint32)_k+2)*_k)/3<<3:1;
362 }
363
364 /*Compute U(5,_k).*/
365 static inline celt_uint32 ucwrs5(int _k){
366   return _k?(((((_k-2)*(unsigned)_k+5)*(celt_uint32)_k-4)*_k)/3<<1)+1:0;
367 }
368
369 /*Compute V(5,_k).*/
370 static inline celt_uint32 ncwrs5(int _k){
371   return _k?(((_k*(unsigned)_k+5)*(celt_uint32)_k*_k)/3<<2)+2:1;
372 }
373
374 #endif /* SMALL_FOOTPRINT */
375
376 /*Computes the next row/column of any recurrence that obeys the relation
377    u[i][j]=u[i-1][j]+u[i][j-1]+u[i-1][j-1].
378   _ui0 is the base case for the new row/column.*/
379 static inline void unext(celt_uint32 *_ui,unsigned _len,celt_uint32 _ui0){
380   celt_uint32 ui1;
381   unsigned      j;
382   /*This do-while will overrun the array if we don't have storage for at least
383      2 values.*/
384   j=1; do {
385     ui1=UADD32(UADD32(_ui[j],_ui[j-1]),_ui0);
386     _ui[j-1]=_ui0;
387     _ui0=ui1;
388   } while (++j<_len);
389   _ui[j-1]=_ui0;
390 }
391
392 /*Computes the previous row/column of any recurrence that obeys the relation
393    u[i-1][j]=u[i][j]-u[i][j-1]-u[i-1][j-1].
394   _ui0 is the base case for the new row/column.*/
395 static inline void uprev(celt_uint32 *_ui,unsigned _n,celt_uint32 _ui0){
396   celt_uint32 ui1;
397   unsigned      j;
398   /*This do-while will overrun the array if we don't have storage for at least
399      2 values.*/
400   j=1; do {
401     ui1=USUB32(USUB32(_ui[j],_ui[j-1]),_ui0);
402     _ui[j-1]=_ui0;
403     _ui0=ui1;
404   } while (++j<_n);
405   _ui[j-1]=_ui0;
406 }
407
408 /*Compute V(_n,_k), as well as U(_n,0..._k+1).
409   _u: On exit, _u[i] contains U(_n,i) for i in [0..._k+1].*/
410 static celt_uint32 ncwrs_urow(unsigned _n,unsigned _k,celt_uint32 *_u){
411   celt_uint32 um2;
412   unsigned      len;
413   unsigned      k;
414   len=_k+2;
415   /*We require storage at least 3 values (e.g., _k>0).*/
416   celt_assert(len>=3);
417   _u[0]=0;
418   _u[1]=um2=1;
419 #ifndef SMALL_FOOTPRINT
420   if(_n<=6 || _k>255)
421 #endif
422  {
423     /*If _n==0, _u[0] should be 1 and the rest should be 0.*/
424     /*If _n==1, _u[i] should be 1 for i>1.*/
425     celt_assert(_n>=2);
426     /*If _k==0, the following do-while loop will overflow the buffer.*/
427     celt_assert(_k>0);
428     k=2;
429     do _u[k]=(k<<1)-1;
430     while(++k<len);
431     for(k=2;k<_n;k++)unext(_u+1,_k+1,1);
432   }
433 #ifndef SMALL_FOOTPRINT
434   else{
435     celt_uint32 um1;
436     celt_uint32 n2m1;
437     _u[2]=n2m1=um1=(_n<<1)-1;
438     for(k=3;k<len;k++){
439       /*U(N,K) = ((2*N-1)*U(N,K-1)-U(N,K-2))/(K-1) + U(N,K-2)*/
440       _u[k]=um2=imusdiv32even(n2m1,um1,um2,k-1)+um2;
441       if(++k>=len)break;
442       _u[k]=um1=imusdiv32odd(n2m1,um2,um1,k-1>>1)+um1;
443     }
444   }
445 #endif /* SMALL_FOOTPRINT */
446   return _u[_k]+_u[_k+1];
447 }
448
449 /*Returns the _i'th combination of _k elements (at most 32767) chosen from a
450    set of size 1 with associated sign bits.
451   _y: Returns the vector of pulses.*/
452 static inline void cwrsi1(int _k,celt_uint32 _i,int *_y){
453   int s;
454   s=-(int)_i;
455   _y[0]=_k+s^s;
456 }
457
458 #ifndef SMALL_FOOTPRINT
459
460 /*Returns the _i'th combination of _k elements (at most 32767) chosen from a
461    set of size 2 with associated sign bits.
462   _y: Returns the vector of pulses.*/
463 static inline void cwrsi2(int _k,celt_uint32 _i,int *_y){
464   celt_uint32 p;
465   int           s;
466   int           yj;
467   p=ucwrs2(_k+1U);
468   s=-(_i>=p);
469   _i-=p&s;
470   yj=_k;
471   _k=_i+1>>1;
472   p=ucwrs2(_k);
473   _i-=p;
474   yj-=_k;
475   _y[0]=yj+s^s;
476   cwrsi1(_k,_i,_y+1);
477 }
478
479 /*Returns the _i'th combination of _k elements (at most 32767) chosen from a
480    set of size 3 with associated sign bits.
481   _y: Returns the vector of pulses.*/
482 static void cwrsi3(int _k,celt_uint32 _i,int *_y){
483   celt_uint32 p;
484   int           s;
485   int           yj;
486   p=ucwrs3(_k+1U);
487   s=-(_i>=p);
488   _i-=p&s;
489   yj=_k;
490   /*Finds the maximum _k such that ucwrs3(_k)<=_i (tested for all
491      _i<2147418113=U(3,32768)).*/
492   _k=_i>0?isqrt32(2*_i-1)+1>>1:0;
493   p=ucwrs3(_k);
494   _i-=p;
495   yj-=_k;
496   _y[0]=yj+s^s;
497   cwrsi2(_k,_i,_y+1);
498 }
499
500 /*Returns the _i'th combination of _k elements (at most 1172) chosen from a set
501    of size 4 with associated sign bits.
502   _y: Returns the vector of pulses.*/
503 static void cwrsi4(int _k,celt_uint32 _i,int *_y){
504   celt_uint32 p;
505   int           s;
506   int           yj;
507   int           kl;
508   int           kr;
509   p=ucwrs4(_k+1);
510   s=-(_i>=p);
511   _i-=p&s;
512   yj=_k;
513   /*We could solve a cubic for k here, but the form of the direct solution does
514      not lend itself well to exact integer arithmetic.
515     Instead we do a binary search on U(4,K).*/
516   kl=0;
517   kr=_k;
518   for(;;){
519     _k=kl+kr>>1;
520     p=ucwrs4(_k);
521     if(p<_i){
522       if(_k>=kr)break;
523       kl=_k+1;
524     }
525     else if(p>_i)kr=_k-1;
526     else break;
527   }
528   _i-=p;
529   yj-=_k;
530   _y[0]=yj+s^s;
531   cwrsi3(_k,_i,_y+1);
532 }
533
534 /*Returns the _i'th combination of _k elements (at most 238) chosen from a set
535    of size 5 with associated sign bits.
536   _y: Returns the vector of pulses.*/
537 static void cwrsi5(int _k,celt_uint32 _i,int *_y){
538   celt_uint32 p;
539   int           s;
540   int           yj;
541   p=ucwrs5(_k+1);
542   s=-(_i>=p);
543   _i-=p&s;
544   yj=_k;
545   /* A binary search on U(5,K) avoids the need for 64-bit arithmetic */
546   {
547     int kl=0;
548     int kr=_k;
549     for(;;){
550       _k=kl+kr>>1;
551       p=ucwrs5(_k);
552       if(p<_i){
553         if(_k>=kr)break;
554         kl=_k+1;
555       }
556       else if(p>_i)kr=_k-1;
557       else break;
558     }  
559   }
560   _i-=p;
561   yj-=_k;
562   _y[0]=yj+s^s;
563   cwrsi4(_k,_i,_y+1);
564 }
565 #endif /* SMALL_FOOTPRINT */
566
567 /*Returns the _i'th combination of _k elements chosen from a set of size _n
568    with associated sign bits.
569   _y: Returns the vector of pulses.
570   _u: Must contain entries [0..._k+1] of row _n of U() on input.
571       Its contents will be destructively modified.*/
572 static void cwrsi(int _n,int _k,celt_uint32 _i,int *_y,celt_uint32 *_u){
573   int j;
574   celt_assert(_n>0);
575   j=0;
576   do{
577     celt_uint32 p;
578     int           s;
579     int           yj;
580     p=_u[_k+1];
581     s=-(_i>=p);
582     _i-=p&s;
583     yj=_k;
584     p=_u[_k];
585     while(p>_i)p=_u[--_k];
586     _i-=p;
587     yj-=_k;
588     _y[j]=yj+s^s;
589     uprev(_u,_k+2,0);
590   }
591   while(++j<_n);
592 }
593
594
595 /*Returns the index of the given combination of K elements chosen from a set
596    of size 1 with associated sign bits.
597   _y: The vector of pulses, whose sum of absolute values is K.
598   _k: Returns K.*/
599 static inline celt_uint32 icwrs1(const int *_y,int *_k){
600   *_k=abs(_y[0]);
601   return _y[0]<0;
602 }
603
604 #ifndef SMALL_FOOTPRINT
605
606 /*Returns the index of the given combination of K elements chosen from a set
607    of size 2 with associated sign bits.
608   _y: The vector of pulses, whose sum of absolute values is K.
609   _k: Returns K.*/
610 static inline celt_uint32 icwrs2(const int *_y,int *_k){
611   celt_uint32 i;
612   int           k;
613   i=icwrs1(_y+1,&k);
614   i+=ucwrs2(k);
615   k+=abs(_y[0]);
616   if(_y[0]<0)i+=ucwrs2(k+1U);
617   *_k=k;
618   return i;
619 }
620
621 /*Returns the index of the given combination of K elements chosen from a set
622    of size 3 with associated sign bits.
623   _y: The vector of pulses, whose sum of absolute values is K.
624   _k: Returns K.*/
625 static inline celt_uint32 icwrs3(const int *_y,int *_k){
626   celt_uint32 i;
627   int           k;
628   i=icwrs2(_y+1,&k);
629   i+=ucwrs3(k);
630   k+=abs(_y[0]);
631   if(_y[0]<0)i+=ucwrs3(k+1U);
632   *_k=k;
633   return i;
634 }
635
636 /*Returns the index of the given combination of K elements chosen from a set
637    of size 4 with associated sign bits.
638   _y: The vector of pulses, whose sum of absolute values is K.
639   _k: Returns K.*/
640 static inline celt_uint32 icwrs4(const int *_y,int *_k){
641   celt_uint32 i;
642   int           k;
643   i=icwrs3(_y+1,&k);
644   i+=ucwrs4(k);
645   k+=abs(_y[0]);
646   if(_y[0]<0)i+=ucwrs4(k+1);
647   *_k=k;
648   return i;
649 }
650
651 /*Returns the index of the given combination of K elements chosen from a set
652    of size 5 with associated sign bits.
653   _y: The vector of pulses, whose sum of absolute values is K.
654   _k: Returns K.*/
655 static inline celt_uint32 icwrs5(const int *_y,int *_k){
656   celt_uint32 i;
657   int           k;
658   i=icwrs4(_y+1,&k);
659   i+=ucwrs5(k);
660   k+=abs(_y[0]);
661   if(_y[0]<0)i+=ucwrs5(k+1);
662   *_k=k;
663   return i;
664 }
665 #endif /* SMALL_FOOTPRINT */
666
667 /*Returns the index of the given combination of K elements chosen from a set
668    of size _n with associated sign bits.
669   _y:  The vector of pulses, whose sum of absolute values must be _k.
670   _nc: Returns V(_n,_k).*/
671 celt_uint32 icwrs(int _n,int _k,celt_uint32 *_nc,const int *_y,
672  celt_uint32 *_u){
673   celt_uint32 i;
674   int           j;
675   int           k;
676   /*We can't unroll the first two iterations of the loop unless _n>=2.*/
677   celt_assert(_n>=2);
678   _u[0]=0;
679   for(k=1;k<=_k+1;k++)_u[k]=(k<<1)-1;
680   i=icwrs1(_y+_n-1,&k);
681   j=_n-2;
682   i+=_u[k];
683   k+=abs(_y[j]);
684   if(_y[j]<0)i+=_u[k+1];
685   while(j-->0){
686     unext(_u,_k+2,0);
687     i+=_u[k];
688     k+=abs(_y[j]);
689     if(_y[j]<0)i+=_u[k+1];
690   }
691   *_nc=_u[k]+_u[k+1];
692   return i;
693 }
694
695
696 /*Computes get_required_bits when splitting is required.
697   _left_bits and _right_bits must contain the required bits for the left and
698    right sides of the split, respectively (which themselves may require
699    splitting).*/
700 static void get_required_split_bits(celt_int16 *_bits,
701  const celt_int16 *_left_bits,const celt_int16 *_right_bits,
702  int _n,int _maxk,int _frac){
703   int k;
704   for(k=_maxk;k-->0;){
705     /*If we've reached a k where everything fits in 32 bits, evaluate the
706        remaining required bits directly.*/
707     if(fits_in32(_n,k)){
708       get_required_bits(_bits,_n,k+1,_frac);
709       break;
710     }
711     else{
712       int worst_bits;
713       int i;
714       /*Due to potentially recursive splitting, it's difficult to derive an
715          analytic expression for the location of the worst-case split index.
716         We simply check them all.*/
717       worst_bits=0;
718       for(i=0;i<=k;i++){
719         int split_bits;
720         split_bits=_left_bits[i]+_right_bits[k-i];
721         if(split_bits>worst_bits)worst_bits=split_bits;
722       }
723       _bits[k]=log2_frac(k+1,_frac)+worst_bits;
724     }
725   }
726 }
727
728 /*Computes get_required_bits for a pair of N values.
729   _n1 and _n2 must either be equal or two consecutive integers.
730   Returns the buffer used to store the required bits for _n2, which is either
731    _bits1 if _n1==_n2 or _bits2 if _n1+1==_n2.*/
732 static celt_int16 *get_required_bits_pair(celt_int16 *_bits1,
733  celt_int16 *_bits2,celt_int16 *_tmp,int _n1,int _n2,int _maxk,int _frac){
734   celt_int16 *tmp2;
735   /*If we only need a single set of required bits...*/
736   if(_n1==_n2){
737     /*Stop recursing if everything fits.*/
738     if(fits_in32(_n1,_maxk-1))get_required_bits(_bits1,_n1,_maxk,_frac);
739     else{
740       _tmp=get_required_bits_pair(_bits2,_tmp,_bits1,
741        _n1>>1,_n1+1>>1,_maxk,_frac);
742       get_required_split_bits(_bits1,_bits2,_tmp,_n1,_maxk,_frac);
743     }
744     return _bits1;
745   }
746   /*Otherwise we need two distinct sets...*/
747   celt_assert(_n1+1==_n2);
748   /*Stop recursing if everything fits.*/
749   if(fits_in32(_n2,_maxk-1)){
750     get_required_bits(_bits1,_n1,_maxk,_frac);
751     get_required_bits(_bits2,_n2,_maxk,_frac);
752   }
753   /*Otherwise choose an evaluation order that doesn't require extra buffers.*/
754   else if(_n1&1){
755     /*This special case isn't really needed, but can save some work.*/
756     if(fits_in32(_n1,_maxk-1)){
757       tmp2=get_required_bits_pair(_tmp,_bits1,_bits2,
758        _n2>>1,_n2>>1,_maxk,_frac);
759       get_required_split_bits(_bits2,_tmp,tmp2,_n2,_maxk,_frac);
760       get_required_bits(_bits1,_n1,_maxk,_frac);
761     }
762     else{
763       _tmp=get_required_bits_pair(_bits2,_tmp,_bits1,
764        _n1>>1,_n1+1>>1,_maxk,_frac);
765       get_required_split_bits(_bits1,_bits2,_tmp,_n1,_maxk,_frac);
766       get_required_split_bits(_bits2,_tmp,_tmp,_n2,_maxk,_frac);
767     }
768   }
769   else{
770     /*There's no need to special case _n1 fitting by itself, since _n2 requires
771        us to recurse for both values anyway.*/
772     tmp2=get_required_bits_pair(_tmp,_bits1,_bits2,
773      _n2>>1,_n2+1>>1,_maxk,_frac);
774     get_required_split_bits(_bits2,_tmp,tmp2,_n2,_maxk,_frac);
775     get_required_split_bits(_bits1,_tmp,_tmp,_n1,_maxk,_frac);
776   }
777   return _bits2;
778 }
779
780 void get_required_bits(celt_int16 *_bits,int _n,int _maxk,int _frac){
781   int k;
782   /*_maxk==0 => there's nothing to do.*/
783   celt_assert(_maxk>0);
784   if (_n==1)
785   {
786     _bits[0] = 0;
787     for (k=1;k<_maxk;k++)
788       _bits[k] = 1<<_frac;
789   }
790   else if(fits_in32(_n,_maxk-1)){
791     _bits[0]=0;
792     if(_maxk>1){
793       VARDECL(celt_uint32,u);
794       SAVE_STACK;
795       ALLOC(u,_maxk+1U,celt_uint32);
796       ncwrs_urow(_n,_maxk-1,u);
797       for(k=1;k<_maxk;k++)_bits[k]=log2_frac(u[k]+u[k+1],_frac);
798       RESTORE_STACK;
799     }
800   }
801   else{
802     VARDECL(celt_int16,n1bits);
803     VARDECL(celt_int16,n2bits_buf);
804     celt_int16 *n2bits;
805     SAVE_STACK;
806     ALLOC(n1bits,_maxk,celt_int16);
807     ALLOC(n2bits_buf,_maxk,celt_int16);
808     n2bits=get_required_bits_pair(n1bits,n2bits_buf,_bits,
809      _n>>1,_n+1>>1,_maxk,_frac);
810     get_required_split_bits(_bits,n1bits,n2bits,_n,_maxk,_frac);
811     RESTORE_STACK;
812   }
813 }
814
815
816 static inline void encode_pulses32(int _n,int _k,const int *_y,ec_enc *_enc){
817   celt_uint32 i;
818   switch(_n){
819     case 1:{
820       i=icwrs1(_y,&_k);
821       celt_assert(ncwrs1(_k)==2);
822       ec_enc_bits(_enc,i,1);
823     }break;
824 #ifndef SMALL_FOOTPRINT
825     case 2:{
826       i=icwrs2(_y,&_k);
827       ec_enc_uint(_enc,i,ncwrs2(_k));
828     }break;
829     case 3:{
830       i=icwrs3(_y,&_k);
831       ec_enc_uint(_enc,i,ncwrs3(_k));
832     }break;
833     case 4:{
834       i=icwrs4(_y,&_k);
835       ec_enc_uint(_enc,i,ncwrs4(_k));
836     }break;
837     case 5:{
838       i=icwrs5(_y,&_k);
839       ec_enc_uint(_enc,i,ncwrs5(_k));
840     }break;
841 #endif
842      default:
843     {
844       VARDECL(celt_uint32,u);
845       celt_uint32 nc;
846       SAVE_STACK;
847       ALLOC(u,_k+2U,celt_uint32);
848       i=icwrs(_n,_k,&nc,_y,u);
849       ec_enc_uint(_enc,i,nc);
850       RESTORE_STACK;
851     };
852   }
853 }
854
855 void encode_pulses(int *_y, int N, int K, ec_enc *enc)
856 {
857    if (K==0) {
858    } else if(fits_in32(N,K))
859    {
860       encode_pulses32(N, K, _y, enc);
861    } else {
862      int i;
863      int count=0;
864      int split;
865      split = (N+1)/2;
866      for (i=0;i<split;i++)
867         count += abs(_y[i]);
868      ec_enc_uint(enc,count,K+1);
869      encode_pulses(_y, split, count, enc);
870      encode_pulses(_y+split, N-split, K-count, enc);
871    }
872 }
873
874 static inline void decode_pulses32(int _n,int _k,int *_y,ec_dec *_dec){
875    switch(_n){
876     case 1:{
877       celt_assert(ncwrs1(_k)==2);
878       cwrsi1(_k,ec_dec_bits(_dec,1),_y);
879     }break;
880 #ifndef SMALL_FOOTPRINT
881     case 2:cwrsi2(_k,ec_dec_uint(_dec,ncwrs2(_k)),_y);break;
882     case 3:cwrsi3(_k,ec_dec_uint(_dec,ncwrs3(_k)),_y);break;
883     case 4:cwrsi4(_k,ec_dec_uint(_dec,ncwrs4(_k)),_y);break;
884     case 5:cwrsi5(_k,ec_dec_uint(_dec,ncwrs5(_k)),_y);break;
885 #endif
886       default:
887     {
888       VARDECL(celt_uint32,u);
889       SAVE_STACK;
890       ALLOC(u,_k+2U,celt_uint32);
891       cwrsi(_n,_k,ec_dec_uint(_dec,ncwrs_urow(_n,_k,u)),_y,u);
892       RESTORE_STACK;
893     }
894   }
895 }
896
897 void decode_pulses(int *_y, int N, int K, ec_dec *dec)
898 {
899    if (K==0) {
900       int i;
901       for (i=0;i<N;i++)
902          _y[i] = 0;
903    } else if(fits_in32(N,K))
904    {
905       decode_pulses32(N, K, _y, dec);
906    } else {
907      int split;
908      int count = ec_dec_uint(dec,K+1);
909      split = (N+1)/2;
910      decode_pulses(_y, split, count, dec);
911      decode_pulses(_y+split, N-split, K-count, dec);
912    }
913 }