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