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