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