]> git.tdb.fi Git - ext/vorbisfile.git/blob - lib/sharedbook.c
Add headers to the library component so dependencies work correctly
[ext/vorbisfile.git] / 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-2015             *
9  * by the Xiph.Org Foundation https://xiph.org/                     *
10  *                                                                  *
11  ********************************************************************
12
13  function: basic shared codebook operations
14
15  ********************************************************************/
16
17 #include <stdlib.h>
18 #include <limits.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
30 int ov_ilog(ogg_uint32_t v){
31   int ret;
32   for(ret=0;v;ret++)v>>=1;
33   return ret;
34 }
35
36 /* 32 bit float (not IEEE; nonnormalized mantissa +
37    biased exponent) : neeeeeee eeemmmmm mmmmmmmm mmmmmmmm
38    Why not IEEE?  It's just not that important here. */
39
40 #define VQ_FEXP 10
41 #define VQ_FMAN 21
42 #define VQ_FEXP_BIAS 768 /* bias toward values smaller than 1. */
43
44 /* doesn't currently guard under/overflow */
45 long _float32_pack(float val){
46   int sign=0;
47   long exp;
48   long mant;
49   if(val<0){
50     sign=0x80000000;
51     val= -val;
52   }
53   exp= floor(log(val)/log(2.f)+.001); /* +epsilon */
54   mant=rint(ldexp(val,(VQ_FMAN-1)-exp));
55   exp=(exp+VQ_FEXP_BIAS)<<VQ_FMAN;
56
57   return(sign|exp|mant);
58 }
59
60 float _float32_unpack(long val){
61   double mant=val&0x1fffff;
62   int    sign=val&0x80000000;
63   long   exp =(val&0x7fe00000L)>>VQ_FMAN;
64   if(sign)mant= -mant;
65   exp=exp-(VQ_FMAN-1)-VQ_FEXP_BIAS;
66   /* clamp excessive exponent values */
67   if (exp>63){
68     exp=63;
69   }
70   if (exp<-63){
71     exp=-63;
72   }
73   return(ldexp(mant,exp));
74 }
75
76 /* given a list of word lengths, generate a list of codewords.  Works
77    for length ordered or unordered, always assigns the lowest valued
78    codewords first.  Extended to handle unused entries (length 0) */
79 ogg_uint32_t *_make_words(char *l,long n,long sparsecount){
80   long i,j,count=0;
81   ogg_uint32_t marker[33];
82   ogg_uint32_t *r=_ogg_malloc((sparsecount?sparsecount:n)*sizeof(*r));
83   memset(marker,0,sizeof(marker));
84
85   for(i=0;i<n;i++){
86     long length=l[i];
87     if(length>0){
88       ogg_uint32_t entry=marker[length];
89
90       /* when we claim a node for an entry, we also claim the nodes
91          below it (pruning off the imagined tree that may have dangled
92          from it) as well as blocking the use of any nodes directly
93          above for leaves */
94
95       /* update ourself */
96       if(length<32 && (entry>>length)){
97         /* error condition; the lengths must specify an overpopulated tree */
98         _ogg_free(r);
99         return(NULL);
100       }
101       r[count++]=entry;
102
103       /* Look to see if the next shorter marker points to the node
104          above. if so, update it and repeat.  */
105       {
106         for(j=length;j>0;j--){
107
108           if(marker[j]&1){
109             /* have to jump branches */
110             if(j==1)
111               marker[1]++;
112             else
113               marker[j]=marker[j-1]<<1;
114             break; /* invariant says next upper marker would already
115                       have been moved if it was on the same path */
116           }
117           marker[j]++;
118         }
119       }
120
121       /* prune the tree; the implicit invariant says all the longer
122          markers were dangling from our just-taken node.  Dangle them
123          from our *new* node. */
124       for(j=length+1;j<33;j++)
125         if((marker[j]>>1) == entry){
126           entry=marker[j];
127           marker[j]=marker[j-1]<<1;
128         }else
129           break;
130     }else
131       if(sparsecount==0)count++;
132   }
133
134   /* any underpopulated tree must be rejected. */
135   /* Single-entry codebooks are a retconned extension to the spec.
136      They have a single codeword '0' of length 1 that results in an
137      underpopulated tree.  Shield that case from the underformed tree check. */
138   if(!(count==1 && marker[2]==2)){
139     for(i=1;i<33;i++)
140       if(marker[i] & (0xffffffffUL>>(32-i))){
141         _ogg_free(r);
142         return(NULL);
143       }
144   }
145
146   /* bitreverse the words because our bitwise packer/unpacker is LSb
147      endian */
148   for(i=0,count=0;i<n;i++){
149     ogg_uint32_t temp=0;
150     for(j=0;j<l[i];j++){
151       temp<<=1;
152       temp|=(r[count]>>j)&1;
153     }
154
155     if(sparsecount){
156       if(l[i])
157         r[count++]=temp;
158     }else
159       r[count++]=temp;
160   }
161
162   return(r);
163 }
164
165 /* there might be a straightforward one-line way to do the below
166    that's portable and totally safe against roundoff, but I haven't
167    thought of it.  Therefore, we opt on the side of caution */
168 long _book_maptype1_quantvals(const static_codebook *b){
169   long vals;
170   if(b->entries<1){
171     return(0);
172   }
173   vals=floor(pow((float)b->entries,1.f/b->dim));
174
175   /* the above *should* be reliable, but we'll not assume that FP is
176      ever reliable when bitstream sync is at stake; verify via integer
177      means that vals really is the greatest value of dim for which
178      vals^b->bim <= b->entries */
179   /* treat the above as an initial guess */
180   if(vals<1){
181     vals=1;
182   }
183   while(1){
184     long acc=1;
185     long acc1=1;
186     int i;
187     for(i=0;i<b->dim;i++){
188       if(b->entries/vals<acc)break;
189       acc*=vals;
190       if(LONG_MAX/(vals+1)<acc1)acc1=LONG_MAX;
191       else acc1*=vals+1;
192     }
193     if(i>=b->dim && acc<=b->entries && acc1>b->entries){
194       return(vals);
195     }else{
196       if(i<b->dim || acc>b->entries){
197         vals--;
198       }else{
199         vals++;
200       }
201     }
202   }
203 }
204
205 /* unpack the quantized list of values for encode/decode ***********/
206 /* we need to deal with two map types: in map type 1, the values are
207    generated algorithmically (each column of the vector counts through
208    the values in the quant vector). in map type 2, all the values came
209    in in an explicit list.  Both value lists must be unpacked */
210 float *_book_unquantize(const static_codebook *b,int n,int *sparsemap){
211   long j,k,count=0;
212   if(b->maptype==1 || b->maptype==2){
213     int quantvals;
214     float mindel=_float32_unpack(b->q_min);
215     float delta=_float32_unpack(b->q_delta);
216     float *r=_ogg_calloc(n*b->dim,sizeof(*r));
217
218     /* maptype 1 and 2 both use a quantized value vector, but
219        different sizes */
220     switch(b->maptype){
221     case 1:
222       /* most of the time, entries%dimensions == 0, but we need to be
223          well defined.  We define that the possible vales at each
224          scalar is values == entries/dim.  If entries%dim != 0, we'll
225          have 'too few' values (values*dim<entries), which means that
226          we'll have 'left over' entries; left over entries use zeroed
227          values (and are wasted).  So don't generate codebooks like
228          that */
229       quantvals=_book_maptype1_quantvals(b);
230       for(j=0;j<b->entries;j++){
231         if((sparsemap && b->lengthlist[j]) || !sparsemap){
232           float last=0.f;
233           int indexdiv=1;
234           for(k=0;k<b->dim;k++){
235             int index= (j/indexdiv)%quantvals;
236             float val=b->quantlist[index];
237             val=fabs(val)*delta+mindel+last;
238             if(b->q_sequencep)last=val;
239             if(sparsemap)
240               r[sparsemap[count]*b->dim+k]=val;
241             else
242               r[count*b->dim+k]=val;
243             indexdiv*=quantvals;
244           }
245           count++;
246         }
247
248       }
249       break;
250     case 2:
251       for(j=0;j<b->entries;j++){
252         if((sparsemap && b->lengthlist[j]) || !sparsemap){
253           float last=0.f;
254
255           for(k=0;k<b->dim;k++){
256             float val=b->quantlist[j*b->dim+k];
257             val=fabs(val)*delta+mindel+last;
258             if(b->q_sequencep)last=val;
259             if(sparsemap)
260               r[sparsemap[count]*b->dim+k]=val;
261             else
262               r[count*b->dim+k]=val;
263           }
264           count++;
265         }
266       }
267       break;
268     }
269
270     return(r);
271   }
272   return(NULL);
273 }
274
275 void vorbis_staticbook_destroy(static_codebook *b){
276   if(b->allocedp){
277     if(b->quantlist)_ogg_free(b->quantlist);
278     if(b->lengthlist)_ogg_free(b->lengthlist);
279     memset(b,0,sizeof(*b));
280     _ogg_free(b);
281   } /* otherwise, it is in static memory */
282 }
283
284 void vorbis_book_clear(codebook *b){
285   /* static book is not cleared; we're likely called on the lookup and
286      the static codebook belongs to the info struct */
287   if(b->valuelist)_ogg_free(b->valuelist);
288   if(b->codelist)_ogg_free(b->codelist);
289
290   if(b->dec_index)_ogg_free(b->dec_index);
291   if(b->dec_codelengths)_ogg_free(b->dec_codelengths);
292   if(b->dec_firsttable)_ogg_free(b->dec_firsttable);
293
294   memset(b,0,sizeof(*b));
295 }
296
297 int vorbis_book_init_encode(codebook *c,const static_codebook *s){
298
299   memset(c,0,sizeof(*c));
300   c->c=s;
301   c->entries=s->entries;
302   c->used_entries=s->entries;
303   c->dim=s->dim;
304   c->codelist=_make_words(s->lengthlist,s->entries,0);
305   /* c->valuelist=_book_unquantize(s,s->entries,NULL); */
306   c->quantvals=_book_maptype1_quantvals(s);
307   c->minval=(int)rint(_float32_unpack(s->q_min));
308   c->delta=(int)rint(_float32_unpack(s->q_delta));
309
310   return(0);
311 }
312
313 static ogg_uint32_t bitreverse(ogg_uint32_t x){
314   x=    ((x>>16)&0x0000ffffUL) | ((x<<16)&0xffff0000UL);
315   x=    ((x>> 8)&0x00ff00ffUL) | ((x<< 8)&0xff00ff00UL);
316   x=    ((x>> 4)&0x0f0f0f0fUL) | ((x<< 4)&0xf0f0f0f0UL);
317   x=    ((x>> 2)&0x33333333UL) | ((x<< 2)&0xccccccccUL);
318   return((x>> 1)&0x55555555UL) | ((x<< 1)&0xaaaaaaaaUL);
319 }
320
321 static int sort32a(const void *a,const void *b){
322   return ( **(ogg_uint32_t **)a>**(ogg_uint32_t **)b)-
323     ( **(ogg_uint32_t **)a<**(ogg_uint32_t **)b);
324 }
325
326 /* decode codebook arrangement is more heavily optimized than encode */
327 int vorbis_book_init_decode(codebook *c,const static_codebook *s){
328   int i,j,n=0,tabn;
329   int *sortindex;
330
331   memset(c,0,sizeof(*c));
332
333   /* count actually used entries and find max length */
334   for(i=0;i<s->entries;i++)
335     if(s->lengthlist[i]>0)
336       n++;
337
338   c->entries=s->entries;
339   c->used_entries=n;
340   c->dim=s->dim;
341
342   if(n>0){
343     /* two different remappings go on here.
344
345     First, we collapse the likely sparse codebook down only to
346     actually represented values/words.  This collapsing needs to be
347     indexed as map-valueless books are used to encode original entry
348     positions as integers.
349
350     Second, we reorder all vectors, including the entry index above,
351     by sorted bitreversed codeword to allow treeless decode. */
352
353     /* perform sort */
354     ogg_uint32_t *codes=_make_words(s->lengthlist,s->entries,c->used_entries);
355     ogg_uint32_t **codep=alloca(sizeof(*codep)*n);
356
357     if(codes==NULL)goto err_out;
358
359     for(i=0;i<n;i++){
360       codes[i]=bitreverse(codes[i]);
361       codep[i]=codes+i;
362     }
363
364     qsort(codep,n,sizeof(*codep),sort32a);
365
366     sortindex=alloca(n*sizeof(*sortindex));
367     c->codelist=_ogg_malloc(n*sizeof(*c->codelist));
368     /* the index is a reverse index */
369     for(i=0;i<n;i++){
370       int position=codep[i]-codes;
371       sortindex[position]=i;
372     }
373
374     for(i=0;i<n;i++)
375       c->codelist[sortindex[i]]=codes[i];
376     _ogg_free(codes);
377
378     c->valuelist=_book_unquantize(s,n,sortindex);
379     c->dec_index=_ogg_malloc(n*sizeof(*c->dec_index));
380
381     for(n=0,i=0;i<s->entries;i++)
382       if(s->lengthlist[i]>0)
383         c->dec_index[sortindex[n++]]=i;
384
385     c->dec_codelengths=_ogg_malloc(n*sizeof(*c->dec_codelengths));
386     c->dec_maxlength=0;
387     for(n=0,i=0;i<s->entries;i++)
388       if(s->lengthlist[i]>0){
389         c->dec_codelengths[sortindex[n++]]=s->lengthlist[i];
390         if(s->lengthlist[i]>c->dec_maxlength)
391           c->dec_maxlength=s->lengthlist[i];
392       }
393
394     if(n==1 && c->dec_maxlength==1){
395       /* special case the 'single entry codebook' with a single bit
396        fastpath table (that always returns entry 0 )in order to use
397        unmodified decode paths. */
398       c->dec_firsttablen=1;
399       c->dec_firsttable=_ogg_calloc(2,sizeof(*c->dec_firsttable));
400       c->dec_firsttable[0]=c->dec_firsttable[1]=1;
401
402     }else{
403       c->dec_firsttablen=ov_ilog(c->used_entries)-4; /* this is magic */
404       if(c->dec_firsttablen<5)c->dec_firsttablen=5;
405       if(c->dec_firsttablen>8)c->dec_firsttablen=8;
406
407       tabn=1<<c->dec_firsttablen;
408       c->dec_firsttable=_ogg_calloc(tabn,sizeof(*c->dec_firsttable));
409
410       for(i=0;i<n;i++){
411         if(c->dec_codelengths[i]<=c->dec_firsttablen){
412           ogg_uint32_t orig=bitreverse(c->codelist[i]);
413           for(j=0;j<(1<<(c->dec_firsttablen-c->dec_codelengths[i]));j++)
414             c->dec_firsttable[orig|(j<<c->dec_codelengths[i])]=i+1;
415         }
416       }
417
418       /* now fill in 'unused' entries in the firsttable with hi/lo search
419          hints for the non-direct-hits */
420       {
421         ogg_uint32_t mask=0xfffffffeUL<<(31-c->dec_firsttablen);
422         long lo=0,hi=0;
423
424         for(i=0;i<tabn;i++){
425           ogg_uint32_t word=i<<(32-c->dec_firsttablen);
426           if(c->dec_firsttable[bitreverse(word)]==0){
427             while((lo+1)<n && c->codelist[lo+1]<=word)lo++;
428             while(    hi<n && word>=(c->codelist[hi]&mask))hi++;
429
430             /* we only actually have 15 bits per hint to play with here.
431                In order to overflow gracefully (nothing breaks, efficiency
432                just drops), encode as the difference from the extremes. */
433             {
434               unsigned long loval=lo;
435               unsigned long hival=n-hi;
436
437               if(loval>0x7fff)loval=0x7fff;
438               if(hival>0x7fff)hival=0x7fff;
439               c->dec_firsttable[bitreverse(word)]=
440                 0x80000000UL | (loval<<15) | hival;
441             }
442           }
443         }
444       }
445     }
446   }
447
448   return(0);
449  err_out:
450   vorbis_book_clear(c);
451   return(-1);
452 }
453
454 long vorbis_book_codeword(codebook *book,int entry){
455   if(book->c) /* only use with encode; decode optimizations are
456                  allowed to break this */
457     return book->codelist[entry];
458   return -1;
459 }
460
461 long vorbis_book_codelen(codebook *book,int entry){
462   if(book->c) /* only use with encode; decode optimizations are
463                  allowed to break this */
464     return book->c->lengthlist[entry];
465   return -1;
466 }
467
468 #ifdef _V_SELFTEST
469
470 /* Unit tests of the dequantizer; this stuff will be OK
471    cross-platform, I simply want to be sure that special mapping cases
472    actually work properly; a bug could go unnoticed for a while */
473
474 #include <stdio.h>
475
476 /* cases:
477
478    no mapping
479    full, explicit mapping
480    algorithmic mapping
481
482    nonsequential
483    sequential
484 */
485
486 static long full_quantlist1[]={0,1,2,3,    4,5,6,7, 8,3,6,1};
487 static long partial_quantlist1[]={0,7,2};
488
489 /* no mapping */
490 static_codebook test1={
491   4,16,
492   NULL,
493   0,
494   0,0,0,0,
495   NULL,
496   0
497 };
498 static float *test1_result=NULL;
499
500 /* linear, full mapping, nonsequential */
501 static_codebook test2={
502   4,3,
503   NULL,
504   2,
505   -533200896,1611661312,4,0,
506   full_quantlist1,
507   0
508 };
509 static float test2_result[]={-3,-2,-1,0, 1,2,3,4, 5,0,3,-2};
510
511 /* linear, full mapping, sequential */
512 static_codebook test3={
513   4,3,
514   NULL,
515   2,
516   -533200896,1611661312,4,1,
517   full_quantlist1,
518   0
519 };
520 static float test3_result[]={-3,-5,-6,-6, 1,3,6,10, 5,5,8,6};
521
522 /* linear, algorithmic mapping, nonsequential */
523 static_codebook test4={
524   3,27,
525   NULL,
526   1,
527   -533200896,1611661312,4,0,
528   partial_quantlist1,
529   0
530 };
531 static float test4_result[]={-3,-3,-3, 4,-3,-3, -1,-3,-3,
532                               -3, 4,-3, 4, 4,-3, -1, 4,-3,
533                               -3,-1,-3, 4,-1,-3, -1,-1,-3,
534                               -3,-3, 4, 4,-3, 4, -1,-3, 4,
535                               -3, 4, 4, 4, 4, 4, -1, 4, 4,
536                               -3,-1, 4, 4,-1, 4, -1,-1, 4,
537                               -3,-3,-1, 4,-3,-1, -1,-3,-1,
538                               -3, 4,-1, 4, 4,-1, -1, 4,-1,
539                               -3,-1,-1, 4,-1,-1, -1,-1,-1};
540
541 /* linear, algorithmic mapping, sequential */
542 static_codebook test5={
543   3,27,
544   NULL,
545   1,
546   -533200896,1611661312,4,1,
547   partial_quantlist1,
548   0
549 };
550 static float test5_result[]={-3,-6,-9, 4, 1,-2, -1,-4,-7,
551                               -3, 1,-2, 4, 8, 5, -1, 3, 0,
552                               -3,-4,-7, 4, 3, 0, -1,-2,-5,
553                               -3,-6,-2, 4, 1, 5, -1,-4, 0,
554                               -3, 1, 5, 4, 8,12, -1, 3, 7,
555                               -3,-4, 0, 4, 3, 7, -1,-2, 2,
556                               -3,-6,-7, 4, 1, 0, -1,-4,-5,
557                               -3, 1, 0, 4, 8, 7, -1, 3, 2,
558                               -3,-4,-5, 4, 3, 2, -1,-2,-3};
559
560 void run_test(static_codebook *b,float *comp){
561   float *out=_book_unquantize(b,b->entries,NULL);
562   int i;
563
564   if(comp){
565     if(!out){
566       fprintf(stderr,"_book_unquantize incorrectly returned NULL\n");
567       exit(1);
568     }
569
570     for(i=0;i<b->entries*b->dim;i++)
571       if(fabs(out[i]-comp[i])>.0001){
572         fprintf(stderr,"disagreement in unquantized and reference data:\n"
573                 "position %d, %g != %g\n",i,out[i],comp[i]);
574         exit(1);
575       }
576
577   }else{
578     if(out){
579       fprintf(stderr,"_book_unquantize returned a value array: \n"
580               " correct result should have been NULL\n");
581       exit(1);
582     }
583   }
584   free(out);
585 }
586
587 int main(){
588   /* run the nine dequant tests, and compare to the hand-rolled results */
589   fprintf(stderr,"Dequant test 1... ");
590   run_test(&test1,test1_result);
591   fprintf(stderr,"OK\nDequant test 2... ");
592   run_test(&test2,test2_result);
593   fprintf(stderr,"OK\nDequant test 3... ");
594   run_test(&test3,test3_result);
595   fprintf(stderr,"OK\nDequant test 4... ");
596   run_test(&test4,test4_result);
597   fprintf(stderr,"OK\nDequant test 5... ");
598   run_test(&test5,test5_result);
599   fprintf(stderr,"OK\n\n");
600
601   return(0);
602 }
603
604 #endif