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