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