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