lib/sharedbook.c

   1 /********************************************************************
   2  *                                                                  *
   3  * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE.   *
   4  * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS     *
   5  * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
   6  * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.       *
   7  *                                                                  *
   8  * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2015             *
   9  * by the Xiph.Org Foundation http://www.xiph.org/                  *
  10  *                                                                  *
  11  ********************************************************************
  12
  13  function: basic shared codebook operations
  14  last mod: $Id$
  15
  16  ********************************************************************/
  17
  18 #include <stdlib.h>
  19 #include <limits.h>
  20 #include <math.h>
  21 #include <string.h>
  22 #include <ogg/ogg.h>
  23 #include "os.h"
  24 #include "misc.h"
  25 #include "vorbis/codec.h"
  26 #include "codebook.h"
  27 #include "scales.h"
  28
  29 /**** pack/unpack helpers ******************************************/
  30
  31 int ov_ilog(ogg_uint32_t v){
  32   int ret;
  33   for(ret=0;v;ret++)v>>=1;
  34   return ret;
  35 }
  36
  37 /* 32 bit float (not IEEE; nonnormalized mantissa +
  38    biased exponent) : neeeeeee eeemmmmm mmmmmmmm mmmmmmmm
  39    Why not IEEE?  It's just not that important here. */
  40
  41 #define VQ_FEXP 10
  42 #define VQ_FMAN 21
  43 #define VQ_FEXP_BIAS 768 /* bias toward values smaller than 1. */
  44
  45 /* doesn't currently guard under/overflow */
  46 long _float32_pack(float val){
  47   int sign=0;
  48   long exp;
  49   long mant;
  50   if(val<0){
  51     sign=0x80000000;
  52     val= -val;
  53   }
  54   exp= floor(log(val)/log(2.f)+.001); //+epsilon
  55   mant=rint(ldexp(val,(VQ_FMAN-1)-exp));
  56   exp=(exp+VQ_FEXP_BIAS)<<VQ_FMAN;
  57
  58   return(sign|exp|mant);
  59 }
  60
  61 float _float32_unpack(long val){
  62   double mant=val&0x1fffff;
  63   int    sign=val&0x80000000;
  64   long   exp =(val&0x7fe00000L)>>VQ_FMAN;
  65   if(sign)mant= -mant;
  66   return(ldexp(mant,exp-(VQ_FMAN-1)-VQ_FEXP_BIAS));
  67 }
  68
  69 /* given a list of word lengths, generate a list of codewords.  Works
  70    for length ordered or unordered, always assigns the lowest valued
  71    codewords first.  Extended to handle unused entries (length 0) */
  72 ogg_uint32_t *_make_words(char *l,long n,long sparsecount){
  73   long i,j,count=0;
  74   ogg_uint32_t marker[33];
  75   ogg_uint32_t *r=_ogg_malloc((sparsecount?sparsecount:n)*sizeof(*r));
  76   memset(marker,0,sizeof(marker));
  77
  78   for(i=0;i<n;i++){
  79     long length=l[i];
  80     if(length>0){
  81       ogg_uint32_t entry=marker[length];
  82
  83       /* when we claim a node for an entry, we also claim the nodes
  84          below it (pruning off the imagined tree that may have dangled
  85          from it) as well as blocking the use of any nodes directly
  86          above for leaves */
  87
  88       /* update ourself */
  89       if(length<32 && (entry>>length)){
  90         /* error condition; the lengths must specify an overpopulated tree */
  91         _ogg_free(r);
  92         return(NULL);
  93       }
  94       r[count++]=entry;
  95
  96       /* Look to see if the next shorter marker points to the node
  97          above. if so, update it and repeat.  */
  98       {
  99         for(j=length;j>0;j--){
 100
 101           if(marker[j]&1){
 102             /* have to jump branches */
 103             if(j==1)
 104               marker[1]++;
 105             else
 106               marker[j]=marker[j-1]<<1;
 107             break; /* invariant says next upper marker would already
 108                       have been moved if it was on the same path */
 109           }
 110           marker[j]++;
 111         }
 112       }
 113
 114       /* prune the tree; the implicit invariant says all the longer
 115          markers were dangling from our just-taken node.  Dangle them
 116          from our *new* node. */
 117       for(j=length+1;j<33;j++)
 118         if((marker[j]>>1) == entry){
 119           entry=marker[j];
 120           marker[j]=marker[j-1]<<1;
 121         }else
 122           break;
 123     }else
 124       if(sparsecount==0)count++;
 125   }
 126
 127   /* any underpopulated tree must be rejected. */
 128   /* Single-entry codebooks are a retconned extension to the spec.
 129      They have a single codeword '0' of length 1 that results in an
 130      underpopulated tree.  Shield that case from the underformed tree check. */
 131   if(!(count==1 && marker[2]==2)){
 132     for(i=1;i<33;i++)
 133       if(marker[i] & (0xffffffffUL>>(32-i))){
 134         _ogg_free(r);
 135         return(NULL);
 136       }
 137   }
 138
 139   /* bitreverse the words because our bitwise packer/unpacker is LSb
 140      endian */
 141   for(i=0,count=0;i<n;i++){
 142     ogg_uint32_t temp=0;
 143     for(j=0;j<l[i];j++){
 144       temp<<=1;
 145       temp|=(r[count]>>j)&1;
 146     }
 147
 148     if(sparsecount){
 149       if(l[i])
 150         r[count++]=temp;
 151     }else
 152       r[count++]=temp;
 153   }
 154
 155   return(r);
 156 }
 157
 158 /* there might be a straightforward one-line way to do the below
 159    that's portable and totally safe against roundoff, but I haven't
 160    thought of it.  Therefore, we opt on the side of caution */
 161 long _book_maptype1_quantvals(const static_codebook *b){
 162   long vals;
 163   if(b->entries<1){
 164     return(0);
 165   }
 166   vals=floor(pow((float)b->entries,1.f/b->dim));
 167
 168   /* the above *should* be reliable, but we'll not assume that FP is
 169      ever reliable when bitstream sync is at stake; verify via integer
 170      means that vals really is the greatest value of dim for which
 171      vals^b->bim <= b->entries */
 172   /* treat the above as an initial guess */
 173   if(vals<1){
 174     vals=1;
 175   }
 176   while(1){
 177     long acc=1;
 178     long acc1=1;
 179     int i;
 180     for(i=0;i<b->dim;i++){
 181       if(b->entries/vals<acc)break;
 182       acc*=vals;
 183       if(LONG_MAX/(vals+1)<acc1)acc1=LONG_MAX;
 184       else acc1*=vals+1;
 185     }
 186     if(i>=b->dim && acc<=b->entries && acc1>b->entries){
 187       return(vals);
 188     }else{
 189       if(i<b->dim || acc>b->entries){
 190         vals--;
 191       }else{
 192         vals++;
 193       }
 194     }
 195   }
 196 }
 197
 198 /* unpack the quantized list of values for encode/decode ***********/
 199 /* we need to deal with two map types: in map type 1, the values are
 200    generated algorithmically (each column of the vector counts through
 201    the values in the quant vector). in map type 2, all the values came
 202    in in an explicit list.  Both value lists must be unpacked */
 203 float *_book_unquantize(const static_codebook *b,int n,int *sparsemap){
 204   long j,k,count=0;
 205   if(b->maptype==1 || b->maptype==2){
 206     int quantvals;
 207     float mindel=_float32_unpack(b->q_min);
 208     float delta=_float32_unpack(b->q_delta);
 209     float *r=_ogg_calloc(n*b->dim,sizeof(*r));
 210
 211     /* maptype 1 and 2 both use a quantized value vector, but
 212        different sizes */
 213     switch(b->maptype){
 214     case 1:
 215       /* most of the time, entries%dimensions == 0, but we need to be
 216          well defined.  We define that the possible vales at each
 217          scalar is values == entries/dim.  If entries%dim != 0, we'll
 218          have 'too few' values (values*dim<entries), which means that
 219          we'll have 'left over' entries; left over entries use zeroed
 220          values (and are wasted).  So don't generate codebooks like
 221          that */
 222       quantvals=_book_maptype1_quantvals(b);
 223       for(j=0;j<b->entries;j++){
 224         if((sparsemap && b->lengthlist[j]) || !sparsemap){
 225           float last=0.f;
 226           int indexdiv=1;
 227           for(k=0;k<b->dim;k++){
 228             int index= (j/indexdiv)%quantvals;
 229             float val=b->quantlist[index];
 230             val=fabs(val)*delta+mindel+last;
 231             if(b->q_sequencep)last=val;
 232             if(sparsemap)
 233               r[sparsemap[count]*b->dim+k]=val;
 234             else
 235               r[count*b->dim+k]=val;
 236             indexdiv*=quantvals;
 237           }
 238           count++;
 239         }
 240
 241       }
 242       break;
 243     case 2:
 244       for(j=0;j<b->entries;j++){
 245         if((sparsemap && b->lengthlist[j]) || !sparsemap){
 246           float last=0.f;
 247
 248           for(k=0;k<b->dim;k++){
 249             float val=b->quantlist[j*b->dim+k];
 250             val=fabs(val)*delta+mindel+last;
 251             if(b->q_sequencep)last=val;
 252             if(sparsemap)
 253               r[sparsemap[count]*b->dim+k]=val;
 254             else
 255               r[count*b->dim+k]=val;
 256           }
 257           count++;
 258         }
 259       }
 260       break;
 261     }
 262
 263     return(r);
 264   }
 265   return(NULL);
 266 }
 267
 268 void vorbis_staticbook_destroy(static_codebook *b){
 269   if(b->allocedp){
 270     if(b->quantlist)_ogg_free(b->quantlist);
 271     if(b->lengthlist)_ogg_free(b->lengthlist);
 272     memset(b,0,sizeof(*b));
 273     _ogg_free(b);
 274   } /* otherwise, it is in static memory */
 275 }
 276
 277 void vorbis_book_clear(codebook *b){
 278   /* static book is not cleared; we're likely called on the lookup and
 279      the static codebook belongs to the info struct */
 280   if(b->valuelist)_ogg_free(b->valuelist);
 281   if(b->codelist)_ogg_free(b->codelist);
 282
 283   if(b->dec_index)_ogg_free(b->dec_index);
 284   if(b->dec_codelengths)_ogg_free(b->dec_codelengths);
 285   if(b->dec_firsttable)_ogg_free(b->dec_firsttable);
 286
 287   memset(b,0,sizeof(*b));
 288 }
 289
 290 int vorbis_book_init_encode(codebook *c,const static_codebook *s){
 291
 292   memset(c,0,sizeof(*c));
 293   c->c=s;
 294   c->entries=s->entries;
 295   c->used_entries=s->entries;
 296   c->dim=s->dim;
 297   c->codelist=_make_words(s->lengthlist,s->entries,0);
 298   //c->valuelist=_book_unquantize(s,s->entries,NULL);
 299   c->quantvals=_book_maptype1_quantvals(s);
 300   c->minval=(int)rint(_float32_unpack(s->q_min));
 301   c->delta=(int)rint(_float32_unpack(s->q_delta));
 302
 303   return(0);
 304 }
 305
 306 static ogg_uint32_t bitreverse(ogg_uint32_t x){
 307   x=    ((x>>16)&0x0000ffffUL) | ((x<<16)&0xffff0000UL);
 308   x=    ((x>> 8)&0x00ff00ffUL) | ((x<< 8)&0xff00ff00UL);
 309   x=    ((x>> 4)&0x0f0f0f0fUL) | ((x<< 4)&0xf0f0f0f0UL);
 310   x=    ((x>> 2)&0x33333333UL) | ((x<< 2)&0xccccccccUL);
 311   return((x>> 1)&0x55555555UL) | ((x<< 1)&0xaaaaaaaaUL);
 312 }
 313
 314 static int sort32a(const void *a,const void *b){
 315   return ( **(ogg_uint32_t **)a>**(ogg_uint32_t **)b)-
 316     ( **(ogg_uint32_t **)a<**(ogg_uint32_t **)b);
 317 }
 318
 319 /* decode codebook arrangement is more heavily optimized than encode */
 320 int vorbis_book_init_decode(codebook *c,const static_codebook *s){
 321   int i,j,n=0,tabn;
 322   int *sortindex;
 323
 324   memset(c,0,sizeof(*c));
 325
 326   /* count actually used entries and find max length */
 327   for(i=0;i<s->entries;i++)
 328     if(s->lengthlist[i]>0)
 329       n++;
 330
 331   c->entries=s->entries;
 332   c->used_entries=n;
 333   c->dim=s->dim;
 334
 335   if(n>0){
 336     /* two different remappings go on here.
 337
 338     First, we collapse the likely sparse codebook down only to
 339     actually represented values/words.  This collapsing needs to be
 340     indexed as map-valueless books are used to encode original entry
 341     positions as integers.
 342
 343     Second, we reorder all vectors, including the entry index above,
 344     by sorted bitreversed codeword to allow treeless decode. */
 345
 346     /* perform sort */
 347     ogg_uint32_t *codes=_make_words(s->lengthlist,s->entries,c->used_entries);
 348     ogg_uint32_t **codep=alloca(sizeof(*codep)*n);
 349
 350     if(codes==NULL)goto err_out;
 351
 352     for(i=0;i<n;i++){
 353       codes[i]=bitreverse(codes[i]);
 354       codep[i]=codes+i;
 355     }
 356
 357     qsort(codep,n,sizeof(*codep),sort32a);
 358
 359     sortindex=alloca(n*sizeof(*sortindex));
 360     c->codelist=_ogg_malloc(n*sizeof(*c->codelist));
 361     /* the index is a reverse index */
 362     for(i=0;i<n;i++){
 363       int position=codep[i]-codes;
 364       sortindex[position]=i;
 365     }
 366
 367     for(i=0;i<n;i++)
 368       c->codelist[sortindex[i]]=codes[i];
 369     _ogg_free(codes);
 370
 371     c->valuelist=_book_unquantize(s,n,sortindex);
 372     c->dec_index=_ogg_malloc(n*sizeof(*c->dec_index));
 373
 374     for(n=0,i=0;i<s->entries;i++)
 375       if(s->lengthlist[i]>0)
 376         c->dec_index[sortindex[n++]]=i;
 377
 378     c->dec_codelengths=_ogg_malloc(n*sizeof(*c->dec_codelengths));
 379     c->dec_maxlength=0;
 380     for(n=0,i=0;i<s->entries;i++)
 381       if(s->lengthlist[i]>0){
 382         c->dec_codelengths[sortindex[n++]]=s->lengthlist[i];
 383         if(s->lengthlist[i]>c->dec_maxlength)
 384           c->dec_maxlength=s->lengthlist[i];
 385       }
 386
 387     if(n==1 && c->dec_maxlength==1){
 388       /* special case the 'single entry codebook' with a single bit
 389        fastpath table (that always returns entry 0 )in order to use
 390        unmodified decode paths. */
 391       c->dec_firsttablen=1;
 392       c->dec_firsttable=_ogg_calloc(2,sizeof(*c->dec_firsttable));
 393       c->dec_firsttable[0]=c->dec_firsttable[1]=1;
 394
 395     }else{
 396       c->dec_firsttablen=ov_ilog(c->used_entries)-4; /* this is magic */
 397       if(c->dec_firsttablen<5)c->dec_firsttablen=5;
 398       if(c->dec_firsttablen>8)c->dec_firsttablen=8;
 399
 400       tabn=1<<c->dec_firsttablen;
 401       c->dec_firsttable=_ogg_calloc(tabn,sizeof(*c->dec_firsttable));
 402
 403       for(i=0;i<n;i++){
 404         if(c->dec_codelengths[i]<=c->dec_firsttablen){
 405           ogg_uint32_t orig=bitreverse(c->codelist[i]);
 406           for(j=0;j<(1<<(c->dec_firsttablen-c->dec_codelengths[i]));j++)
 407             c->dec_firsttable[orig|(j<<c->dec_codelengths[i])]=i+1;
 408         }
 409       }
 410
 411       /* now fill in 'unused' entries in the firsttable with hi/lo search
 412          hints for the non-direct-hits */
 413       {
 414         ogg_uint32_t mask=0xfffffffeUL<<(31-c->dec_firsttablen);
 415         long lo=0,hi=0;
 416
 417         for(i=0;i<tabn;i++){
 418           ogg_uint32_t word=i<<(32-c->dec_firsttablen);
 419           if(c->dec_firsttable[bitreverse(word)]==0){
 420             while((lo+1)<n && c->codelist[lo+1]<=word)lo++;
 421             while(    hi<n && word>=(c->codelist[hi]&mask))hi++;
 422
 423             /* we only actually have 15 bits per hint to play with here.
 424                In order to overflow gracefully (nothing breaks, efficiency
 425                just drops), encode as the difference from the extremes. */
 426             {
 427               unsigned long loval=lo;
 428               unsigned long hival=n-hi;
 429
 430               if(loval>0x7fff)loval=0x7fff;
 431               if(hival>0x7fff)hival=0x7fff;
 432               c->dec_firsttable[bitreverse(word)]=
 433                 0x80000000UL | (loval<<15) | hival;
 434             }
 435           }
 436         }
 437       }
 438     }
 439   }
 440
 441   return(0);
 442  err_out:
 443   vorbis_book_clear(c);
 444   return(-1);
 445 }
 446
 447 long vorbis_book_codeword(codebook *book,int entry){
 448   if(book->c) /* only use with encode; decode optimizations are
 449                  allowed to break this */
 450     return book->codelist[entry];
 451   return -1;
 452 }
 453
 454 long vorbis_book_codelen(codebook *book,int entry){
 455   if(book->c) /* only use with encode; decode optimizations are
 456                  allowed to break this */
 457     return book->c->lengthlist[entry];
 458   return -1;
 459 }
 460
 461 #ifdef _V_SELFTEST
 462
 463 /* Unit tests of the dequantizer; this stuff will be OK
 464    cross-platform, I simply want to be sure that special mapping cases
 465    actually work properly; a bug could go unnoticed for a while */
 466
 467 #include <stdio.h>
 468
 469 /* cases:
 470
 471    no mapping
 472    full, explicit mapping
 473    algorithmic mapping
 474
 475    nonsequential
 476    sequential
 477 */
 478
 479 static long full_quantlist1[]={0,1,2,3,    4,5,6,7, 8,3,6,1};
 480 static long partial_quantlist1[]={0,7,2};
 481
 482 /* no mapping */
 483 static_codebook test1={
 484   4,16,
 485   NULL,
 486   0,
 487   0,0,0,0,
 488   NULL,
 489   0
 490 };
 491 static float *test1_result=NULL;
 492
 493 /* linear, full mapping, nonsequential */
 494 static_codebook test2={
 495   4,3,
 496   NULL,
 497   2,
 498   -533200896,1611661312,4,0,
 499   full_quantlist1,
 500   0
 501 };
 502 static float test2_result[]={-3,-2,-1,0, 1,2,3,4, 5,0,3,-2};
 503
 504 /* linear, full mapping, sequential */
 505 static_codebook test3={
 506   4,3,
 507   NULL,
 508   2,
 509   -533200896,1611661312,4,1,
 510   full_quantlist1,
 511   0
 512 };
 513 static float test3_result[]={-3,-5,-6,-6, 1,3,6,10, 5,5,8,6};
 514
 515 /* linear, algorithmic mapping, nonsequential */
 516 static_codebook test4={
 517   3,27,
 518   NULL,
 519   1,
 520   -533200896,1611661312,4,0,
 521   partial_quantlist1,
 522   0
 523 };
 524 static float test4_result[]={-3,-3,-3, 4,-3,-3, -1,-3,-3,
 525                               -3, 4,-3, 4, 4,-3, -1, 4,-3,
 526                               -3,-1,-3, 4,-1,-3, -1,-1,-3,
 527                               -3,-3, 4, 4,-3, 4, -1,-3, 4,
 528                               -3, 4, 4, 4, 4, 4, -1, 4, 4,
 529                               -3,-1, 4, 4,-1, 4, -1,-1, 4,
 530                               -3,-3,-1, 4,-3,-1, -1,-3,-1,
 531                               -3, 4,-1, 4, 4,-1, -1, 4,-1,
 532                               -3,-1,-1, 4,-1,-1, -1,-1,-1};
 533
 534 /* linear, algorithmic mapping, sequential */
 535 static_codebook test5={
 536   3,27,
 537   NULL,
 538   1,
 539   -533200896,1611661312,4,1,
 540   partial_quantlist1,
 541   0
 542 };
 543 static float test5_result[]={-3,-6,-9, 4, 1,-2, -1,-4,-7,
 544                               -3, 1,-2, 4, 8, 5, -1, 3, 0,
 545                               -3,-4,-7, 4, 3, 0, -1,-2,-5,
 546                               -3,-6,-2, 4, 1, 5, -1,-4, 0,
 547                               -3, 1, 5, 4, 8,12, -1, 3, 7,
 548                               -3,-4, 0, 4, 3, 7, -1,-2, 2,
 549                               -3,-6,-7, 4, 1, 0, -1,-4,-5,
 550                               -3, 1, 0, 4, 8, 7, -1, 3, 2,
 551                               -3,-4,-5, 4, 3, 2, -1,-2,-3};
 552
 553 void run_test(static_codebook *b,float *comp){
 554   float *out=_book_unquantize(b,b->entries,NULL);
 555   int i;
 556
 557   if(comp){
 558     if(!out){
 559       fprintf(stderr,"_book_unquantize incorrectly returned NULL\n");
 560       exit(1);
 561     }
 562
 563     for(i=0;i<b->entries*b->dim;i++)
 564       if(fabs(out[i]-comp[i])>.0001){
 565         fprintf(stderr,"disagreement in unquantized and reference data:\n"
 566                 "position %d, %g != %g\n",i,out[i],comp[i]);
 567         exit(1);
 568       }
 569
 570   }else{
 571     if(out){
 572       fprintf(stderr,"_book_unquantize returned a value array: \n"
 573               " correct result should have been NULL\n");
 574       exit(1);
 575     }
 576   }
 577 }
 578
 579 int main(){
 580   /* run the nine dequant tests, and compare to the hand-rolled results */
 581   fprintf(stderr,"Dequant test 1... ");
 582   run_test(&test1,test1_result);
 583   fprintf(stderr,"OK\nDequant test 2... ");
 584   run_test(&test2,test2_result);
 585   fprintf(stderr,"OK\nDequant test 3... ");
 586   run_test(&test3,test3_result);
 587   fprintf(stderr,"OK\nDequant test 4... ");
 588   run_test(&test4,test4_result);
 589   fprintf(stderr,"OK\nDequant test 5... ");
 590   run_test(&test5,test5_result);
 591   fprintf(stderr,"OK\n\n");
 592
 593   return(0);
 594 }
 595
 596 #endif