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
  15  ********************************************************************/
  16
  17 #include <stdlib.h>
  18 #include <limits.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;
 162   if(b->entries<1){
 163     return(0);
 164   }
 165   vals=floor(pow((float)b->entries,1.f/b->dim));
 166
 167   /* the above *should* be reliable, but we'll not assume that FP is
 168      ever reliable when bitstream sync is at stake; verify via integer
 169      means that vals really is the greatest value of dim for which
 170      vals^b->bim <= b->entries */
 171   /* treat the above as an initial guess */
 172   if(vals<1){
 173     vals=1;
 174   }
 175   while(1){
 176     long acc=1;
 177     long acc1=1;
 178     int i;
 179     for(i=0;i<b->dim;i++){
 180       if(b->entries/vals<acc)break;
 181       acc*=vals;
 182       if(LONG_MAX/(vals+1)<acc1)acc1=LONG_MAX;
 183       else acc1*=vals+1;
 184     }
 185     if(i>=b->dim && acc<=b->entries && acc1>b->entries){
 186       return(vals);
 187     }else{
 188       if(i<b->dim || acc>b->entries){
 189         vals--;
 190       }else{
 191         vals++;
 192       }
 193     }
 194   }
 195 }
 196
 197 /* unpack the quantized list of values for encode/decode ***********/
 198 /* we need to deal with two map types: in map type 1, the values are
 199    generated algorithmically (each column of the vector counts through
 200    the values in the quant vector). in map type 2, all the values came
 201    in in an explicit list.  Both value lists must be unpacked */
 202 float *_book_unquantize(const static_codebook *b,int n,int *sparsemap){
 203   long j,k,count=0;
 204   if(b->maptype==1 || b->maptype==2){
 205     int quantvals;
 206     float mindel=_float32_unpack(b->q_min);
 207     float delta=_float32_unpack(b->q_delta);
 208     float *r=_ogg_calloc(n*b->dim,sizeof(*r));
 209
 210     /* maptype 1 and 2 both use a quantized value vector, but
 211        different sizes */
 212     switch(b->maptype){
 213     case 1:
 214       /* most of the time, entries%dimensions == 0, but we need to be
 215          well defined.  We define that the possible vales at each
 216          scalar is values == entries/dim.  If entries%dim != 0, we'll
 217          have 'too few' values (values*dim<entries), which means that
 218          we'll have 'left over' entries; left over entries use zeroed
 219          values (and are wasted).  So don't generate codebooks like
 220          that */
 221       quantvals=_book_maptype1_quantvals(b);
 222       for(j=0;j<b->entries;j++){
 223         if((sparsemap && b->lengthlist[j]) || !sparsemap){
 224           float last=0.f;
 225           int indexdiv=1;
 226           for(k=0;k<b->dim;k++){
 227             int index= (j/indexdiv)%quantvals;
 228             float val=b->quantlist[index];
 229             val=fabs(val)*delta+mindel+last;
 230             if(b->q_sequencep)last=val;
 231             if(sparsemap)
 232               r[sparsemap[count]*b->dim+k]=val;
 233             else
 234               r[count*b->dim+k]=val;
 235             indexdiv*=quantvals;
 236           }
 237           count++;
 238         }
 239
 240       }
 241       break;
 242     case 2:
 243       for(j=0;j<b->entries;j++){
 244         if((sparsemap && b->lengthlist[j]) || !sparsemap){
 245           float last=0.f;
 246
 247           for(k=0;k<b->dim;k++){
 248             float val=b->quantlist[j*b->dim+k];
 249             val=fabs(val)*delta+mindel+last;
 250             if(b->q_sequencep)last=val;
 251             if(sparsemap)
 252               r[sparsemap[count]*b->dim+k]=val;
 253             else
 254               r[count*b->dim+k]=val;
 255           }
 256           count++;
 257         }
 258       }
 259       break;
 260     }
 261
 262     return(r);
 263   }
 264   return(NULL);
 265 }
 266
 267 void vorbis_staticbook_destroy(static_codebook *b){
 268   if(b->allocedp){
 269     if(b->quantlist)_ogg_free(b->quantlist);
 270     if(b->lengthlist)_ogg_free(b->lengthlist);
 271     memset(b,0,sizeof(*b));
 272     _ogg_free(b);
 273   } /* otherwise, it is in static memory */
 274 }
 275
 276 void vorbis_book_clear(codebook *b){
 277   /* static book is not cleared; we're likely called on the lookup and
 278      the static codebook belongs to the info struct */
 279   if(b->valuelist)_ogg_free(b->valuelist);
 280   if(b->codelist)_ogg_free(b->codelist);
 281
 282   if(b->dec_index)_ogg_free(b->dec_index);
 283   if(b->dec_codelengths)_ogg_free(b->dec_codelengths);
 284   if(b->dec_firsttable)_ogg_free(b->dec_firsttable);
 285
 286   memset(b,0,sizeof(*b));
 287 }
 288
 289 int vorbis_book_init_encode(codebook *c,const static_codebook *s){
 290
 291   memset(c,0,sizeof(*c));
 292   c->c=s;
 293   c->entries=s->entries;
 294   c->used_entries=s->entries;
 295   c->dim=s->dim;
 296   c->codelist=_make_words(s->lengthlist,s->entries,0);
 297   //c->valuelist=_book_unquantize(s,s->entries,NULL);
 298   c->quantvals=_book_maptype1_quantvals(s);
 299   c->minval=(int)rint(_float32_unpack(s->q_min));
 300   c->delta=(int)rint(_float32_unpack(s->q_delta));
 301
 302   return(0);
 303 }
 304
 305 static ogg_uint32_t bitreverse(ogg_uint32_t x){
 306   x=    ((x>>16)&0x0000ffffUL) | ((x<<16)&0xffff0000UL);
 307   x=    ((x>> 8)&0x00ff00ffUL) | ((x<< 8)&0xff00ff00UL);
 308   x=    ((x>> 4)&0x0f0f0f0fUL) | ((x<< 4)&0xf0f0f0f0UL);
 309   x=    ((x>> 2)&0x33333333UL) | ((x<< 2)&0xccccccccUL);
 310   return((x>> 1)&0x55555555UL) | ((x<< 1)&0xaaaaaaaaUL);
 311 }
 312
 313 static int sort32a(const void *a,const void *b){
 314   return ( **(ogg_uint32_t **)a>**(ogg_uint32_t **)b)-
 315     ( **(ogg_uint32_t **)a<**(ogg_uint32_t **)b);
 316 }
 317
 318 /* decode codebook arrangement is more heavily optimized than encode */
 319 int vorbis_book_init_decode(codebook *c,const static_codebook *s){
 320   int i,j,n=0,tabn;
 321   int *sortindex;
 322
 323   memset(c,0,sizeof(*c));
 324
 325   /* count actually used entries and find max length */
 326   for(i=0;i<s->entries;i++)
 327     if(s->lengthlist[i]>0)
 328       n++;
 329
 330   c->entries=s->entries;
 331   c->used_entries=n;
 332   c->dim=s->dim;
 333
 334   if(n>0){
 335     /* two different remappings go on here.
 336
 337     First, we collapse the likely sparse codebook down only to
 338     actually represented values/words.  This collapsing needs to be
 339     indexed as map-valueless books are used to encode original entry
 340     positions as integers.
 341
 342     Second, we reorder all vectors, including the entry index above,
 343     by sorted bitreversed codeword to allow treeless decode. */
 344
 345     /* perform sort */
 346     ogg_uint32_t *codes=_make_words(s->lengthlist,s->entries,c->used_entries);
 347     ogg_uint32_t **codep=alloca(sizeof(*codep)*n);
 348
 349     if(codes==NULL)goto err_out;
 350
 351     for(i=0;i<n;i++){
 352       codes[i]=bitreverse(codes[i]);
 353       codep[i]=codes+i;
 354     }
 355
 356     qsort(codep,n,sizeof(*codep),sort32a);
 357
 358     sortindex=alloca(n*sizeof(*sortindex));
 359     c->codelist=_ogg_malloc(n*sizeof(*c->codelist));
 360     /* the index is a reverse index */
 361     for(i=0;i<n;i++){
 362       int position=codep[i]-codes;
 363       sortindex[position]=i;
 364     }
 365
 366     for(i=0;i<n;i++)
 367       c->codelist[sortindex[i]]=codes[i];
 368     _ogg_free(codes);
 369
 370     c->valuelist=_book_unquantize(s,n,sortindex);
 371     c->dec_index=_ogg_malloc(n*sizeof(*c->dec_index));
 372
 373     for(n=0,i=0;i<s->entries;i++)
 374       if(s->lengthlist[i]>0)
 375         c->dec_index[sortindex[n++]]=i;
 376
 377     c->dec_codelengths=_ogg_malloc(n*sizeof(*c->dec_codelengths));
 378     c->dec_maxlength=0;
 379     for(n=0,i=0;i<s->entries;i++)
 380       if(s->lengthlist[i]>0){
 381         c->dec_codelengths[sortindex[n++]]=s->lengthlist[i];
 382         if(s->lengthlist[i]>c->dec_maxlength)
 383           c->dec_maxlength=s->lengthlist[i];
 384       }
 385
 386     if(n==1 && c->dec_maxlength==1){
 387       /* special case the 'single entry codebook' with a single bit
 388        fastpath table (that always returns entry 0 )in order to use
 389        unmodified decode paths. */
 390       c->dec_firsttablen=1;
 391       c->dec_firsttable=_ogg_calloc(2,sizeof(*c->dec_firsttable));
 392       c->dec_firsttable[0]=c->dec_firsttable[1]=1;
 393
 394     }else{
 395       c->dec_firsttablen=ov_ilog(c->used_entries)-4; /* this is magic */
 396       if(c->dec_firsttablen<5)c->dec_firsttablen=5;
 397       if(c->dec_firsttablen>8)c->dec_firsttablen=8;
 398
 399       tabn=1<<c->dec_firsttablen;
 400       c->dec_firsttable=_ogg_calloc(tabn,sizeof(*c->dec_firsttable));
 401
 402       for(i=0;i<n;i++){
 403         if(c->dec_codelengths[i]<=c->dec_firsttablen){
 404           ogg_uint32_t orig=bitreverse(c->codelist[i]);
 405           for(j=0;j<(1<<(c->dec_firsttablen-c->dec_codelengths[i]));j++)
 406             c->dec_firsttable[orig|(j<<c->dec_codelengths[i])]=i+1;
 407         }
 408       }
 409
 410       /* now fill in 'unused' entries in the firsttable with hi/lo search
 411          hints for the non-direct-hits */
 412       {
 413         ogg_uint32_t mask=0xfffffffeUL<<(31-c->dec_firsttablen);
 414         long lo=0,hi=0;
 415
 416         for(i=0;i<tabn;i++){
 417           ogg_uint32_t word=i<<(32-c->dec_firsttablen);
 418           if(c->dec_firsttable[bitreverse(word)]==0){
 419             while((lo+1)<n && c->codelist[lo+1]<=word)lo++;
 420             while(    hi<n && word>=(c->codelist[hi]&mask))hi++;
 421
 422             /* we only actually have 15 bits per hint to play with here.
 423                In order to overflow gracefully (nothing breaks, efficiency
 424                just drops), encode as the difference from the extremes. */
 425             {
 426               unsigned long loval=lo;
 427               unsigned long hival=n-hi;
 428
 429               if(loval>0x7fff)loval=0x7fff;
 430               if(hival>0x7fff)hival=0x7fff;
 431               c->dec_firsttable[bitreverse(word)]=
 432                 0x80000000UL | (loval<<15) | hival;
 433             }
 434           }
 435         }
 436       }
 437     }
 438   }
 439
 440   return(0);
 441  err_out:
 442   vorbis_book_clear(c);
 443   return(-1);
 444 }
 445
 446 long vorbis_book_codeword(codebook *book,int entry){
 447   if(book->c) /* only use with encode; decode optimizations are
 448                  allowed to break this */
 449     return book->codelist[entry];
 450   return -1;
 451 }
 452
 453 long vorbis_book_codelen(codebook *book,int entry){
 454   if(book->c) /* only use with encode; decode optimizations are
 455                  allowed to break this */
 456     return book->c->lengthlist[entry];
 457   return -1;
 458 }
 459
 460 #ifdef _V_SELFTEST
 461
 462 /* Unit tests of the dequantizer; this stuff will be OK
 463    cross-platform, I simply want to be sure that special mapping cases
 464    actually work properly; a bug could go unnoticed for a while */
 465
 466 #include <stdio.h>
 467
 468 /* cases:
 469
 470    no mapping
 471    full, explicit mapping
 472    algorithmic mapping
 473
 474    nonsequential
 475    sequential
 476 */
 477
 478 static long full_quantlist1[]={0,1,2,3,    4,5,6,7, 8,3,6,1};
 479 static long partial_quantlist1[]={0,7,2};
 480
 481 /* no mapping */
 482 static_codebook test1={
 483   4,16,
 484   NULL,
 485   0,
 486   0,0,0,0,
 487   NULL,
 488   0
 489 };
 490 static float *test1_result=NULL;
 491
 492 /* linear, full mapping, nonsequential */
 493 static_codebook test2={
 494   4,3,
 495   NULL,
 496   2,
 497   -533200896,1611661312,4,0,
 498   full_quantlist1,
 499   0
 500 };
 501 static float test2_result[]={-3,-2,-1,0, 1,2,3,4, 5,0,3,-2};
 502
 503 /* linear, full mapping, sequential */
 504 static_codebook test3={
 505   4,3,
 506   NULL,
 507   2,
 508   -533200896,1611661312,4,1,
 509   full_quantlist1,
 510   0
 511 };
 512 static float test3_result[]={-3,-5,-6,-6, 1,3,6,10, 5,5,8,6};
 513
 514 /* linear, algorithmic mapping, nonsequential */
 515 static_codebook test4={
 516   3,27,
 517   NULL,
 518   1,
 519   -533200896,1611661312,4,0,
 520   partial_quantlist1,
 521   0
 522 };
 523 static float test4_result[]={-3,-3,-3, 4,-3,-3, -1,-3,-3,
 524                               -3, 4,-3, 4, 4,-3, -1, 4,-3,
 525                               -3,-1,-3, 4,-1,-3, -1,-1,-3,
 526                               -3,-3, 4, 4,-3, 4, -1,-3, 4,
 527                               -3, 4, 4, 4, 4, 4, -1, 4, 4,
 528                               -3,-1, 4, 4,-1, 4, -1,-1, 4,
 529                               -3,-3,-1, 4,-3,-1, -1,-3,-1,
 530                               -3, 4,-1, 4, 4,-1, -1, 4,-1,
 531                               -3,-1,-1, 4,-1,-1, -1,-1,-1};
 532
 533 /* linear, algorithmic mapping, sequential */
 534 static_codebook test5={
 535   3,27,
 536   NULL,
 537   1,
 538   -533200896,1611661312,4,1,
 539   partial_quantlist1,
 540   0
 541 };
 542 static float test5_result[]={-3,-6,-9, 4, 1,-2, -1,-4,-7,
 543                               -3, 1,-2, 4, 8, 5, -1, 3, 0,
 544                               -3,-4,-7, 4, 3, 0, -1,-2,-5,
 545                               -3,-6,-2, 4, 1, 5, -1,-4, 0,
 546                               -3, 1, 5, 4, 8,12, -1, 3, 7,
 547                               -3,-4, 0, 4, 3, 7, -1,-2, 2,
 548                               -3,-6,-7, 4, 1, 0, -1,-4,-5,
 549                               -3, 1, 0, 4, 8, 7, -1, 3, 2,
 550                               -3,-4,-5, 4, 3, 2, -1,-2,-3};
 551
 552 void run_test(static_codebook *b,float *comp){
 553   float *out=_book_unquantize(b,b->entries,NULL);
 554   int i;
 555
 556   if(comp){
 557     if(!out){
 558       fprintf(stderr,"_book_unquantize incorrectly returned NULL\n");
 559       exit(1);
 560     }
 561
 562     for(i=0;i<b->entries*b->dim;i++)
 563       if(fabs(out[i]-comp[i])>.0001){
 564         fprintf(stderr,"disagreement in unquantized and reference data:\n"
 565                 "position %d, %g != %g\n",i,out[i],comp[i]);
 566         exit(1);
 567       }
 568
 569   }else{
 570     if(out){
 571       fprintf(stderr,"_book_unquantize returned a value array: \n"
 572               " correct result should have been NULL\n");
 573       exit(1);
 574     }
 575   }
 576 }
 577
 578 int main(){
 579   /* run the nine dequant tests, and compare to the hand-rolled results */
 580   fprintf(stderr,"Dequant test 1... ");
 581   run_test(&test1,test1_result);
 582   fprintf(stderr,"OK\nDequant test 2... ");
 583   run_test(&test2,test2_result);
 584   fprintf(stderr,"OK\nDequant test 3... ");
 585   run_test(&test3,test3_result);
 586   fprintf(stderr,"OK\nDequant test 4... ");
 587   run_test(&test4,test4_result);
 588   fprintf(stderr,"OK\nDequant test 5... ");
 589   run_test(&test5,test5_result);
 590   fprintf(stderr,"OK\n\n");
 591
 592   return(0);
 593 }
 594
 595 #endif