sys/src/cmd/audio/libvorbis/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;
 356
 357     if(codes==NULL)goto err_out;
 358     codep=malloc(sizeof(*codep)*n);
 359
 360     for(i=0;i<n;i++){
 361       codes[i]=bitreverse(codes[i]);
 362       codep[i]=codes+i;
 363     }
 364
 365     qsort(codep,n,sizeof(*codep),sort32a);
 366
 367     sortindex=malloc(n*sizeof(*sortindex));
 368     c->codelist=_ogg_malloc(n*sizeof(*c->codelist));
 369     /* the index is a reverse index */
 370     for(i=0;i<n;i++){
 371       int position=codep[i]-codes;
 372       sortindex[position]=i;
 373     }
 374
 375     for(i=0;i<n;i++)
 376       c->codelist[sortindex[i]]=codes[i];
 377     _ogg_free(codes);
 378
 379     c->valuelist=_book_unquantize(s,n,sortindex);
 380     c->dec_index=_ogg_malloc(n*sizeof(*c->dec_index));
 381
 382     for(n=0,i=0;i<s->entries;i++)
 383       if(s->lengthlist[i]>0)
 384         c->dec_index[sortindex[n++]]=i;
 385
 386     c->dec_codelengths=_ogg_malloc(n*sizeof(*c->dec_codelengths));
 387     c->dec_maxlength=0;
 388     for(n=0,i=0;i<s->entries;i++)
 389       if(s->lengthlist[i]>0){
 390         c->dec_codelengths[sortindex[n++]]=s->lengthlist[i];
 391         if(s->lengthlist[i]>c->dec_maxlength)
 392           c->dec_maxlength=s->lengthlist[i];
 393       }
 394
 395     if(n==1 && c->dec_maxlength==1){
 396       /* special case the 'single entry codebook' with a single bit
 397        fastpath table (that always returns entry 0 )in order to use
 398        unmodified decode paths. */
 399       c->dec_firsttablen=1;
 400       c->dec_firsttable=_ogg_calloc(2,sizeof(*c->dec_firsttable));
 401       c->dec_firsttable[0]=c->dec_firsttable[1]=1;
 402
 403     }else{
 404       c->dec_firsttablen=ov_ilog(c->used_entries)-4; /* this is magic */
 405       if(c->dec_firsttablen<5)c->dec_firsttablen=5;
 406       if(c->dec_firsttablen>8)c->dec_firsttablen=8;
 407
 408       tabn=1<<c->dec_firsttablen;
 409       c->dec_firsttable=_ogg_calloc(tabn,sizeof(*c->dec_firsttable));
 410
 411       for(i=0;i<n;i++){
 412         if(c->dec_codelengths[i]<=c->dec_firsttablen){
 413           ogg_uint32_t orig=bitreverse(c->codelist[i]);
 414           for(j=0;j<(1<<(c->dec_firsttablen-c->dec_codelengths[i]));j++)
 415             c->dec_firsttable[orig|(j<<c->dec_codelengths[i])]=i+1;
 416         }
 417       }
 418
 419       /* now fill in 'unused' entries in the firsttable with hi/lo search
 420          hints for the non-direct-hits */
 421       {
 422         ogg_uint32_t mask=0xfffffffeUL<<(31-c->dec_firsttablen);
 423         long lo=0,hi=0;
 424
 425         for(i=0;i<tabn;i++){
 426           ogg_uint32_t word=i<<(32-c->dec_firsttablen);
 427           if(c->dec_firsttable[bitreverse(word)]==0){
 428             while((lo+1)<n && c->codelist[lo+1]<=word)lo++;
 429             while(    hi<n && word>=(c->codelist[hi]&mask))hi++;
 430
 431             /* we only actually have 15 bits per hint to play with here.
 432                In order to overflow gracefully (nothing breaks, efficiency
 433                just drops), encode as the difference from the extremes. */
 434             {
 435               unsigned long loval=lo;
 436               unsigned long hival=n-hi;
 437
 438               if(loval>0x7fff)loval=0x7fff;
 439               if(hival>0x7fff)hival=0x7fff;
 440               c->dec_firsttable[bitreverse(word)]=
 441                 0x80000000UL | (loval<<15) | hival;
 442             }
 443           }
 444         }
 445       }
 446     }
 447     free(sortindex);
 448     free(codep);
 449   }
 450
 451   return(0);
 452  err_out:
 453   vorbis_book_clear(c);
 454   return(-1);
 455 }
 456
 457 long vorbis_book_codeword(codebook *book,int entry){
 458   if(book->c) /* only use with encode; decode optimizations are
 459                  allowed to break this */
 460     return book->codelist[entry];
 461   return -1;
 462 }
 463
 464 long vorbis_book_codelen(codebook *book,int entry){
 465   if(book->c) /* only use with encode; decode optimizations are
 466                  allowed to break this */
 467     return book->c->lengthlist[entry];
 468   return -1;
 469 }
 470
 471 #ifdef _V_SELFTEST
 472
 473 /* Unit tests of the dequantizer; this stuff will be OK
 474    cross-platform, I simply want to be sure that special mapping cases
 475    actually work properly; a bug could go unnoticed for a while */
 476
 477 #include <stdio.h>
 478
 479 /* cases:
 480
 481    no mapping
 482    full, explicit mapping
 483    algorithmic mapping
 484
 485    nonsequential
 486    sequential
 487 */
 488
 489 static long full_quantlist1[]={0,1,2,3,    4,5,6,7, 8,3,6,1};
 490 static long partial_quantlist1[]={0,7,2};
 491
 492 /* no mapping */
 493 static_codebook test1={
 494   4,16,
 495   NULL,
 496   0,
 497   0,0,0,0,
 498   NULL,
 499   0
 500 };
 501 static float *test1_result=NULL;
 502
 503 /* linear, full mapping, nonsequential */
 504 static_codebook test2={
 505   4,3,
 506   NULL,
 507   2,
 508   -533200896,1611661312,4,0,
 509   full_quantlist1,
 510   0
 511 };
 512 static float test2_result[]={-3,-2,-1,0, 1,2,3,4, 5,0,3,-2};
 513
 514 /* linear, full mapping, sequential */
 515 static_codebook test3={
 516   4,3,
 517   NULL,
 518   2,
 519   -533200896,1611661312,4,1,
 520   full_quantlist1,
 521   0
 522 };
 523 static float test3_result[]={-3,-5,-6,-6, 1,3,6,10, 5,5,8,6};
 524
 525 /* linear, algorithmic mapping, nonsequential */
 526 static_codebook test4={
 527   3,27,
 528   NULL,
 529   1,
 530   -533200896,1611661312,4,0,
 531   partial_quantlist1,
 532   0
 533 };
 534 static float test4_result[]={-3,-3,-3, 4,-3,-3, -1,-3,-3,
 535                               -3, 4,-3, 4, 4,-3, -1, 4,-3,
 536                               -3,-1,-3, 4,-1,-3, -1,-1,-3,
 537                               -3,-3, 4, 4,-3, 4, -1,-3, 4,
 538                               -3, 4, 4, 4, 4, 4, -1, 4, 4,
 539                               -3,-1, 4, 4,-1, 4, -1,-1, 4,
 540                               -3,-3,-1, 4,-3,-1, -1,-3,-1,
 541                               -3, 4,-1, 4, 4,-1, -1, 4,-1,
 542                               -3,-1,-1, 4,-1,-1, -1,-1,-1};
 543
 544 /* linear, algorithmic mapping, sequential */
 545 static_codebook test5={
 546   3,27,
 547   NULL,
 548   1,
 549   -533200896,1611661312,4,1,
 550   partial_quantlist1,
 551   0
 552 };
 553 static float test5_result[]={-3,-6,-9, 4, 1,-2, -1,-4,-7,
 554                               -3, 1,-2, 4, 8, 5, -1, 3, 0,
 555                               -3,-4,-7, 4, 3, 0, -1,-2,-5,
 556                               -3,-6,-2, 4, 1, 5, -1,-4, 0,
 557                               -3, 1, 5, 4, 8,12, -1, 3, 7,
 558                               -3,-4, 0, 4, 3, 7, -1,-2, 2,
 559                               -3,-6,-7, 4, 1, 0, -1,-4,-5,
 560                               -3, 1, 0, 4, 8, 7, -1, 3, 2,
 561                               -3,-4,-5, 4, 3, 2, -1,-2,-3};
 562
 563 void run_test(static_codebook *b,float *comp){
 564   float *out=_book_unquantize(b,b->entries,NULL);
 565   int i;
 566
 567   if(comp){
 568     if(!out){
 569       fprintf(stderr,"_book_unquantize incorrectly returned NULL\n");
 570       exit(1);
 571     }
 572
 573     for(i=0;i<b->entries*b->dim;i++)
 574       if(fabs(out[i]-comp[i])>.0001){
 575         fprintf(stderr,"disagreement in unquantized and reference data:\n"
 576                 "position %d, %g != %g\n",i,out[i],comp[i]);
 577         exit(1);
 578       }
 579
 580   }else{
 581     if(out){
 582       fprintf(stderr,"_book_unquantize returned a value array: \n"
 583               " correct result should have been NULL\n");
 584       exit(1);
 585     }
 586   }
 587   free(out);
 588 }
 589
 590 int main(){
 591   /* run the nine dequant tests, and compare to the hand-rolled results */
 592   fprintf(stderr,"Dequant test 1... ");
 593   run_test(&test1,test1_result);
 594   fprintf(stderr,"OK\nDequant test 2... ");
 595   run_test(&test2,test2_result);
 596   fprintf(stderr,"OK\nDequant test 3... ");
 597   run_test(&test3,test3_result);
 598   fprintf(stderr,"OK\nDequant test 4... ");
 599   run_test(&test4,test4_result);
 600   fprintf(stderr,"OK\nDequant test 5... ");
 601   run_test(&test5,test5_result);
 602   fprintf(stderr,"OK\n\n");
 603
 604   return(0);
 605 }
 606
 607 #endif