Further optimization to cwrsi()
[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 COPYRIGHT OWNER
21    OR 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        Doesn't work for val=0, but that case fails the test above.*/
53     if(l>16)val=((val-1)>>(l-16))+1;
54     else val<<=16-l;
55     l=(l-1)<<frac;
56     /*Note that we always need one iteration, since the rounding up above means
57        that we might need to adjust the integer part of the logarithm.*/
58     do{
59       int b;
60       b=(int)(val>>16);
61       l+=b<<frac;
62       val=(val+b)>>b;
63       val=(val*val+0x7FFF)>>15;
64     }
65     while(frac-->0);
66     /*If val is not exactly 0x8000, then we have to round up the remainder.*/
67     return l+(val>0x8000);
68   }
69   /*Exact powers of two require no rounding.*/
70   else return (l-1)<<frac;
71 }
72 #endif
73
74 /*Although derived separately, the pulse vector coding scheme is equivalent to
75    a Pyramid Vector Quantizer \cite{Fis86}.
76   Some additional notes about an early version appear at
77    http://people.xiph.org/~tterribe/notes/cwrs.html, but the codebook ordering
78    and the definitions of some terms have evolved since that was written.
79
80   The conversion from a pulse vector to an integer index (encoding) and back
81    (decoding) is governed by two related functions, V(N,K) and U(N,K).
82
83   V(N,K) = the number of combinations, with replacement, of N items, taken K
84    at a time, when a sign bit is added to each item taken at least once (i.e.,
85    the number of N-dimensional unit pulse vectors with K pulses).
86   One way to compute this is via
87     V(N,K) = K>0 ? sum(k=1...K,2**k*choose(N,k)*choose(K-1,k-1)) : 1,
88    where choose() is the binomial function.
89   A table of values for N<10 and K<10 looks like:
90   V[10][10] = {
91     {1,  0,   0,    0,    0,     0,     0,      0,      0,       0},
92     {1,  2,   2,    2,    2,     2,     2,      2,      2,       2},
93     {1,  4,   8,   12,   16,    20,    24,     28,     32,      36},
94     {1,  6,  18,   38,   66,   102,   146,    198,    258,     326},
95     {1,  8,  32,   88,  192,   360,   608,    952,   1408,    1992},
96     {1, 10,  50,  170,  450,  1002,  1970,   3530,   5890,    9290},
97     {1, 12,  72,  292,  912,  2364,  5336,  10836,  20256,   35436},
98     {1, 14,  98,  462, 1666,  4942, 12642,  28814,  59906,  115598},
99     {1, 16, 128,  688, 2816,  9424, 27008,  68464, 157184,  332688},
100     {1, 18, 162,  978, 4482, 16722, 53154, 148626, 374274,  864146}
101   };
102
103   U(N,K) = the number of such combinations wherein N-1 objects are taken at
104    most K-1 at a time.
105   This is given by
106     U(N,K) = sum(k=0...K-1,V(N-1,k))
107            = K>0 ? (V(N-1,K-1) + V(N,K-1))/2 : 0.
108   The latter expression also makes clear that U(N,K) is half the number of such
109    combinations wherein the first object is taken at least once.
110   Although it may not be clear from either of these definitions, U(N,K) is the
111    natural function to work with when enumerating the pulse vector codebooks,
112    not V(N,K).
113   U(N,K) is not well-defined for N=0, but with the extension
114     U(0,K) = K>0 ? 0 : 1,
115    the function becomes symmetric: U(N,K) = U(K,N), with a similar table:
116   U[10][10] = {
117     {1, 0,  0,   0,    0,    0,     0,     0,      0,      0},
118     {0, 1,  1,   1,    1,    1,     1,     1,      1,      1},
119     {0, 1,  3,   5,    7,    9,    11,    13,     15,     17},
120     {0, 1,  5,  13,   25,   41,    61,    85,    113,    145},
121     {0, 1,  7,  25,   63,  129,   231,   377,    575,    833},
122     {0, 1,  9,  41,  129,  321,   681,  1289,   2241,   3649},
123     {0, 1, 11,  61,  231,  681,  1683,  3653,   7183,  13073},
124     {0, 1, 13,  85,  377, 1289,  3653,  8989,  19825,  40081},
125     {0, 1, 15, 113,  575, 2241,  7183, 19825,  48639, 108545},
126     {0, 1, 17, 145,  833, 3649, 13073, 40081, 108545, 265729}
127   };
128
129   With this extension, V(N,K) may be written in terms of U(N,K):
130     V(N,K) = U(N,K) + U(N,K+1)
131    for all N>=0, K>=0.
132   Thus U(N,K+1) represents the number of combinations where the first element
133    is positive or zero, and U(N,K) represents the number of combinations where
134    it is negative.
135   With a large enough table of U(N,K) values, we could write O(N) encoding
136    and O(min(N*log(K),N+K)) decoding routines, but such a table would be
137    prohibitively large for small embedded devices (K may be as large as 32767
138    for small N, and N may be as large as 200).
139
140   Both functions obey the same recurrence relation:
141     V(N,K) = V(N-1,K) + V(N,K-1) + V(N-1,K-1),
142     U(N,K) = U(N-1,K) + U(N,K-1) + U(N-1,K-1),
143    for all N>0, K>0, with different initial conditions at N=0 or K=0.
144   This allows us to construct a row of one of the tables above given the
145    previous row or the next row.
146   Thus we can derive O(NK) encoding and decoding routines with O(K) memory
147    using only addition and subtraction.
148
149   When encoding, we build up from the U(2,K) row and work our way forwards.
150   When decoding, we need to start at the U(N,K) row and work our way backwards,
151    which requires a means of computing U(N,K).
152   U(N,K) may be computed from two previous values with the same N:
153     U(N,K) = ((2*N-1)*U(N,K-1) - U(N,K-2))/(K-1) + U(N,K-2)
154    for all N>1, and since U(N,K) is symmetric, a similar relation holds for two
155    previous values with the same K:
156     U(N,K>1) = ((2*K-1)*U(N-1,K) - U(N-2,K))/(N-1) + U(N-2,K)
157    for all K>1.
158   This allows us to construct an arbitrary row of the U(N,K) table by starting
159    with the first two values, which are constants.
160   This saves roughly 2/3 the work in our O(NK) decoding routine, but costs O(K)
161    multiplications.
162   Similar relations can be derived for V(N,K), but are not used here.
163
164   For N>0 and K>0, U(N,K) and V(N,K) take on the form of an (N-1)-degree
165    polynomial for fixed N.
166   The first few are
167     U(1,K) = 1,
168     U(2,K) = 2*K-1,
169     U(3,K) = (2*K-2)*K+1,
170     U(4,K) = (((4*K-6)*K+8)*K-3)/3,
171     U(5,K) = ((((2*K-4)*K+10)*K-8)*K+3)/3,
172    and
173     V(1,K) = 2,
174     V(2,K) = 4*K,
175     V(3,K) = 4*K*K+2,
176     V(4,K) = 8*(K*K+2)*K/3,
177     V(5,K) = ((4*K*K+20)*K*K+6)/3,
178    for all K>0.
179   This allows us to derive O(N) encoding and O(N*log(K)) decoding routines for
180    small N (and indeed decoding is also O(N) for N<3).
181
182   @ARTICLE{Fis86,
183     author="Thomas R. Fischer",
184     title="A Pyramid Vector Quantizer",
185     journal="IEEE Transactions on Information Theory",
186     volume="IT-32",
187     number=4,
188     pages="568--583",
189     month=Jul,
190     year=1986
191   }*/
192
193 #if !defined(SMALL_FOOTPRINT)
194
195 /*U(N,K) = U(K,N) := N>0?K>0?U(N-1,K)+U(N,K-1)+U(N-1,K-1):0:K>0?1:0*/
196 # define CELT_PVQ_U(_n,_k) (CELT_PVQ_U_ROW[IMIN(_n,_k)][IMAX(_n,_k)])
197 /*V(N,K) := U(N,K)+U(N,K+1) = the number of PVQ codewords for a band of size N
198    with K pulses allocated to it.*/
199 # define CELT_PVQ_V(_n,_k) (CELT_PVQ_U(_n,_k)+CELT_PVQ_U(_n,(_k)+1))
200
201 /*For each V(N,K) supported, we will access element U(min(N,K+1),max(N,K+1)).
202   Thus, the number of entries in row I is the larger of the maximum number of
203    pulses we will ever allocate for a given N=I (K=128, or however many fit in
204    32 bits, whichever is smaller), plus one, and the maximum N for which
205    K=I-1 pulses fit in 32 bits.
206   The largest band size in an Opus Custom mode is 208.
207   Otherwise, we can limit things to the set of N which can be achieved by
208    splitting a band from a standard Opus mode: 176, 144, 96, 88, 72, 64, 48,
209    44, 36, 32, 24, 22, 18, 16, 8, 4, 2).*/
210 #if defined(CUSTOM_MODES)
211 static const opus_uint32 CELT_PVQ_U_DATA[1488]={
212 #else
213 static const opus_uint32 CELT_PVQ_U_DATA[1272]={
214 #endif
215   /*N=0, K=0...176:*/
216   1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
217   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
218   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
219   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
220   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
221   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
222   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
223 #if defined(CUSTOM_MODES)
224   /*...208:*/
225   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
226   0, 0, 0, 0, 0, 0,
227 #endif
228   /*N=1, K=1...176:*/
229   1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
230   1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
231   1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
232   1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
233   1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
234   1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
235   1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
236 #if defined(CUSTOM_MODES)
237   /*...208:*/
238   1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
239   1, 1, 1, 1, 1, 1,
240 #endif
241   /*N=2, K=2...176:*/
242   3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, 31, 33, 35, 37, 39, 41,
243   43, 45, 47, 49, 51, 53, 55, 57, 59, 61, 63, 65, 67, 69, 71, 73, 75, 77, 79,
244   81, 83, 85, 87, 89, 91, 93, 95, 97, 99, 101, 103, 105, 107, 109, 111, 113,
245   115, 117, 119, 121, 123, 125, 127, 129, 131, 133, 135, 137, 139, 141, 143,
246   145, 147, 149, 151, 153, 155, 157, 159, 161, 163, 165, 167, 169, 171, 173,
247   175, 177, 179, 181, 183, 185, 187, 189, 191, 193, 195, 197, 199, 201, 203,
248   205, 207, 209, 211, 213, 215, 217, 219, 221, 223, 225, 227, 229, 231, 233,
249   235, 237, 239, 241, 243, 245, 247, 249, 251, 253, 255, 257, 259, 261, 263,
250   265, 267, 269, 271, 273, 275, 277, 279, 281, 283, 285, 287, 289, 291, 293,
251   295, 297, 299, 301, 303, 305, 307, 309, 311, 313, 315, 317, 319, 321, 323,
252   325, 327, 329, 331, 333, 335, 337, 339, 341, 343, 345, 347, 349, 351,
253 #if defined(CUSTOM_MODES)
254   /*...208:*/
255   353, 355, 357, 359, 361, 363, 365, 367, 369, 371, 373, 375, 377, 379, 381,
256   383, 385, 387, 389, 391, 393, 395, 397, 399, 401, 403, 405, 407, 409, 411,
257   413, 415,
258 #endif
259   /*N=3, K=3...176:*/
260   13, 25, 41, 61, 85, 113, 145, 181, 221, 265, 313, 365, 421, 481, 545, 613,
261   685, 761, 841, 925, 1013, 1105, 1201, 1301, 1405, 1513, 1625, 1741, 1861,
262   1985, 2113, 2245, 2381, 2521, 2665, 2813, 2965, 3121, 3281, 3445, 3613, 3785,
263   3961, 4141, 4325, 4513, 4705, 4901, 5101, 5305, 5513, 5725, 5941, 6161, 6385,
264   6613, 6845, 7081, 7321, 7565, 7813, 8065, 8321, 8581, 8845, 9113, 9385, 9661,
265   9941, 10225, 10513, 10805, 11101, 11401, 11705, 12013, 12325, 12641, 12961,
266   13285, 13613, 13945, 14281, 14621, 14965, 15313, 15665, 16021, 16381, 16745,
267   17113, 17485, 17861, 18241, 18625, 19013, 19405, 19801, 20201, 20605, 21013,
268   21425, 21841, 22261, 22685, 23113, 23545, 23981, 24421, 24865, 25313, 25765,
269   26221, 26681, 27145, 27613, 28085, 28561, 29041, 29525, 30013, 30505, 31001,
270   31501, 32005, 32513, 33025, 33541, 34061, 34585, 35113, 35645, 36181, 36721,
271   37265, 37813, 38365, 38921, 39481, 40045, 40613, 41185, 41761, 42341, 42925,
272   43513, 44105, 44701, 45301, 45905, 46513, 47125, 47741, 48361, 48985, 49613,
273   50245, 50881, 51521, 52165, 52813, 53465, 54121, 54781, 55445, 56113, 56785,
274   57461, 58141, 58825, 59513, 60205, 60901, 61601,
275 #if defined(CUSTOM_MODES)
276   /*...208:*/
277   62305, 63013, 63725, 64441, 65161, 65885, 66613, 67345, 68081, 68821, 69565,
278   70313, 71065, 71821, 72581, 73345, 74113, 74885, 75661, 76441, 77225, 78013,
279   78805, 79601, 80401, 81205, 82013, 82825, 83641, 84461, 85285, 86113,
280 #endif
281   /*N=4, K=4...176:*/
282   63, 129, 231, 377, 575, 833, 1159, 1561, 2047, 2625, 3303, 4089, 4991, 6017,
283   7175, 8473, 9919, 11521, 13287, 15225, 17343, 19649, 22151, 24857, 27775,
284   30913, 34279, 37881, 41727, 45825, 50183, 54809, 59711, 64897, 70375, 76153,
285   82239, 88641, 95367, 102425, 109823, 117569, 125671, 134137, 142975, 152193,
286   161799, 171801, 182207, 193025, 204263, 215929, 228031, 240577, 253575,
287   267033, 280959, 295361, 310247, 325625, 341503, 357889, 374791, 392217,
288   410175, 428673, 447719, 467321, 487487, 508225, 529543, 551449, 573951,
289   597057, 620775, 645113, 670079, 695681, 721927, 748825, 776383, 804609,
290   833511, 863097, 893375, 924353, 956039, 988441, 1021567, 1055425, 1090023,
291   1125369, 1161471, 1198337, 1235975, 1274393, 1313599, 1353601, 1394407,
292   1436025, 1478463, 1521729, 1565831, 1610777, 1656575, 1703233, 1750759,
293   1799161, 1848447, 1898625, 1949703, 2001689, 2054591, 2108417, 2163175,
294   2218873, 2275519, 2333121, 2391687, 2451225, 2511743, 2573249, 2635751,
295   2699257, 2763775, 2829313, 2895879, 2963481, 3032127, 3101825, 3172583,
296   3244409, 3317311, 3391297, 3466375, 3542553, 3619839, 3698241, 3777767,
297   3858425, 3940223, 4023169, 4107271, 4192537, 4278975, 4366593, 4455399,
298   4545401, 4636607, 4729025, 4822663, 4917529, 5013631, 5110977, 5209575,
299   5309433, 5410559, 5512961, 5616647, 5721625, 5827903, 5935489, 6044391,
300   6154617, 6266175, 6379073, 6493319, 6608921, 6725887, 6844225, 6963943,
301   7085049, 7207551,
302 #if defined(CUSTOM_MODES)
303   /*...208:*/
304   7331457, 7456775, 7583513, 7711679, 7841281, 7972327, 8104825, 8238783,
305   8374209, 8511111, 8649497, 8789375, 8930753, 9073639, 9218041, 9363967,
306   9511425, 9660423, 9810969, 9963071, 10116737, 10271975, 10428793, 10587199,
307   10747201, 10908807, 11072025, 11236863, 11403329, 11571431, 11741177,
308   11912575,
309 #endif
310   /*N=5, K=5...176:*/
311   321, 681, 1289, 2241, 3649, 5641, 8361, 11969, 16641, 22569, 29961, 39041,
312   50049, 63241, 78889, 97281, 118721, 143529, 172041, 204609, 241601, 283401,
313   330409, 383041, 441729, 506921, 579081, 658689, 746241, 842249, 947241,
314   1061761, 1186369, 1321641, 1468169, 1626561, 1797441, 1981449, 2179241,
315   2391489, 2618881, 2862121, 3121929, 3399041, 3694209, 4008201, 4341801,
316   4695809, 5071041, 5468329, 5888521, 6332481, 6801089, 7295241, 7815849,
317   8363841, 8940161, 9545769, 10181641, 10848769, 11548161, 12280841, 13047849,
318   13850241, 14689089, 15565481, 16480521, 17435329, 18431041, 19468809,
319   20549801, 21675201, 22846209, 24064041, 25329929, 26645121, 28010881,
320   29428489, 30899241, 32424449, 34005441, 35643561, 37340169, 39096641,
321   40914369, 42794761, 44739241, 46749249, 48826241, 50971689, 53187081,
322   55473921, 57833729, 60268041, 62778409, 65366401, 68033601, 70781609,
323   73612041, 76526529, 79526721, 82614281, 85790889, 89058241, 92418049,
324   95872041, 99421961, 103069569, 106816641, 110664969, 114616361, 118672641,
325   122835649, 127107241, 131489289, 135983681, 140592321, 145317129, 150160041,
326   155123009, 160208001, 165417001, 170752009, 176215041, 181808129, 187533321,
327   193392681, 199388289, 205522241, 211796649, 218213641, 224775361, 231483969,
328   238341641, 245350569, 252512961, 259831041, 267307049, 274943241, 282741889,
329   290705281, 298835721, 307135529, 315607041, 324252609, 333074601, 342075401,
330   351257409, 360623041, 370174729, 379914921, 389846081, 399970689, 410291241,
331   420810249, 431530241, 442453761, 453583369, 464921641, 476471169, 488234561,
332   500214441, 512413449, 524834241, 537479489, 550351881, 563454121, 576788929,
333   590359041, 604167209, 618216201, 632508801,
334 #if defined(CUSTOM_MODES)
335   /*...208:*/
336   647047809, 661836041, 676876329, 692171521, 707724481, 723538089, 739615241,
337   755958849, 772571841, 789457161, 806617769, 824056641, 841776769, 859781161,
338   878072841, 896654849, 915530241, 934702089, 954173481, 973947521, 994027329,
339   1014416041, 1035116809, 1056132801, 1077467201, 1099123209, 1121104041,
340   1143412929, 1166053121, 1189027881, 1212340489, 1235994241,
341 #endif
342   /*N=6, K=6...96:*/
343   1683, 3653, 7183, 13073, 22363, 36365, 56695, 85305, 124515, 177045, 246047,
344   335137, 448427, 590557, 766727, 982729, 1244979, 1560549, 1937199, 2383409,
345   2908411, 3522221, 4235671, 5060441, 6009091, 7095093, 8332863, 9737793,
346   11326283, 13115773, 15124775, 17372905, 19880915, 22670725, 25765455,
347   29189457, 32968347, 37129037, 41699767, 46710137, 52191139, 58175189,
348   64696159, 71789409, 79491819, 87841821, 96879431, 106646281, 117185651,
349   128542501, 140763503, 153897073, 167993403, 183104493, 199284183, 216588185,
350   235074115, 254801525, 275831935, 298228865, 322057867, 347386557, 374284647,
351   402823977, 433078547, 465124549, 499040399, 534906769, 572806619, 612825229,
352   655050231, 699571641, 746481891, 795875861, 847850911, 902506913, 959946283,
353   1020274013, 1083597703, 1150027593, 1219676595, 1292660325, 1369097135,
354   1449108145, 1532817275, 1620351277, 1711839767, 1807415257, 1907213187,
355   2011371957, 2120032959,
356 #if defined(CUSTOM_MODES)
357   /*...109:*/
358   2233340609U, 2351442379U, 2474488829U, 2602633639U, 2736033641U, 2874848851U,
359   3019242501U, 3169381071U, 3325434321U, 3487575323U, 3655980493U, 3830829623U,
360   4012305913U,
361 #endif
362   /*N=7, K=7...54*/
363   8989, 19825, 40081, 75517, 134245, 227305, 369305, 579125, 880685, 1303777,
364   1884961, 2668525, 3707509, 5064793, 6814249, 9041957, 11847485, 15345233,
365   19665841, 24957661, 31388293, 39146185, 48442297, 59511829, 72616013,
366   88043969, 106114625, 127178701, 151620757, 179861305, 212358985, 249612805,
367   292164445, 340600625, 395555537, 457713341, 527810725, 606639529, 695049433,
368   793950709, 904317037, 1027188385, 1163673953, 1314955181, 1482288821,
369   1667010073, 1870535785, 2094367717,
370 #if defined(CUSTOM_MODES)
371   /*...60:*/
372   2340095869U, 2609401873U, 2904062449U, 3225952925U, 3577050821U, 3959439497U,
373 #endif
374   /*N=8, K=8...37*/
375   48639, 108545, 224143, 433905, 795455, 1392065, 2340495, 3800305, 5984767,
376   9173505, 13726991, 20103025, 28875327, 40754369, 56610575, 77500017,
377   104692735, 139703809, 184327311, 240673265, 311207743, 398796225, 506750351,
378   638878193, 799538175, 993696769, 1226990095, 1505789553, 1837271615,
379   2229491905U,
380 #if defined(CUSTOM_MODES)
381   /*...40:*/
382   2691463695U, 3233240945U, 3866006015U,
383 #endif
384   /*N=9, K=9...28:*/
385   265729, 598417, 1256465, 2485825, 4673345, 8405905, 14546705, 24331777,
386   39490049, 62390545, 96220561, 145198913, 214828609, 312193553, 446304145,
387   628496897, 872893441, 1196924561, 1621925137, 2173806145U,
388 #if defined(CUSTOM_MODES)
389   /*...29:*/
390   2883810113U,
391 #endif
392   /*N=10, K=10...24:*/
393   1462563, 3317445, 7059735, 14218905, 27298155, 50250765, 89129247, 152951073,
394   254831667, 413442773, 654862247, 1014889769, 1541911931, 2300409629U,
395   3375210671U,
396   /*N=11, K=11...19:*/
397   8097453, 18474633, 39753273, 81270333, 158819253, 298199265, 540279585,
398   948062325, 1616336765,
399 #if defined(CUSTOM_MODES)
400   /*...20:*/
401   2684641785U,
402 #endif
403   /*N=12, K=12...18:*/
404   45046719, 103274625, 224298231, 464387817, 921406335, 1759885185,
405   3248227095U,
406   /*N=13, K=13...16:*/
407   251595969, 579168825, 1267854873, 2653649025U,
408   /*N=14, K=14:*/
409   1409933619
410 };
411
412 #if defined(CUSTOM_MODES)
413 const opus_uint32 *const CELT_PVQ_U_ROW[15]={
414   CELT_PVQ_U_DATA+   0,CELT_PVQ_U_DATA+ 208,CELT_PVQ_U_DATA+ 415,
415   CELT_PVQ_U_DATA+ 621,CELT_PVQ_U_DATA+ 826,CELT_PVQ_U_DATA+1030,
416   CELT_PVQ_U_DATA+1233,CELT_PVQ_U_DATA+1336,CELT_PVQ_U_DATA+1389,
417   CELT_PVQ_U_DATA+1421,CELT_PVQ_U_DATA+1441,CELT_PVQ_U_DATA+1455,
418   CELT_PVQ_U_DATA+1464,CELT_PVQ_U_DATA+1470,CELT_PVQ_U_DATA+1473
419 };
420 #else
421 const opus_uint32 *const CELT_PVQ_U_ROW[15]={
422   CELT_PVQ_U_DATA+   0,CELT_PVQ_U_DATA+ 176,CELT_PVQ_U_DATA+ 351,
423   CELT_PVQ_U_DATA+ 525,CELT_PVQ_U_DATA+ 698,CELT_PVQ_U_DATA+ 870,
424   CELT_PVQ_U_DATA+1041,CELT_PVQ_U_DATA+1131,CELT_PVQ_U_DATA+1178,
425   CELT_PVQ_U_DATA+1207,CELT_PVQ_U_DATA+1226,CELT_PVQ_U_DATA+1240,
426   CELT_PVQ_U_DATA+1248,CELT_PVQ_U_DATA+1254,CELT_PVQ_U_DATA+1257
427 };
428 #endif
429
430 #if defined(CUSTOM_MODES)
431 void get_required_bits(opus_int16 *_bits,int _n,int _maxk,int _frac){
432   int k;
433   /*_maxk==0 => there's nothing to do.*/
434   celt_assert(_maxk>0);
435   _bits[0]=0;
436   for(k=1;k<=_maxk;k++)_bits[k]=log2_frac(CELT_PVQ_V(_n,k),_frac);
437 }
438 #endif
439
440 static opus_uint32 icwrs(int _n,const int *_y){
441   opus_uint32 i;
442   int         j;
443   int         k;
444   celt_assert(_n>=2);
445   j=_n-1;
446   i=_y[j]<0;
447   k=abs(_y[j]);
448   do{
449     j--;
450     i+=CELT_PVQ_U(_n-j,k);
451     k+=abs(_y[j]);
452     if(_y[j]<0)i+=CELT_PVQ_U(_n-j,k+1);
453   }
454   while(j>0);
455   return i;
456 }
457
458 void encode_pulses(const int *_y,int _n,int _k,ec_enc *_enc){
459   celt_assert(_k>0);
460   ec_enc_uint(_enc,icwrs(_n,_y),CELT_PVQ_V(_n,_k));
461 }
462
463 static void cwrsi(int _n,int _k,opus_uint32 _i,int *_y){
464   int s;
465   celt_assert(_k>0);
466   celt_assert(_n>1);
467   do{
468     opus_uint32 p;
469     int         k0;
470     /*Are the pulses in this dimension negative?*/
471     p=CELT_PVQ_U(_n,_k+1);
472     s=-(_i>=p);
473     _i-=p&s;
474     /*Count how many pulses were placed in this dimension.*/
475     k0=_k;
476     p=CELT_PVQ_U(_n,_k);
477     if(_k>_n){
478       const opus_uint32 *row;
479       opus_uint32        q;
480       row=CELT_PVQ_U_ROW[_n];
481       q=row[_n];
482       if(q>_i){
483         celt_assert(p>q);
484         /*Setting p=q is unnecessary, but it helps the optimizer prove p>_i,
485            allowing it to jump straight past the initial test in the second
486            loop below.
487           Once it's removed that first comparison, a smart compiler should be
488            able to figure out that the result of this assignment isn't used and
489            optimize it away anyway.*/
490         p=q;
491         _k=_n;
492       }
493       else for(;p>_i;p=row[_k])_k--;
494     }
495     for(;p>_i;p=CELT_PVQ_U_ROW[_k][_n])_k--;
496     _i-=p;
497     *_y++=(k0-_k+s)^s;
498   }
499   while(--_n>1);
500   s=-(_i>=1);
501   *_y=(_k+s)^s;
502 }
503
504 void decode_pulses(int *_y,int _n,int _k,ec_dec *_dec){
505   cwrsi(_n,_k,ec_dec_uint(_dec,CELT_PVQ_V(_n,_k)),_y);
506 }
507
508 #else /* SMALL_FOOTPRINT */
509
510 /*Computes the next row/column of any recurrence that obeys the relation
511    u[i][j]=u[i-1][j]+u[i][j-1]+u[i-1][j-1].
512   _ui0 is the base case for the new row/column.*/
513 static inline void unext(opus_uint32 *_ui,unsigned _len,opus_uint32 _ui0){
514   opus_uint32 ui1;
515   unsigned      j;
516   /*This do-while will overrun the array if we don't have storage for at least
517      2 values.*/
518   j=1; do {
519     ui1=UADD32(UADD32(_ui[j],_ui[j-1]),_ui0);
520     _ui[j-1]=_ui0;
521     _ui0=ui1;
522   } while (++j<_len);
523   _ui[j-1]=_ui0;
524 }
525
526 /*Computes the previous row/column of any recurrence that obeys the relation
527    u[i-1][j]=u[i][j]-u[i][j-1]-u[i-1][j-1].
528   _ui0 is the base case for the new row/column.*/
529 static inline void uprev(opus_uint32 *_ui,unsigned _n,opus_uint32 _ui0){
530   opus_uint32 ui1;
531   unsigned      j;
532   /*This do-while will overrun the array if we don't have storage for at least
533      2 values.*/
534   j=1; do {
535     ui1=USUB32(USUB32(_ui[j],_ui[j-1]),_ui0);
536     _ui[j-1]=_ui0;
537     _ui0=ui1;
538   } while (++j<_n);
539   _ui[j-1]=_ui0;
540 }
541
542 /*Compute V(_n,_k), as well as U(_n,0..._k+1).
543   _u: On exit, _u[i] contains U(_n,i) for i in [0..._k+1].*/
544 static opus_uint32 ncwrs_urow(unsigned _n,unsigned _k,opus_uint32 *_u){
545   opus_uint32 um2;
546   unsigned      len;
547   unsigned      k;
548   len=_k+2;
549   /*We require storage at least 3 values (e.g., _k>0).*/
550   celt_assert(len>=3);
551   _u[0]=0;
552   _u[1]=um2=1;
553   /*If _n==0, _u[0] should be 1 and the rest should be 0.*/
554   /*If _n==1, _u[i] should be 1 for i>1.*/
555   celt_assert(_n>=2);
556   /*If _k==0, the following do-while loop will overflow the buffer.*/
557   celt_assert(_k>0);
558   k=2;
559   do _u[k]=(k<<1)-1;
560   while(++k<len);
561   for(k=2;k<_n;k++)unext(_u+1,_k+1,1);
562   return _u[_k]+_u[_k+1];
563 }
564
565 /*Returns the _i'th combination of _k elements chosen from a set of size _n
566    with associated sign bits.
567   _y: Returns the vector of pulses.
568   _u: Must contain entries [0..._k+1] of row _n of U() on input.
569       Its contents will be destructively modified.*/
570 static void cwrsi(int _n,int _k,opus_uint32 _i,int *_y,opus_uint32 *_u){
571   int j;
572   celt_assert(_n>0);
573   j=0;
574   do{
575     opus_uint32 p;
576     int           s;
577     int           yj;
578     p=_u[_k+1];
579     s=-(_i>=p);
580     _i-=p&s;
581     yj=_k;
582     p=_u[_k];
583     while(p>_i)p=_u[--_k];
584     _i-=p;
585     yj-=_k;
586     _y[j]=(yj+s)^s;
587     uprev(_u,_k+2,0);
588   }
589   while(++j<_n);
590 }
591
592 /*Returns the index of the given combination of K elements chosen from a set
593    of size 1 with associated sign bits.
594   _y: The vector of pulses, whose sum of absolute values is K.
595   _k: Returns K.*/
596 static inline opus_uint32 icwrs1(const int *_y,int *_k){
597   *_k=abs(_y[0]);
598   return _y[0]<0;
599 }
600
601 /*Returns the index of the given combination of K elements chosen from a set
602    of size _n with associated sign bits.
603   _y:  The vector of pulses, whose sum of absolute values must be _k.
604   _nc: Returns V(_n,_k).*/
605 static inline opus_uint32 icwrs(int _n,int _k,opus_uint32 *_nc,const int *_y,
606  opus_uint32 *_u){
607   opus_uint32 i;
608   int         j;
609   int         k;
610   /*We can't unroll the first two iterations of the loop unless _n>=2.*/
611   celt_assert(_n>=2);
612   _u[0]=0;
613   for(k=1;k<=_k+1;k++)_u[k]=(k<<1)-1;
614   i=icwrs1(_y+_n-1,&k);
615   j=_n-2;
616   i+=_u[k];
617   k+=abs(_y[j]);
618   if(_y[j]<0)i+=_u[k+1];
619   while(j-->0){
620     unext(_u,_k+2,0);
621     i+=_u[k];
622     k+=abs(_y[j]);
623     if(_y[j]<0)i+=_u[k+1];
624   }
625   *_nc=_u[k]+_u[k+1];
626   return i;
627 }
628
629 #ifdef CUSTOM_MODES
630 void get_required_bits(opus_int16 *_bits,int _n,int _maxk,int _frac){
631   int k;
632   /*_maxk==0 => there's nothing to do.*/
633   celt_assert(_maxk>0);
634   _bits[0]=0;
635   if (_n==1)
636   {
637     for (k=1;k<=_maxk;k++)
638       _bits[k] = 1<<_frac;
639   }
640   else {
641     VARDECL(opus_uint32,u);
642     SAVE_STACK;
643     ALLOC(u,_maxk+2U,opus_uint32);
644     ncwrs_urow(_n,_maxk,u);
645     for(k=1;k<=_maxk;k++)
646       _bits[k]=log2_frac(u[k]+u[k+1],_frac);
647     RESTORE_STACK;
648   }
649 }
650 #endif /* CUSTOM_MODES */
651
652 void encode_pulses(const int *_y,int _n,int _k,ec_enc *_enc){
653   opus_uint32 i;
654   VARDECL(opus_uint32,u);
655   opus_uint32 nc;
656   SAVE_STACK;
657   celt_assert(_k>0);
658   ALLOC(u,_k+2U,opus_uint32);
659   i=icwrs(_n,_k,&nc,_y,u);
660   ec_enc_uint(_enc,i,nc);
661   RESTORE_STACK;
662 }
663
664 void decode_pulses(int *_y,int _n,int _k,ec_dec *_dec){
665   VARDECL(opus_uint32,u);
666   SAVE_STACK;
667   celt_assert(_k>0);
668   ALLOC(u,_k+2U,opus_uint32);
669   cwrsi(_n,_k,ec_dec_uint(_dec,ncwrs_urow(_n,_k,u)),_y,u);
670   RESTORE_STACK;
671 }
672
673 #endif /* SMALL_FOOTPRINT */