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