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
30int 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 */
45long _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
60float _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) */
79ogg_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 */
168long _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 */
210float *_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
275void 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
284void 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
297int 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
313static 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
321static 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 */
327int 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
454long 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
461long 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
486static long full_quantlist1[]={0,1,2,3, 4,5,6,7, 8,3,6,1};
487static long partial_quantlist1[]={0,7,2};
488
489/* no mapping */
490static_codebook test1={
491 4,16,
492 NULL,
493 0,
494 0,0,0,0,
495 NULL,
496 0
497};
498static float *test1_result=NULL;
499
500/* linear, full mapping, nonsequential */
501static_codebook test2={
502 4,3,
503 NULL,
504 2,
505 -533200896,1611661312,4,0,
506 full_quantlist1,
507 0
508};
509static float test2_result[]={-3,-2,-1,0, 1,2,3,4, 5,0,3,-2};
510
511/* linear, full mapping, sequential */
512static_codebook test3={
513 4,3,
514 NULL,
515 2,
516 -533200896,1611661312,4,1,
517 full_quantlist1,
518 0
519};
520static float test3_result[]={-3,-5,-6,-6, 1,3,6,10, 5,5,8,6};
521
522/* linear, algorithmic mapping, nonsequential */
523static_codebook test4={
524 3,27,
525 NULL,
526 1,
527 -533200896,1611661312,4,0,
528 partial_quantlist1,
529 0
530};
531static 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 */
542static_codebook test5={
543 3,27,
544 NULL,
545 1,
546 -533200896,1611661312,4,1,
547 partial_quantlist1,
548 0
549};
550static 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
560void 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
587int 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
605