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