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-2007             *
   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 int _ilog(unsigned int v){
  30   int ret=0;
  31   while(v){
  32     ret++;
  33     v>>=1;
  34   }
  35   return(ret);
  36 }
  37
  38 /* 32 bit float (not IEEE; nonnormalized mantissa +
  39    biased exponent) : neeeeeee eeemmmmm mmmmmmmm mmmmmmmm
  40    Why not IEEE?  It's just not that important here. */
  41
  42 #define VQ_FEXP 10
  43 #define VQ_FMAN 21
  44 #define VQ_FEXP_BIAS 768 /* bias toward values smaller than 1. */
  45
  46 /* doesn't currently guard under/overflow */
  47 long _float32_pack(float val){
  48   int sign=0;
  49   long exp;
  50   long mant;
  51   if(val<0){
  52     sign=0x80000000;
  53     val= -val;
  54   }
  55   exp= floor(log(val)/log(2.f));
  56   mant=rint(ldexp(val,(VQ_FMAN-1)-exp));
  57   exp=(exp+VQ_FEXP_BIAS)<<VQ_FMAN;
  58
  59   return(sign|exp|mant);
  60 }
  61
  62 float _float32_unpack(long val){
  63   double mant=val&0x1fffff;
  64   int    sign=val&0x80000000;
  65   long   exp =(val&0x7fe00000L)>>VQ_FMAN;
  66   if(sign)mant= -mant;
  67   return(ldexp(mant,exp-(VQ_FMAN-1)-VQ_FEXP_BIAS));
  68 }
  69
  70 /* given a list of word lengths, generate a list of codewords.  Works
  71    for length ordered or unordered, always assigns the lowest valued
  72    codewords first.  Extended to handle unused entries (length 0) */
  73 ogg_uint32_t *_make_words(long *l,long n,long sparsecount){
  74   long i,j,count=0;
  75   ogg_uint32_t marker[33];
  76   ogg_uint32_t *r=_ogg_malloc((sparsecount?sparsecount:n)*sizeof(*r));
  77   memset(marker,0,sizeof(marker));
  78
  79   for(i=0;i<n;i++){
  80     long length=l[i];
  81     if(length>0){
  82       ogg_uint32_t entry=marker[length];
  83
  84       /* when we claim a node for an entry, we also claim the nodes
  85          below it (pruning off the imagined tree that may have dangled
  86          from it) as well as blocking the use of any nodes directly
  87          above for leaves */
  88
  89       /* update ourself */
  90       if(length<32 && (entry>>length)){
  91         /* error condition; the lengths must specify an overpopulated tree */
  92         _ogg_free(r);
  93         return(NULL);
  94       }
  95       r[count++]=entry;
  96
  97       /* Look to see if the next shorter marker points to the node
  98          above. if so, update it and repeat.  */
  99       {
 100         for(j=length;j>0;j--){
 101
 102           if(marker[j]&1){
 103             /* have to jump branches */
 104             if(j==1)
 105               marker[1]++;
 106             else
 107               marker[j]=marker[j-1]<<1;
 108             break; /* invariant says next upper marker would already
 109                       have been moved if it was on the same path */
 110           }
 111           marker[j]++;
 112         }
 113       }
 114
 115       /* prune the tree; the implicit invariant says all the longer
 116          markers were dangling from our just-taken node.  Dangle them
 117          from our *new* node. */
 118       for(j=length+1;j<33;j++)
 119         if((marker[j]>>1) == entry){
 120           entry=marker[j];
 121           marker[j]=marker[j-1]<<1;
 122         }else
 123           break;
 124     }else
 125       if(sparsecount==0)count++;
 126   }
 127
 128   /* sanity check the huffman tree; an underpopulated tree must be rejected. */
 129   for(i=1;i<33;i++)
 130     if(marker[i] & (0xffffffffUL>>(32-i))){
 131       _ogg_free(r);
 132       return(NULL);
 133     }
 134
 135   /* bitreverse the words because our bitwise packer/unpacker is LSb
 136      endian */
 137   for(i=0,count=0;i<n;i++){
 138     ogg_uint32_t temp=0;
 139     for(j=0;j<l[i];j++){
 140       temp<<=1;
 141       temp|=(r[count]>>j)&1;
 142     }
 143
 144     if(sparsecount){
 145       if(l[i])
 146         r[count++]=temp;
 147     }else
 148       r[count++]=temp;
 149   }
 150
 151   return(r);
 152 }
 153
 154 /* there might be a straightforward one-line way to do the below
 155    that's portable and totally safe against roundoff, but I haven't
 156    thought of it.  Therefore, we opt on the side of caution */
 157 long _book_maptype1_quantvals(const static_codebook *b){
 158   long vals=floor(pow((float)b->entries,1.f/b->dim));
 159
 160   /* the above *should* be reliable, but we'll not assume that FP is
 161      ever reliable when bitstream sync is at stake; verify via integer
 162      means that vals really is the greatest value of dim for which
 163      vals^b->bim <= b->entries */
 164   /* treat the above as an initial guess */
 165   while(1){
 166     long acc=1;
 167     long acc1=1;
 168     int i;
 169     for(i=0;i<b->dim;i++){
 170       acc*=vals;
 171       acc1*=vals+1;
 172     }
 173     if(acc<=b->entries && acc1>b->entries){
 174       return(vals);
 175     }else{
 176       if(acc>b->entries){
 177         vals--;
 178       }else{
 179         vals++;
 180       }
 181     }
 182   }
 183 }
 184
 185 /* unpack the quantized list of values for encode/decode ***********/
 186 /* we need to deal with two map types: in map type 1, the values are
 187    generated algorithmically (each column of the vector counts through
 188    the values in the quant vector). in map type 2, all the values came
 189    in in an explicit list.  Both value lists must be unpacked */
 190 float *_book_unquantize(const static_codebook *b,int n,int *sparsemap){
 191   long j,k,count=0;
 192   if(b->maptype==1 || b->maptype==2){
 193     int quantvals;
 194     float mindel=_float32_unpack(b->q_min);
 195     float delta=_float32_unpack(b->q_delta);
 196     float *r=_ogg_calloc(n*b->dim,sizeof(*r));
 197
 198     /* maptype 1 and 2 both use a quantized value vector, but
 199        different sizes */
 200     switch(b->maptype){
 201     case 1:
 202       /* most of the time, entries%dimensions == 0, but we need to be
 203          well defined.  We define that the possible vales at each
 204          scalar is values == entries/dim.  If entries%dim != 0, we'll
 205          have 'too few' values (values*dim<entries), which means that
 206          we'll have 'left over' entries; left over entries use zeroed
 207          values (and are wasted).  So don't generate codebooks like
 208          that */
 209       quantvals=_book_maptype1_quantvals(b);
 210       for(j=0;j<b->entries;j++){
 211         if((sparsemap && b->lengthlist[j]) || !sparsemap){
 212           float last=0.f;
 213           int indexdiv=1;
 214           for(k=0;k<b->dim;k++){
 215             int index= (j/indexdiv)%quantvals;
 216             float val=b->quantlist[index];
 217             val=fabs(val)*delta+mindel+last;
 218             if(b->q_sequencep)last=val;
 219             if(sparsemap)
 220               r[sparsemap[count]*b->dim+k]=val;
 221             else
 222               r[count*b->dim+k]=val;
 223             indexdiv*=quantvals;
 224           }
 225           count++;
 226         }
 227
 228       }
 229       break;
 230     case 2:
 231       for(j=0;j<b->entries;j++){
 232         if((sparsemap && b->lengthlist[j]) || !sparsemap){
 233           float last=0.f;
 234
 235           for(k=0;k<b->dim;k++){
 236             float val=b->quantlist[j*b->dim+k];
 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           }
 244           count++;
 245         }
 246       }
 247       break;
 248     }
 249
 250     return(r);
 251   }
 252   return(NULL);
 253 }
 254
 255 void vorbis_staticbook_clear(static_codebook *b){
 256   if(b->allocedp){
 257     if(b->quantlist)_ogg_free(b->quantlist);
 258     if(b->lengthlist)_ogg_free(b->lengthlist);
 259     if(b->nearest_tree){
 260       _ogg_free(b->nearest_tree->ptr0);
 261       _ogg_free(b->nearest_tree->ptr1);
 262       _ogg_free(b->nearest_tree->p);
 263       _ogg_free(b->nearest_tree->q);
 264       memset(b->nearest_tree,0,sizeof(*b->nearest_tree));
 265       _ogg_free(b->nearest_tree);
 266     }
 267     if(b->thresh_tree){
 268       _ogg_free(b->thresh_tree->quantthresh);
 269       _ogg_free(b->thresh_tree->quantmap);
 270       memset(b->thresh_tree,0,sizeof(*b->thresh_tree));
 271       _ogg_free(b->thresh_tree);
 272     }
 273
 274     memset(b,0,sizeof(*b));
 275   }
 276 }
 277
 278 void vorbis_staticbook_destroy(static_codebook *b){
 279   if(b->allocedp){
 280     vorbis_staticbook_clear(b);
 281     _ogg_free(b);
 282   }
 283 }
 284
 285 void vorbis_book_clear(codebook *b){
 286   /* static book is not cleared; we're likely called on the lookup and
 287      the static codebook belongs to the info struct */
 288   if(b->valuelist)_ogg_free(b->valuelist);
 289   if(b->codelist)_ogg_free(b->codelist);
 290
 291   if(b->dec_index)_ogg_free(b->dec_index);
 292   if(b->dec_codelengths)_ogg_free(b->dec_codelengths);
 293   if(b->dec_firsttable)_ogg_free(b->dec_firsttable);
 294
 295   memset(b,0,sizeof(*b));
 296 }
 297
 298 int vorbis_book_init_encode(codebook *c,const static_codebook *s){
 299
 300   memset(c,0,sizeof(*c));
 301   c->c=s;
 302   c->entries=s->entries;
 303   c->used_entries=s->entries;
 304   c->dim=s->dim;
 305   c->codelist=_make_words(s->lengthlist,s->entries,0);
 306   c->valuelist=_book_unquantize(s,s->entries,NULL);
 307
 308   return(0);
 309 }
 310
 311 static ogg_uint32_t bitreverse(ogg_uint32_t x){
 312   x=    ((x>>16)&0x0000ffffUL) | ((x<<16)&0xffff0000UL);
 313   x=    ((x>> 8)&0x00ff00ffUL) | ((x<< 8)&0xff00ff00UL);
 314   x=    ((x>> 4)&0x0f0f0f0fUL) | ((x<< 4)&0xf0f0f0f0UL);
 315   x=    ((x>> 2)&0x33333333UL) | ((x<< 2)&0xccccccccUL);
 316   return((x>> 1)&0x55555555UL) | ((x<< 1)&0xaaaaaaaaUL);
 317 }
 318
 319 static int sort32a(const void *a,const void *b){
 320   return ( **(ogg_uint32_t **)a>**(ogg_uint32_t **)b)-
 321     ( **(ogg_uint32_t **)a<**(ogg_uint32_t **)b);
 322 }
 323
 324 /* decode codebook arrangement is more heavily optimized than encode */
 325 int vorbis_book_init_decode(codebook *c,const static_codebook *s){
 326   int i,j,n=0,tabn;
 327   int *sortindex;
 328   memset(c,0,sizeof(*c));
 329
 330   /* count actually used entries */
 331   for(i=0;i<s->entries;i++)
 332     if(s->lengthlist[i]>0)
 333       n++;
 334
 335   c->entries=s->entries;
 336   c->used_entries=n;
 337   c->dim=s->dim;
 338
 339   if(n>0){
 340
 341     /* two different remappings go on here.
 342
 343     First, we collapse the likely sparse codebook down only to
 344     actually represented values/words.  This collapsing needs to be
 345     indexed as map-valueless books are used to encode original entry
 346     positions as integers.
 347
 348     Second, we reorder all vectors, including the entry index above,
 349     by sorted bitreversed codeword to allow treeless decode. */
 350
 351     /* perform sort */
 352     ogg_uint32_t *codes=_make_words(s->lengthlist,s->entries,c->used_entries);
 353     ogg_uint32_t **codep=alloca(sizeof(*codep)*n);
 354
 355     if(codes==NULL)goto err_out;
 356
 357     for(i=0;i<n;i++){
 358       codes[i]=bitreverse(codes[i]);
 359       codep[i]=codes+i;
 360     }
 361
 362     qsort(codep,n,sizeof(*codep),sort32a);
 363
 364     sortindex=alloca(n*sizeof(*sortindex));
 365     c->codelist=_ogg_malloc(n*sizeof(*c->codelist));
 366     /* the index is a reverse index */
 367     for(i=0;i<n;i++){
 368       int position=codep[i]-codes;
 369       sortindex[position]=i;
 370     }
 371
 372     for(i=0;i<n;i++)
 373       c->codelist[sortindex[i]]=codes[i];
 374     _ogg_free(codes);
 375
 376
 377     c->valuelist=_book_unquantize(s,n,sortindex);
 378     c->dec_index=_ogg_malloc(n*sizeof(*c->dec_index));
 379
 380     for(n=0,i=0;i<s->entries;i++)
 381       if(s->lengthlist[i]>0)
 382         c->dec_index[sortindex[n++]]=i;
 383
 384     c->dec_codelengths=_ogg_malloc(n*sizeof(*c->dec_codelengths));
 385     for(n=0,i=0;i<s->entries;i++)
 386       if(s->lengthlist[i]>0)
 387         c->dec_codelengths[sortindex[n++]]=s->lengthlist[i];
 388
 389     c->dec_firsttablen=_ilog(c->used_entries)-4; /* this is magic */
 390     if(c->dec_firsttablen<5)c->dec_firsttablen=5;
 391     if(c->dec_firsttablen>8)c->dec_firsttablen=8;
 392
 393     tabn=1<<c->dec_firsttablen;
 394     c->dec_firsttable=_ogg_calloc(tabn,sizeof(*c->dec_firsttable));
 395     c->dec_maxlength=0;
 396
 397     for(i=0;i<n;i++){
 398       if(c->dec_maxlength<c->dec_codelengths[i])
 399         c->dec_maxlength=c->dec_codelengths[i];
 400       if(c->dec_codelengths[i]<=c->dec_firsttablen){
 401         ogg_uint32_t orig=bitreverse(c->codelist[i]);
 402         for(j=0;j<(1<<(c->dec_firsttablen-c->dec_codelengths[i]));j++)
 403           c->dec_firsttable[orig|(j<<c->dec_codelengths[i])]=i+1;
 404       }
 405     }
 406
 407     /* now fill in 'unused' entries in the firsttable with hi/lo search
 408        hints for the non-direct-hits */
 409     {
 410       ogg_uint32_t mask=0xfffffffeUL<<(31-c->dec_firsttablen);
 411       long lo=0,hi=0;
 412
 413       for(i=0;i<tabn;i++){
 414         ogg_uint32_t word=i<<(32-c->dec_firsttablen);
 415         if(c->dec_firsttable[bitreverse(word)]==0){
 416           while((lo+1)<n && c->codelist[lo+1]<=word)lo++;
 417           while(    hi<n && word>=(c->codelist[hi]&mask))hi++;
 418
 419           /* we only actually have 15 bits per hint to play with here.
 420              In order to overflow gracefully (nothing breaks, efficiency
 421              just drops), encode as the difference from the extremes. */
 422           {
 423             unsigned long loval=lo;
 424             unsigned long hival=n-hi;
 425
 426             if(loval>0x7fff)loval=0x7fff;
 427             if(hival>0x7fff)hival=0x7fff;
 428             c->dec_firsttable[bitreverse(word)]=
 429               0x80000000UL | (loval<<15) | hival;
 430           }
 431         }
 432       }
 433     }
 434   }
 435
 436   return(0);
 437  err_out:
 438   vorbis_book_clear(c);
 439   return(-1);
 440 }
 441
 442 static float _dist(int el,float *ref, float *b,int step){
 443   int i;
 444   float acc=0.f;
 445   for(i=0;i<el;i++){
 446     float val=(ref[i]-b[i*step]);
 447     acc+=val*val;
 448   }
 449   return(acc);
 450 }
 451
 452 int _best(codebook *book, float *a, int step){
 453   encode_aux_threshmatch *tt=book->c->thresh_tree;
 454
 455 #if 0
 456   encode_aux_nearestmatch *nt=book->c->nearest_tree;
 457   encode_aux_pigeonhole *pt=book->c->pigeon_tree;
 458 #endif
 459   int dim=book->dim;
 460   int k,o;
 461   /*int savebest=-1;
 462     float saverr;*/
 463
 464   /* do we have a threshhold encode hint? */
 465   if(tt){
 466     int index=0,i;
 467     /* find the quant val of each scalar */
 468     for(k=0,o=step*(dim-1);k<dim;k++,o-=step){
 469
 470       i=tt->threshvals>>1;
 471       if(a[o]<tt->quantthresh[i]){
 472
 473         for(;i>0;i--)
 474           if(a[o]>=tt->quantthresh[i-1])
 475             break;
 476
 477       }else{
 478
 479         for(i++;i<tt->threshvals-1;i++)
 480           if(a[o]<tt->quantthresh[i])break;
 481
 482       }
 483
 484       index=(index*tt->quantvals)+tt->quantmap[i];
 485     }
 486     /* regular lattices are easy :-) */
 487     if(book->c->lengthlist[index]>0) /* is this unused?  If so, we'll
 488                                         use a decision tree after all
 489                                         and fall through*/
 490       return(index);
 491   }
 492
 493 #if 0
 494   /* do we have a pigeonhole encode hint? */
 495   if(pt){
 496     const static_codebook *c=book->c;
 497     int i,besti=-1;
 498     float best=0.f;
 499     int entry=0;
 500
 501     /* dealing with sequentialness is a pain in the ass */
 502     if(c->q_sequencep){
 503       int pv;
 504       long mul=1;
 505       float qlast=0;
 506       for(k=0,o=0;k<dim;k++,o+=step){
 507         pv=(int)((a[o]-qlast-pt->min)/pt->del);
 508         if(pv<0 || pv>=pt->mapentries)break;
 509         entry+=pt->pigeonmap[pv]*mul;
 510         mul*=pt->quantvals;
 511         qlast+=pv*pt->del+pt->min;
 512       }
 513     }else{
 514       for(k=0,o=step*(dim-1);k<dim;k++,o-=step){
 515         int pv=(int)((a[o]-pt->min)/pt->del);
 516         if(pv<0 || pv>=pt->mapentries)break;
 517         entry=entry*pt->quantvals+pt->pigeonmap[pv];
 518       }
 519     }
 520
 521     /* must be within the pigeonholable range; if we quant outside (or
 522        in an entry that we define no list for), brute force it */
 523     if(k==dim && pt->fitlength[entry]){
 524       /* search the abbreviated list */
 525       long *list=pt->fitlist+pt->fitmap[entry];
 526       for(i=0;i<pt->fitlength[entry];i++){
 527         float this=_dist(dim,book->valuelist+list[i]*dim,a,step);
 528         if(besti==-1 || this<best){
 529           best=this;
 530           besti=list[i];
 531         }
 532       }
 533
 534       return(besti);
 535     }
 536   }
 537
 538   if(nt){
 539     /* optimized using the decision tree */
 540     while(1){
 541       float c=0.f;
 542       float *p=book->valuelist+nt->p[ptr];
 543       float *q=book->valuelist+nt->q[ptr];
 544
 545       for(k=0,o=0;k<dim;k++,o+=step)
 546         c+=(p[k]-q[k])*(a[o]-(p[k]+q[k])*.5);
 547
 548       if(c>0.f) /* in A */
 549         ptr= -nt->ptr0[ptr];
 550       else     /* in B */
 551         ptr= -nt->ptr1[ptr];
 552       if(ptr<=0)break;
 553     }
 554     return(-ptr);
 555   }
 556 #endif
 557
 558   /* brute force it! */
 559   {
 560     const static_codebook *c=book->c;
 561     int i,besti=-1;
 562     float best=0.f;
 563     float *e=book->valuelist;
 564     for(i=0;i<book->entries;i++){
 565       if(c->lengthlist[i]>0){
 566         float this=_dist(dim,e,a,step);
 567         if(besti==-1 || this<best){
 568           best=this;
 569           besti=i;
 570         }
 571       }
 572       e+=dim;
 573     }
 574
 575     /*if(savebest!=-1 && savebest!=besti){
 576       fprintf(stderr,"brute force/pigeonhole disagreement:\n"
 577               "original:");
 578       for(i=0;i<dim*step;i+=step)fprintf(stderr,"%g,",a[i]);
 579       fprintf(stderr,"\n"
 580               "pigeonhole (entry %d, err %g):",savebest,saverr);
 581       for(i=0;i<dim;i++)fprintf(stderr,"%g,",
 582                                 (book->valuelist+savebest*dim)[i]);
 583       fprintf(stderr,"\n"
 584               "bruteforce (entry %d, err %g):",besti,best);
 585       for(i=0;i<dim;i++)fprintf(stderr,"%g,",
 586                                 (book->valuelist+besti*dim)[i]);
 587       fprintf(stderr,"\n");
 588       }*/
 589     return(besti);
 590   }
 591 }
 592
 593 long vorbis_book_codeword(codebook *book,int entry){
 594   if(book->c) /* only use with encode; decode optimizations are
 595                  allowed to break this */
 596     return book->codelist[entry];
 597   return -1;
 598 }
 599
 600 long vorbis_book_codelen(codebook *book,int entry){
 601   if(book->c) /* only use with encode; decode optimizations are
 602                  allowed to break this */
 603     return book->c->lengthlist[entry];
 604   return -1;
 605 }
 606
 607 #ifdef _V_SELFTEST
 608
 609 /* Unit tests of the dequantizer; this stuff will be OK
 610    cross-platform, I simply want to be sure that special mapping cases
 611    actually work properly; a bug could go unnoticed for a while */
 612
 613 #include <stdio.h>
 614
 615 /* cases:
 616
 617    no mapping
 618    full, explicit mapping
 619    algorithmic mapping
 620
 621    nonsequential
 622    sequential
 623 */
 624
 625 static long full_quantlist1[]={0,1,2,3,    4,5,6,7, 8,3,6,1};
 626 static long partial_quantlist1[]={0,7,2};
 627
 628 /* no mapping */
 629 static_codebook test1={
 630   4,16,
 631   NULL,
 632   0,
 633   0,0,0,0,
 634   NULL,
 635   NULL,NULL,NULL,
 636   0
 637 };
 638 static float *test1_result=NULL;
 639
 640 /* linear, full mapping, nonsequential */
 641 static_codebook test2={
 642   4,3,
 643   NULL,
 644   2,
 645   -533200896,1611661312,4,0,
 646   full_quantlist1,
 647   NULL,NULL,NULL,
 648   0
 649 };
 650 static float test2_result[]={-3,-2,-1,0, 1,2,3,4, 5,0,3,-2};
 651
 652 /* linear, full mapping, sequential */
 653 static_codebook test3={
 654   4,3,
 655   NULL,
 656   2,
 657   -533200896,1611661312,4,1,
 658   full_quantlist1,
 659   NULL,NULL,NULL,
 660   0
 661 };
 662 static float test3_result[]={-3,-5,-6,-6, 1,3,6,10, 5,5,8,6};
 663
 664 /* linear, algorithmic mapping, nonsequential */
 665 static_codebook test4={
 666   3,27,
 667   NULL,
 668   1,
 669   -533200896,1611661312,4,0,
 670   partial_quantlist1,
 671   NULL,NULL,NULL,
 672   0
 673 };
 674 static float test4_result[]={-3,-3,-3, 4,-3,-3, -1,-3,-3,
 675                               -3, 4,-3, 4, 4,-3, -1, 4,-3,
 676                               -3,-1,-3, 4,-1,-3, -1,-1,-3,
 677                               -3,-3, 4, 4,-3, 4, -1,-3, 4,
 678                               -3, 4, 4, 4, 4, 4, -1, 4, 4,
 679                               -3,-1, 4, 4,-1, 4, -1,-1, 4,
 680                               -3,-3,-1, 4,-3,-1, -1,-3,-1,
 681                               -3, 4,-1, 4, 4,-1, -1, 4,-1,
 682                               -3,-1,-1, 4,-1,-1, -1,-1,-1};
 683
 684 /* linear, algorithmic mapping, sequential */
 685 static_codebook test5={
 686   3,27,
 687   NULL,
 688   1,
 689   -533200896,1611661312,4,1,
 690   partial_quantlist1,
 691   NULL,NULL,NULL,
 692   0
 693 };
 694 static float test5_result[]={-3,-6,-9, 4, 1,-2, -1,-4,-7,
 695                               -3, 1,-2, 4, 8, 5, -1, 3, 0,
 696                               -3,-4,-7, 4, 3, 0, -1,-2,-5,
 697                               -3,-6,-2, 4, 1, 5, -1,-4, 0,
 698                               -3, 1, 5, 4, 8,12, -1, 3, 7,
 699                               -3,-4, 0, 4, 3, 7, -1,-2, 2,
 700                               -3,-6,-7, 4, 1, 0, -1,-4,-5,
 701                               -3, 1, 0, 4, 8, 7, -1, 3, 2,
 702                               -3,-4,-5, 4, 3, 2, -1,-2,-3};
 703
 704 void run_test(static_codebook *b,float *comp){
 705   float *out=_book_unquantize(b,b->entries,NULL);
 706   int i;
 707
 708   if(comp){
 709     if(!out){
 710       fprintf(stderr,"_book_unquantize incorrectly returned NULL\n");
 711       exit(1);
 712     }
 713
 714     for(i=0;i<b->entries*b->dim;i++)
 715       if(fabs(out[i]-comp[i])>.0001){
 716         fprintf(stderr,"disagreement in unquantized and reference data:\n"
 717                 "position %d, %g != %g\n",i,out[i],comp[i]);
 718         exit(1);
 719       }
 720
 721   }else{
 722     if(out){
 723       fprintf(stderr,"_book_unquantize returned a value array: \n"
 724               " correct result should have been NULL\n");
 725       exit(1);
 726     }
 727   }
 728 }
 729
 730 int main(){
 731   /* run the nine dequant tests, and compare to the hand-rolled results */
 732   fprintf(stderr,"Dequant test 1... ");
 733   run_test(&test1,test1_result);
 734   fprintf(stderr,"OK\nDequant test 2... ");
 735   run_test(&test2,test2_result);
 736   fprintf(stderr,"OK\nDequant test 3... ");
 737   run_test(&test3,test3_result);
 738   fprintf(stderr,"OK\nDequant test 4... ");
 739   run_test(&test4,test4_result);
 740   fprintf(stderr,"OK\nDequant test 5... ");
 741   run_test(&test5,test5_result);
 742   fprintf(stderr,"OK\n\n");
 743
 744   return(0);
 745 }
 746
 747 #endif