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