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 17030 2010-03-25 06:52:55Z 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)+.001); //+epsilon
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_destroy(static_codebook *b){
262  if(b->allocedp){
263    if(b->quantlist)_ogg_free(b->quantlist);
264    if(b->lengthlist)_ogg_free(b->lengthlist);
265    memset(b,0,sizeof(*b));
266    _ogg_free(b);
267  } /* otherwise, it is in static memory */
268}
269
270void vorbis_book_clear(codebook *b){
271  /* static book is not cleared; we're likely called on the lookup and
272     the static codebook belongs to the info struct */
273  if(b->valuelist)_ogg_free(b->valuelist);
274  if(b->codelist)_ogg_free(b->codelist);
275
276  if(b->dec_index)_ogg_free(b->dec_index);
277  if(b->dec_codelengths)_ogg_free(b->dec_codelengths);
278  if(b->dec_firsttable)_ogg_free(b->dec_firsttable);
279
280  memset(b,0,sizeof(*b));
281}
282
283int vorbis_book_init_encode(codebook *c,const static_codebook *s){
284
285  memset(c,0,sizeof(*c));
286  c->c=s;
287  c->entries=s->entries;
288  c->used_entries=s->entries;
289  c->dim=s->dim;
290  c->codelist=_make_words(s->lengthlist,s->entries,0);
291  //c->valuelist=_book_unquantize(s,s->entries,NULL);
292  c->quantvals=_book_maptype1_quantvals(s);
293  c->minval=(int)rint(_float32_unpack(s->q_min));
294  c->delta=(int)rint(_float32_unpack(s->q_delta));
295
296  return(0);
297}
298
299static ogg_uint32_t bitreverse(ogg_uint32_t x){
300  x=    ((x>>16)&0x0000ffffUL) | ((x<<16)&0xffff0000UL);
301  x=    ((x>> 8)&0x00ff00ffUL) | ((x<< 8)&0xff00ff00UL);
302  x=    ((x>> 4)&0x0f0f0f0fUL) | ((x<< 4)&0xf0f0f0f0UL);
303  x=    ((x>> 2)&0x33333333UL) | ((x<< 2)&0xccccccccUL);
304  return((x>> 1)&0x55555555UL) | ((x<< 1)&0xaaaaaaaaUL);
305}
306
307static int sort32a(const void *a,const void *b){
308  return ( **(ogg_uint32_t **)a>**(ogg_uint32_t **)b)-
309    ( **(ogg_uint32_t **)a<**(ogg_uint32_t **)b);
310}
311
312/* decode codebook arrangement is more heavily optimized than encode */
313int vorbis_book_init_decode(codebook *c,const static_codebook *s){
314  int i,j,n=0,tabn;
315  int *sortindex;
316  memset(c,0,sizeof(*c));
317
318  /* count actually used entries */
319  for(i=0;i<s->entries;i++)
320    if(s->lengthlist[i]>0)
321      n++;
322
323  c->entries=s->entries;
324  c->used_entries=n;
325  c->dim=s->dim;
326
327  if(n>0){
328
329    /* two different remappings go on here.
330
331    First, we collapse the likely sparse codebook down only to
332    actually represented values/words.  This collapsing needs to be
333    indexed as map-valueless books are used to encode original entry
334    positions as integers.
335
336    Second, we reorder all vectors, including the entry index above,
337    by sorted bitreversed codeword to allow treeless decode. */
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_calloc(tabn,sizeof(*c->dec_firsttable));
383    c->dec_maxlength=0;
384
385    for(i=0;i<n;i++){
386      if(c->dec_maxlength<c->dec_codelengths[i])
387        c->dec_maxlength=c->dec_codelengths[i];
388      if(c->dec_codelengths[i]<=c->dec_firsttablen){
389        ogg_uint32_t orig=bitreverse(c->codelist[i]);
390        for(j=0;j<(1<<(c->dec_firsttablen-c->dec_codelengths[i]));j++)
391          c->dec_firsttable[orig|(j<<c->dec_codelengths[i])]=i+1;
392      }
393    }
394
395    /* now fill in 'unused' entries in the firsttable with hi/lo search
396       hints for the non-direct-hits */
397    {
398      ogg_uint32_t mask=0xfffffffeUL<<(31-c->dec_firsttablen);
399      long lo=0,hi=0;
400
401      for(i=0;i<tabn;i++){
402        ogg_uint32_t word=i<<(32-c->dec_firsttablen);
403        if(c->dec_firsttable[bitreverse(word)]==0){
404          while((lo+1)<n && c->codelist[lo+1]<=word)lo++;
405          while(    hi<n && word>=(c->codelist[hi]&mask))hi++;
406
407          /* we only actually have 15 bits per hint to play with here.
408             In order to overflow gracefully (nothing breaks, efficiency
409             just drops), encode as the difference from the extremes. */
410          {
411            unsigned long loval=lo;
412            unsigned long hival=n-hi;
413
414            if(loval>0x7fff)loval=0x7fff;
415            if(hival>0x7fff)hival=0x7fff;
416            c->dec_firsttable[bitreverse(word)]=
417              0x80000000UL | (loval<<15) | hival;
418          }
419        }
420      }
421    }
422  }
423
424  return(0);
425 err_out:
426  vorbis_book_clear(c);
427  return(-1);
428}
429
430long vorbis_book_codeword(codebook *book,int entry){
431  if(book->c) /* only use with encode; decode optimizations are
432                 allowed to break this */
433    return book->codelist[entry];
434  return -1;
435}
436
437long vorbis_book_codelen(codebook *book,int entry){
438  if(book->c) /* only use with encode; decode optimizations are
439                 allowed to break this */
440    return book->c->lengthlist[entry];
441  return -1;
442}
443
444#ifdef _V_SELFTEST
445
446/* Unit tests of the dequantizer; this stuff will be OK
447   cross-platform, I simply want to be sure that special mapping cases
448   actually work properly; a bug could go unnoticed for a while */
449
450#include <stdio.h>
451
452/* cases:
453
454   no mapping
455   full, explicit mapping
456   algorithmic mapping
457
458   nonsequential
459   sequential
460*/
461
462static long full_quantlist1[]={0,1,2,3,    4,5,6,7, 8,3,6,1};
463static long partial_quantlist1[]={0,7,2};
464
465/* no mapping */
466static_codebook test1={
467  4,16,
468  NULL,
469  0,
470  0,0,0,0,
471  NULL,
472  0
473};
474static float *test1_result=NULL;
475
476/* linear, full mapping, nonsequential */
477static_codebook test2={
478  4,3,
479  NULL,
480  2,
481  -533200896,1611661312,4,0,
482  full_quantlist1,
483  0
484};
485static float test2_result[]={-3,-2,-1,0, 1,2,3,4, 5,0,3,-2};
486
487/* linear, full mapping, sequential */
488static_codebook test3={
489  4,3,
490  NULL,
491  2,
492  -533200896,1611661312,4,1,
493  full_quantlist1,
494  0
495};
496static float test3_result[]={-3,-5,-6,-6, 1,3,6,10, 5,5,8,6};
497
498/* linear, algorithmic mapping, nonsequential */
499static_codebook test4={
500  3,27,
501  NULL,
502  1,
503  -533200896,1611661312,4,0,
504  partial_quantlist1,
505  0
506};
507static float test4_result[]={-3,-3,-3, 4,-3,-3, -1,-3,-3,
508                              -3, 4,-3, 4, 4,-3, -1, 4,-3,
509                              -3,-1,-3, 4,-1,-3, -1,-1,-3,
510                              -3,-3, 4, 4,-3, 4, -1,-3, 4,
511                              -3, 4, 4, 4, 4, 4, -1, 4, 4,
512                              -3,-1, 4, 4,-1, 4, -1,-1, 4,
513                              -3,-3,-1, 4,-3,-1, -1,-3,-1,
514                              -3, 4,-1, 4, 4,-1, -1, 4,-1,
515                              -3,-1,-1, 4,-1,-1, -1,-1,-1};
516
517/* linear, algorithmic mapping, sequential */
518static_codebook test5={
519  3,27,
520  NULL,
521  1,
522  -533200896,1611661312,4,1,
523  partial_quantlist1,
524  0
525};
526static float test5_result[]={-3,-6,-9, 4, 1,-2, -1,-4,-7,
527                              -3, 1,-2, 4, 8, 5, -1, 3, 0,
528                              -3,-4,-7, 4, 3, 0, -1,-2,-5,
529                              -3,-6,-2, 4, 1, 5, -1,-4, 0,
530                              -3, 1, 5, 4, 8,12, -1, 3, 7,
531                              -3,-4, 0, 4, 3, 7, -1,-2, 2,
532                              -3,-6,-7, 4, 1, 0, -1,-4,-5,
533                              -3, 1, 0, 4, 8, 7, -1, 3, 2,
534                              -3,-4,-5, 4, 3, 2, -1,-2,-3};
535
536void run_test(static_codebook *b,float *comp){
537  float *out=_book_unquantize(b,b->entries,NULL);
538  int i;
539
540  if(comp){
541    if(!out){
542      fprintf(stderr,"_book_unquantize incorrectly returned NULL\n");
543      exit(1);
544    }
545
546    for(i=0;i<b->entries*b->dim;i++)
547      if(fabs(out[i]-comp[i])>.0001){
548        fprintf(stderr,"disagreement in unquantized and reference data:\n"
549                "position %d, %g != %g\n",i,out[i],comp[i]);
550        exit(1);
551      }
552
553  }else{
554    if(out){
555      fprintf(stderr,"_book_unquantize returned a value array: \n"
556              " correct result should have been NULL\n");
557      exit(1);
558    }
559  }
560}
561
562int main(){
563  /* run the nine dequant tests, and compare to the hand-rolled results */
564  fprintf(stderr,"Dequant test 1... ");
565  run_test(&test1,test1_result);
566  fprintf(stderr,"OK\nDequant test 2... ");
567  run_test(&test2,test2_result);
568  fprintf(stderr,"OK\nDequant test 3... ");
569  run_test(&test3,test3_result);
570  fprintf(stderr,"OK\nDequant test 4... ");
571  run_test(&test4,test4_result);
572  fprintf(stderr,"OK\nDequant test 5... ");
573  run_test(&test5,test5_result);
574  fprintf(stderr,"OK\n\n");
575
576  return(0);
577}
578
579#endif
580