Added shared lib project.
[camargo/neiasound.git] / src / stb_vorbis / stb_vorbis.c
1 // Ogg Vorbis I audio decoder  -- version 0.99996
2 //
3 // Written in April 2007 by Sean Barrett, sponsored by RAD Game Tools.
4 //
5 // Placed in the public domain April 2007 by the author: no copyright is
6 // claimed, and you may use it for any purpose you like.
7 //
8 // No warranty for any purpose is expressed or implied by the author (nor
9 // by RAD Game Tools). Report bugs and send enhancements to the author.
10 //
11 // Get the latest version and other information at:
12 //     http://nothings.org/stb_vorbis/
13
14
15 // Todo:
16 //
17 //   - seeking (note you can seek yourself using the pushdata API)
18 //
19 // Limitations:
20 //
21 //   - floor 0 not supported (used in old ogg vorbis files)
22 //   - lossless sample-truncation at beginning ignored
23 //   - cannot concatenate multiple vorbis streams
24 //   - sample positions are 32-bit, limiting seekable 192Khz
25 //       files to around 6 hours (Ogg supports 64-bit)
26 // 
27 // All of these limitations may be removed in future versions.
28
29
30 //////////////////////////////////////////////////////////////////////////////
31 //
32 //  HEADER BEGINS HERE
33 //
34
35 #ifndef STB_VORBIS_INCLUDE_STB_VORBIS_H
36 #define STB_VORBIS_INCLUDE_STB_VORBIS_H
37
38 #if defined(STB_VORBIS_NO_CRT) && !defined(STB_VORBIS_NO_STDIO)
39 #define STB_VORBIS_NO_STDIO 1
40 #endif
41
42 #ifndef STB_VORBIS_NO_STDIO
43 #include <stdio.h>
44 #endif
45
46 #ifdef __cplusplus
47 extern "C" {
48 #endif
49
50 ///////////   THREAD SAFETY
51
52 // Individual stb_vorbis* handles are not thread-safe; you cannot decode from
53 // them from multiple threads at the same time. However, you can have multiple
54 // stb_vorbis* handles and decode from them independently in multiple thrads.
55
56
57 ///////////   MEMORY ALLOCATION
58
59 // normally stb_vorbis uses malloc() to allocate memory at startup,
60 // and alloca() to allocate temporary memory during a frame on the
61 // stack. (Memory consumption will depend on the amount of setup
62 // data in the file and how you set the compile flags for speed
63 // vs. size. In my test files the maximal-size usage is ~150KB.)
64 //
65 // You can modify the wrapper functions in the source (setup_malloc,
66 // setup_temp_malloc, temp_malloc) to change this behavior, or you
67 // can use a simpler allocation model: you pass in a buffer from
68 // which stb_vorbis will allocate _all_ its memory (including the
69 // temp memory). "open" may fail with a VORBIS_outofmem if you
70 // do not pass in enough data; there is no way to determine how
71 // much you do need except to succeed (at which point you can
72 // query get_info to find the exact amount required. yes I know
73 // this is lame).
74 //
75 // If you pass in a non-NULL buffer of the type below, allocation
76 // will occur from it as described above. Otherwise just pass NULL
77 // to use malloc()/alloca()
78
79 typedef struct
80 {
81    char *alloc_buffer;
82    int   alloc_buffer_length_in_bytes;
83 } stb_vorbis_alloc;
84
85
86 ///////////   FUNCTIONS USEABLE WITH ALL INPUT MODES
87
88 typedef struct stb_vorbis stb_vorbis;
89
90 typedef struct
91 {
92    unsigned int sample_rate;
93    int channels;
94
95    unsigned int setup_memory_required;
96    unsigned int setup_temp_memory_required;
97    unsigned int temp_memory_required;
98
99    int max_frame_size;
100 } stb_vorbis_info;
101
102 // get general information about the file
103 extern stb_vorbis_info stb_vorbis_get_info(stb_vorbis *f);
104
105 // get the last error detected (clears it, too)
106 extern int stb_vorbis_get_error(stb_vorbis *f);
107
108 // close an ogg vorbis file and free all memory in use
109 extern void stb_vorbis_close(stb_vorbis *f);
110
111 // this function returns the offset (in samples) from the beginning of the
112 // file that will be returned by the next decode, if it is known, or -1
113 // otherwise. after a flush_pushdata() call, this may take a while before
114 // it becomes valid again.
115 // NOT WORKING YET after a seek with PULLDATA API
116 extern int stb_vorbis_get_sample_offset(stb_vorbis *f);
117
118 // returns the current seek point within the file, or offset from the beginning
119 // of the memory buffer. In pushdata mode it returns 0.
120 extern unsigned int stb_vorbis_get_file_offset(stb_vorbis *f);
121
122 ///////////   PUSHDATA API
123
124 #ifndef STB_VORBIS_NO_PUSHDATA_API
125
126 // this API allows you to get blocks of data from any source and hand
127 // them to stb_vorbis. you have to buffer them; stb_vorbis will tell
128 // you how much it used, and you have to give it the rest next time;
129 // and stb_vorbis may not have enough data to work with and you will
130 // need to give it the same data again PLUS more. Note that the Vorbis
131 // specification does not bound the size of an individual frame.
132
133 extern stb_vorbis *stb_vorbis_open_pushdata(
134          unsigned char *datablock, int datablock_length_in_bytes,
135          int *datablock_memory_consumed_in_bytes,
136          int *error,
137          stb_vorbis_alloc *alloc_buffer);
138 // create a vorbis decoder by passing in the initial data block containing
139 //    the ogg&vorbis headers (you don't need to do parse them, just provide
140 //    the first N bytes of the file--you're told if it's not enough, see below)
141 // on success, returns an stb_vorbis *, does not set error, returns the amount of
142 //    data parsed/consumed on this call in *datablock_memory_consumed_in_bytes;
143 // on failure, returns NULL on error and sets *error, does not change *datablock_memory_consumed
144 // if returns NULL and *error is VORBIS_need_more_data, then the input block was
145 //       incomplete and you need to pass in a larger block from the start of the file
146
147 extern int stb_vorbis_decode_frame_pushdata(
148          stb_vorbis *f, unsigned char *datablock, int datablock_length_in_bytes,
149          int *channels,             // place to write number of float * buffers
150          float ***output,           // place to write float ** array of float * buffers
151          int *samples               // place to write number of output samples
152      );
153 // decode a frame of audio sample data if possible from the passed-in data block
154 //
155 // return value: number of bytes we used from datablock
156 // possible cases:
157 //     0 bytes used, 0 samples output (need more data)
158 //     N bytes used, 0 samples output (resynching the stream, keep going)
159 //     N bytes used, M samples output (one frame of data)
160 // note that after opening a file, you will ALWAYS get one N-bytes,0-sample
161 // frame, because Vorbis always "discards" the first frame.
162 //
163 // Note that on resynch, stb_vorbis will rarely consume all of the buffer,
164 // instead only datablock_length_in_bytes-3 or less. This is because it wants
165 // to avoid missing parts of a page header if they cross a datablock boundary,
166 // without writing state-machiney code to record a partial detection.
167 //
168 // The number of channels returned are stored in *channels (which can be
169 // NULL--it is always the same as the number of channels reported by
170 // get_info). *output will contain an array of float* buffers, one per
171 // channel. In other words, (*output)[0][0] contains the first sample from
172 // the first channel, and (*output)[1][0] contains the first sample from
173 // the second channel.
174
175 extern void stb_vorbis_flush_pushdata(stb_vorbis *f);
176 // inform stb_vorbis that your next datablock will not be contiguous with
177 // previous ones (e.g. you've seeked in the data); future attempts to decode
178 // frames will cause stb_vorbis to resynchronize (as noted above), and
179 // once it sees a valid Ogg page (typically 4-8KB, as large as 64KB), it
180 // will begin decoding the _next_ frame.
181 //
182 // if you want to seek using pushdata, you need to seek in your file, then
183 // call stb_vorbis_flush_pushdata(), then start calling decoding, then once
184 // decoding is returning you data, call stb_vorbis_get_sample_offset, and
185 // if you don't like the result, seek your file again and repeat.
186 #endif
187
188
189 //////////   PULLING INPUT API
190
191 #ifndef STB_VORBIS_NO_PULLDATA_API
192 // This API assumes stb_vorbis is allowed to pull data from a source--
193 // either a block of memory containing the _entire_ vorbis stream, or a
194 // FILE * that you or it create, or possibly some other reading mechanism
195 // if you go modify the source to replace the FILE * case with some kind
196 // of callback to your code. (But if you don't support seeking, you may
197 // just want to go ahead and use pushdata.)
198
199 #if !defined(STB_VORBIS_NO_STDIO) && !defined(STB_VORBIS_NO_INTEGER_CONVERSION)
200 extern int stb_vorbis_decode_filename(char *filename, int *channels, short **output);
201 #endif
202 extern int stb_vorbis_decode_memory(unsigned char *mem, int len, int *channels, short **output);
203 // decode an entire file and output the data interleaved into a malloc()ed
204 // buffer stored in *output. The return value is the number of samples
205 // decoded, or -1 if the file could not be opened or was not an ogg vorbis file.
206 // When you're done with it, just free() the pointer returned in *output.
207
208 extern stb_vorbis * stb_vorbis_open_memory(unsigned char *data, int len,
209                                   int *error, stb_vorbis_alloc *alloc_buffer);
210 // create an ogg vorbis decoder from an ogg vorbis stream in memory (note
211 // this must be the entire stream!). on failure, returns NULL and sets *error
212
213 #ifndef STB_VORBIS_NO_STDIO
214 extern stb_vorbis * stb_vorbis_open_filename(char *filename,
215                                   int *error, stb_vorbis_alloc *alloc_buffer);
216 // create an ogg vorbis decoder from a filename via fopen(). on failure,
217 // returns NULL and sets *error (possibly to VORBIS_file_open_failure).
218
219 extern stb_vorbis * stb_vorbis_open_file(FILE *f, int close_handle_on_close,
220                                   int *error, stb_vorbis_alloc *alloc_buffer);
221 // create an ogg vorbis decoder from an open FILE *, looking for a stream at
222 // the _current_ seek point (ftell). on failure, returns NULL and sets *error.
223 // note that stb_vorbis must "own" this stream; if you seek it in between
224 // calls to stb_vorbis, it will become confused. Morever, if you attempt to
225 // perform stb_vorbis_seek_*() operations on this file, it will assume it
226 // owns the _entire_ rest of the file after the start point. Use the next
227 // function, stb_vorbis_open_file_section(), to limit it.
228
229 extern stb_vorbis * stb_vorbis_open_file_section(FILE *f, int close_handle_on_close,
230                 int *error, stb_vorbis_alloc *alloc_buffer, unsigned int len);
231 // create an ogg vorbis decoder from an open FILE *, looking for a stream at
232 // the _current_ seek point (ftell); the stream will be of length 'len' bytes.
233 // on failure, returns NULL and sets *error. note that stb_vorbis must "own"
234 // this stream; if you seek it in between calls to stb_vorbis, it will become
235 // confused.
236 #endif
237
238 extern int stb_vorbis_seek_frame(stb_vorbis *f, unsigned int sample_number);
239 extern int stb_vorbis_seek(stb_vorbis *f, unsigned int sample_number);
240 // NOT WORKING YET
241 // these functions seek in the Vorbis file to (approximately) 'sample_number'.
242 // after calling seek_frame(), the next call to get_frame_*() will include
243 // the specified sample. after calling stb_vorbis_seek(), the next call to
244 // stb_vorbis_get_samples_* will start with the specified sample. If you
245 // do not need to seek to EXACTLY the target sample when using get_samples_*,
246 // you can also use seek_frame().
247
248 extern void stb_vorbis_seek_start(stb_vorbis *f);
249 // this function is equivalent to stb_vorbis_seek(f,0), but it
250 // actually works
251
252 extern unsigned int stb_vorbis_stream_length_in_samples(stb_vorbis *f);
253 extern float        stb_vorbis_stream_length_in_seconds(stb_vorbis *f);
254 // these functions return the total length of the vorbis stream
255
256 extern int stb_vorbis_get_frame_float(stb_vorbis *f, int *channels, float ***output);
257 // decode the next frame and return the number of samples. the number of
258 // channels returned are stored in *channels (which can be NULL--it is always
259 // the same as the number of channels reported by get_info). *output will
260 // contain an array of float* buffers, one per channel. These outputs will
261 // be overwritten on the next call to stb_vorbis_get_frame_*.
262 //
263 // You generally should not intermix calls to stb_vorbis_get_frame_*()
264 // and stb_vorbis_get_samples_*(), since the latter calls the former.
265
266 #ifndef STB_VORBIS_NO_INTEGER_CONVERSION
267 extern int stb_vorbis_get_frame_short_interleaved(stb_vorbis *f, int num_c, short *buffer, int num_shorts);
268 extern int stb_vorbis_get_frame_short            (stb_vorbis *f, int num_c, short **buffer, int num_samples);
269 #endif
270 // decode the next frame and return the number of samples per channel. the
271 // data is coerced to the number of channels you request according to the
272 // channel coercion rules (see below). You must pass in the size of your
273 // buffer(s) so that stb_vorbis will not overwrite the end of the buffer.
274 // The maximum buffer size needed can be gotten from get_info(); however,
275 // the Vorbis I specification implies an absolute maximum of 4096 samples
276 // per channel. Note that for interleaved data, you pass in the number of
277 // shorts (the size of your array), but the return value is the number of
278 // samples per channel, not the total number of samples.
279
280 // Channel coercion rules:
281 //    Let M be the number of channels requested, and N the number of channels present,
282 //    and Cn be the nth channel; let stereo L be the sum of all L and center channels,
283 //    and stereo R be the sum of all R and center channels (channel assignment from the
284 //    vorbis spec).
285 //        M    N       output
286 //        1    k      sum(Ck) for all k
287 //        2    *      stereo L, stereo R
288 //        k    l      k > l, the first l channels, then 0s
289 //        k    l      k <= l, the first k channels
290 //    Note that this is not _good_ surround etc. mixing at all! It's just so
291 //    you get something useful.
292
293 extern int stb_vorbis_get_samples_float_interleaved(stb_vorbis *f, int channels, float *buffer, int num_floats);
294 extern int stb_vorbis_get_samples_float(stb_vorbis *f, int channels, float **buffer, int num_samples);
295 // gets num_samples samples, not necessarily on a frame boundary--this requires
296 // buffering so you have to supply the buffers. DOES NOT APPLY THE COERCION RULES.
297 // Returns the number of samples stored per channel; it may be less than requested
298 // at the end of the file. If there are no more samples in the file, returns 0.
299
300 #ifndef STB_VORBIS_NO_INTEGER_CONVERSION
301 extern int stb_vorbis_get_samples_short_interleaved(stb_vorbis *f, int channels, short *buffer, int num_shorts);
302 extern int stb_vorbis_get_samples_short(stb_vorbis *f, int channels, short **buffer, int num_samples);
303 #endif
304 // gets num_samples samples, not necessarily on a frame boundary--this requires
305 // buffering so you have to supply the buffers. Applies the coercion rules above
306 // to produce 'channels' channels. Returns the number of samples stored per channel;
307 // it may be less than requested at the end of the file. If there are no more
308 // samples in the file, returns 0.
309
310 #endif
311
312 ////////   ERROR CODES
313
314 enum STBVorbisError
315 {
316    VORBIS__no_error,
317
318    VORBIS_need_more_data=1,             // not a real error
319
320    VORBIS_invalid_api_mixing,           // can't mix API modes
321    VORBIS_outofmem,                     // not enough memory
322    VORBIS_feature_not_supported,        // uses floor 0
323    VORBIS_too_many_channels,            // STB_VORBIS_MAX_CHANNELS is too small
324    VORBIS_file_open_failure,            // fopen() failed
325    VORBIS_seek_without_length,          // can't seek in unknown-length file
326
327    VORBIS_unexpected_eof=10,            // file is truncated?
328    VORBIS_seek_invalid,                 // seek past EOF
329
330    // decoding errors (corrupt/invalid stream) -- you probably
331    // don't care about the exact details of these
332
333    // vorbis errors:
334    VORBIS_invalid_setup=20,
335    VORBIS_invalid_stream,
336
337    // ogg errors:
338    VORBIS_missing_capture_pattern=30,
339    VORBIS_invalid_stream_structure_version,
340    VORBIS_continued_packet_flag_invalid,
341    VORBIS_incorrect_stream_serial_number,
342    VORBIS_invalid_first_page,
343    VORBIS_bad_packet_type,
344    VORBIS_cant_find_last_page,
345    VORBIS_seek_failed,
346 };
347
348
349 #ifdef __cplusplus
350 }
351 #endif
352
353 #endif // STB_VORBIS_INCLUDE_STB_VORBIS_H
354 //
355 //  HEADER ENDS HERE
356 //
357 //////////////////////////////////////////////////////////////////////////////
358
359 #ifndef STB_VORBIS_HEADER_ONLY
360
361 // global configuration settings (e.g. set these in the project/makefile),
362 // or just set them in this file at the top (although ideally the first few
363 // should be visible when the header file is compiled too, although it's not
364 // crucial)
365
366 // STB_VORBIS_NO_PUSHDATA_API
367 //     does not compile the code for the various stb_vorbis_*_pushdata()
368 //     functions
369 // #define STB_VORBIS_NO_PUSHDATA_API
370
371 // STB_VORBIS_NO_PULLDATA_API
372 //     does not compile the code for the non-pushdata APIs
373 // #define STB_VORBIS_NO_PULLDATA_API
374
375 // STB_VORBIS_NO_STDIO
376 //     does not compile the code for the APIs that use FILE *s internally
377 //     or externally (implied by STB_VORBIS_NO_PULLDATA_API)
378 // #define STB_VORBIS_NO_STDIO
379
380 // STB_VORBIS_NO_INTEGER_CONVERSION
381 //     does not compile the code for converting audio sample data from
382 //     float to integer (implied by STB_VORBIS_NO_PULLDATA_API)
383 // #define STB_VORBIS_NO_INTEGER_CONVERSION
384
385 // STB_VORBIS_NO_FAST_SCALED_FLOAT
386 //      does not use a fast float-to-int trick to accelerate float-to-int on
387 //      most platforms which requires endianness be defined correctly.
388 //#define STB_VORBIS_NO_FAST_SCALED_FLOAT
389
390
391 // STB_VORBIS_MAX_CHANNELS [number]
392 //     globally define this to the maximum number of channels you need.
393 //     The spec does not put a restriction on channels except that
394 //     the count is stored in a byte, so 255 is the hard limit.
395 //     Reducing this saves about 16 bytes per value, so using 16 saves
396 //     (255-16)*16 or around 4KB. Plus anything other memory usage
397 //     I forgot to account for. Can probably go as low as 8 (7.1 audio),
398 //     6 (5.1 audio), or 2 (stereo only).
399 #ifndef STB_VORBIS_MAX_CHANNELS
400 #define STB_VORBIS_MAX_CHANNELS    16  // enough for anyone?
401 #endif
402
403 // STB_VORBIS_PUSHDATA_CRC_COUNT [number]
404 //     after a flush_pushdata(), stb_vorbis begins scanning for the
405 //     next valid page, without backtracking. when it finds something
406 //     that looks like a page, it streams through it and verifies its
407 //     CRC32. Should that validation fail, it keeps scanning. But it's
408 //     possible that _while_ streaming through to check the CRC32 of
409 //     one candidate page, it sees another candidate page. This #define
410 //     determines how many "overlapping" candidate pages it can search
411 //     at once. Note that "real" pages are typically ~4KB to ~8KB, whereas
412 //     garbage pages could be as big as 64KB, but probably average ~16KB.
413 //     So don't hose ourselves by scanning an apparent 64KB page and
414 //     missing a ton of real ones in the interim; so minimum of 2
415 #ifndef STB_VORBIS_PUSHDATA_CRC_COUNT
416 #define STB_VORBIS_PUSHDATA_CRC_COUNT  4
417 #endif
418
419 // STB_VORBIS_FAST_HUFFMAN_LENGTH [number]
420 //     sets the log size of the huffman-acceleration table.  Maximum
421 //     supported value is 24. with larger numbers, more decodings are O(1),
422 //     but the table size is larger so worse cache missing, so you'll have
423 //     to probe (and try multiple ogg vorbis files) to find the sweet spot.
424 #ifndef STB_VORBIS_FAST_HUFFMAN_LENGTH
425 #define STB_VORBIS_FAST_HUFFMAN_LENGTH   10
426 #endif
427
428 // STB_VORBIS_FAST_BINARY_LENGTH [number]
429 //     sets the log size of the binary-search acceleration table. this
430 //     is used in similar fashion to the fast-huffman size to set initial
431 //     parameters for the binary search
432
433 // STB_VORBIS_FAST_HUFFMAN_INT
434 //     The fast huffman tables are much more efficient if they can be
435 //     stored as 16-bit results instead of 32-bit results. This restricts
436 //     the codebooks to having only 65535 possible outcomes, though.
437 //     (At least, accelerated by the huffman table.)
438 #ifndef STB_VORBIS_FAST_HUFFMAN_INT
439 #define STB_VORBIS_FAST_HUFFMAN_SHORT
440 #endif
441
442 // STB_VORBIS_NO_HUFFMAN_BINARY_SEARCH
443 //     If the 'fast huffman' search doesn't succeed, then stb_vorbis falls
444 //     back on binary searching for the correct one. This requires storing
445 //     extra tables with the huffman codes in sorted order. Defining this
446 //     symbol trades off space for speed by forcing a linear search in the
447 //     non-fast case, except for "sparse" codebooks.
448 // #define STB_VORBIS_NO_HUFFMAN_BINARY_SEARCH
449
450 // STB_VORBIS_DIVIDES_IN_RESIDUE
451 //     stb_vorbis precomputes the result of the scalar residue decoding
452 //     that would otherwise require a divide per chunk. you can trade off
453 //     space for time by defining this symbol.
454 // #define STB_VORBIS_DIVIDES_IN_RESIDUE
455
456 // STB_VORBIS_DIVIDES_IN_CODEBOOK
457 //     vorbis VQ codebooks can be encoded two ways: with every case explicitly
458 //     stored, or with all elements being chosen from a small range of values,
459 //     and all values possible in all elements. By default, stb_vorbis expands
460 //     this latter kind out to look like the former kind for ease of decoding,
461 //     because otherwise an integer divide-per-vector-element is required to
462 //     unpack the index. If you define STB_VORBIS_DIVIDES_IN_CODEBOOK, you can
463 //     trade off storage for speed.
464 //#define STB_VORBIS_DIVIDES_IN_CODEBOOK
465
466 // STB_VORBIS_CODEBOOK_SHORTS
467 //     The vorbis file format encodes VQ codebook floats as ax+b where a and
468 //     b are floating point per-codebook constants, and x is a 16-bit int.
469 //     Normally, stb_vorbis decodes them to floats rather than leaving them
470 //     as 16-bit ints and computing ax+b while decoding. This is a speed/space
471 //     tradeoff; you can save space by defining this flag.
472 #ifndef STB_VORBIS_CODEBOOK_SHORTS
473 #define STB_VORBIS_CODEBOOK_FLOATS
474 #endif
475
476 // STB_VORBIS_DIVIDE_TABLE
477 //     this replaces small integer divides in the floor decode loop with
478 //     table lookups. made less than 1% difference, so disabled by default.
479
480 // STB_VORBIS_NO_INLINE_DECODE
481 //     disables the inlining of the scalar codebook fast-huffman decode.
482 //     might save a little codespace; useful for debugging
483 // #define STB_VORBIS_NO_INLINE_DECODE
484
485 // STB_VORBIS_NO_DEFER_FLOOR
486 //     Normally we only decode the floor without synthesizing the actual
487 //     full curve. We can instead synthesize the curve immediately. This
488 //     requires more memory and is very likely slower, so I don't think
489 //     you'd ever want to do it except for debugging.
490 // #define STB_VORBIS_NO_DEFER_FLOOR
491
492
493
494
495 //////////////////////////////////////////////////////////////////////////////
496
497 #ifdef STB_VORBIS_NO_PULLDATA_API
498    #define STB_VORBIS_NO_INTEGER_CONVERSION
499    #define STB_VORBIS_NO_STDIO
500 #endif
501
502 #if defined(STB_VORBIS_NO_CRT) && !defined(STB_VORBIS_NO_STDIO)
503    #define STB_VORBIS_NO_STDIO 1
504 #endif
505
506 #ifndef STB_VORBIS_NO_INTEGER_CONVERSION
507 #ifndef STB_VORBIS_NO_FAST_SCALED_FLOAT
508
509    // only need endianness for fast-float-to-int, which we don't
510    // use for pushdata
511
512    #ifndef STB_VORBIS_BIG_ENDIAN
513      #define STB_VORBIS_ENDIAN  0
514    #else
515      #define STB_VORBIS_ENDIAN  1
516    #endif
517
518 #endif
519 #endif
520
521
522 #ifndef STB_VORBIS_NO_STDIO
523 #include <stdio.h>
524 #endif
525
526 #ifndef STB_VORBIS_NO_CRT
527 #include <stdlib.h>
528 #include <string.h>
529 #include <assert.h>
530 #include <math.h>
531 #if !(defined(__APPLE__) || defined(MACOSX) || defined(macintosh) || defined(Macintosh))
532 #include <malloc.h>
533 #endif
534 #else
535 #define NULL 0
536 #endif
537
538 #ifndef _MSC_VER
539    #if __GNUC__
540       #define __forceinline inline
541    #else
542       #define __forceinline
543    #endif
544 #endif
545
546 #if STB_VORBIS_MAX_CHANNELS > 256
547 #error "Value of STB_VORBIS_MAX_CHANNELS outside of allowed range"
548 #endif
549
550 #if STB_VORBIS_FAST_HUFFMAN_LENGTH > 24
551 #error "Value of STB_VORBIS_FAST_HUFFMAN_LENGTH outside of allowed range"
552 #endif
553
554
555 #define MAX_BLOCKSIZE_LOG  13   // from specification
556 #define MAX_BLOCKSIZE      (1 << MAX_BLOCKSIZE_LOG)
557
558
559 typedef unsigned char  uint8;
560 typedef   signed char   int8;
561 typedef unsigned short uint16;
562 typedef   signed short  int16;
563 typedef unsigned int   uint32;
564 typedef   signed int    int32;
565
566 #ifndef TRUE
567 #define TRUE 1
568 #define FALSE 0
569 #endif
570
571 #ifdef STB_VORBIS_CODEBOOK_FLOATS
572 typedef float codetype;
573 #else
574 typedef uint16 codetype;
575 #endif
576
577 // @NOTE
578 //
579 // Some arrays below are tagged "//varies", which means it's actually
580 // a variable-sized piece of data, but rather than malloc I assume it's
581 // small enough it's better to just allocate it all together with the
582 // main thing
583 //
584 // Most of the variables are specified with the smallest size I could pack
585 // them into. It might give better performance to make them all full-sized
586 // integers. It should be safe to freely rearrange the structures or change
587 // the sizes larger--nothing relies on silently truncating etc., nor the
588 // order of variables.
589
590 #define FAST_HUFFMAN_TABLE_SIZE   (1 << STB_VORBIS_FAST_HUFFMAN_LENGTH)
591 #define FAST_HUFFMAN_TABLE_MASK   (FAST_HUFFMAN_TABLE_SIZE - 1)
592
593 typedef struct
594 {
595    int dimensions, entries;
596    uint8 *codeword_lengths;
597    float  minimum_value;
598    float  delta_value;
599    uint8  value_bits;
600    uint8  lookup_type;
601    uint8  sequence_p;
602    uint8  sparse;
603    uint32 lookup_values;
604    codetype *multiplicands;
605    uint32 *codewords;
606    #ifdef STB_VORBIS_FAST_HUFFMAN_SHORT
607     int16  fast_huffman[FAST_HUFFMAN_TABLE_SIZE];
608    #else
609     int32  fast_huffman[FAST_HUFFMAN_TABLE_SIZE];
610    #endif
611    uint32 *sorted_codewords;
612    int    *sorted_values;
613    int     sorted_entries;
614 } Codebook;
615
616 typedef struct
617 {
618    uint8 order;
619    uint16 rate;
620    uint16 bark_map_size;
621    uint8 amplitude_bits;
622    uint8 amplitude_offset;
623    uint8 number_of_books;
624    uint8 book_list[16]; // varies
625 } Floor0;
626
627 typedef struct
628 {
629    uint8 partitions;
630    uint8 partition_class_list[32]; // varies
631    uint8 class_dimensions[16]; // varies
632    uint8 class_subclasses[16]; // varies
633    uint8 class_masterbooks[16]; // varies
634    int16 subclass_books[16][8]; // varies
635    uint16 Xlist[31*8+2]; // varies
636    uint8 sorted_order[31*8+2];
637    uint8 neighbors[31*8+2][2];
638    uint8 floor1_multiplier;
639    uint8 rangebits;
640    int values;
641 } Floor1;
642
643 typedef union
644 {
645    Floor0 floor0;
646    Floor1 floor1;
647 } Floor;
648
649 typedef struct
650 {
651    uint32 begin, end;
652    uint32 part_size;
653    uint8 classifications;
654    uint8 classbook;
655    uint8 **classdata;
656    int16 (*residue_books)[8];
657 } Residue;
658
659 typedef struct
660 {
661    uint8 magnitude;
662    uint8 angle;
663    uint8 mux;
664 } MappingChannel;
665
666 typedef struct
667 {
668    uint16 coupling_steps;
669    MappingChannel *chan;
670    uint8  submaps;
671    uint8  submap_floor[15]; // varies
672    uint8  submap_residue[15]; // varies
673 } Mapping;
674
675 typedef struct
676 {
677    uint8 blockflag;
678    uint8 mapping;
679    uint16 windowtype;
680    uint16 transformtype;
681 } Mode;
682
683 typedef struct
684 {
685    uint32  goal_crc;    // expected crc if match
686    int     bytes_left;  // bytes left in packet
687    uint32  crc_so_far;  // running crc
688    int     bytes_done;  // bytes processed in _current_ chunk
689    uint32  sample_loc;  // granule pos encoded in page
690 } CRCscan;
691
692 typedef struct
693 {
694    uint32 page_start, page_end;
695    uint32 after_previous_page_start;
696    uint32 first_decoded_sample;
697    uint32 last_decoded_sample;
698 } ProbedPage;
699
700 struct stb_vorbis
701 {
702   // user-accessible info
703    unsigned int sample_rate;
704    int channels;
705
706    unsigned int setup_memory_required;
707    unsigned int temp_memory_required;
708    unsigned int setup_temp_memory_required;
709
710   // input config
711 #ifndef STB_VORBIS_NO_STDIO
712    FILE *f;
713    uint32 f_start;
714    int close_on_free;
715 #endif
716
717    uint8 *stream;
718    uint8 *stream_start;
719    uint8 *stream_end;
720
721    uint32 stream_len;
722
723    uint8  push_mode;
724
725    uint32 first_audio_page_offset;
726
727    ProbedPage p_first, p_last;
728
729   // memory management
730    stb_vorbis_alloc alloc;
731    int setup_offset;
732    int temp_offset;
733
734   // run-time results
735    int eof;
736    enum STBVorbisError error;
737
738   // user-useful data
739
740   // header info
741    int blocksize[2];
742    int blocksize_0, blocksize_1;
743    int codebook_count;
744    Codebook *codebooks;
745    int floor_count;
746    uint16 floor_types[64]; // varies
747    Floor *floor_config;
748    int residue_count;
749    uint16 residue_types[64]; // varies
750    Residue *residue_config;
751    int mapping_count;
752    Mapping *mapping;
753    int mode_count;
754    Mode mode_config[64];  // varies
755
756    uint32 total_samples;
757
758   // decode buffer
759    float *channel_buffers[STB_VORBIS_MAX_CHANNELS];
760    float *outputs        [STB_VORBIS_MAX_CHANNELS];
761
762    float *previous_window[STB_VORBIS_MAX_CHANNELS];
763    int previous_length;
764
765    #ifndef STB_VORBIS_NO_DEFER_FLOOR
766    int16 *finalY[STB_VORBIS_MAX_CHANNELS];
767    #else
768    float *floor_buffers[STB_VORBIS_MAX_CHANNELS];
769    #endif
770
771    uint32 current_loc; // sample location of next frame to decode
772    int    current_loc_valid;
773
774   // per-blocksize precomputed data
775    
776    // twiddle factors
777    float *A[2],*B[2],*C[2];
778    float *window[2];
779    uint16 *bit_reverse[2];
780
781   // current page/packet/segment streaming info
782    uint32 serial; // stream serial number for verification
783    int last_page;
784    int segment_count;
785    uint8 segments[255];
786    uint8 page_flag;
787    uint8 bytes_in_seg;
788    uint8 first_decode;
789    int next_seg;
790    int last_seg;  // flag that we're on the last segment
791    int last_seg_which; // what was the segment number of the last seg?
792    uint32 acc;
793    int valid_bits;
794    int packet_bytes;
795    int end_seg_with_known_loc;
796    uint32 known_loc_for_packet;
797    int discard_samples_deferred;
798    uint32 samples_output;
799
800   // push mode scanning
801    int page_crc_tests; // only in push_mode: number of tests active; -1 if not searching
802 #ifndef STB_VORBIS_NO_PUSHDATA_API
803    CRCscan scan[STB_VORBIS_PUSHDATA_CRC_COUNT];
804 #endif
805
806   // sample-access
807    int channel_buffer_start;
808    int channel_buffer_end;
809 };
810
811 extern int my_prof(int slot);
812 //#define stb_prof my_prof
813
814 #ifndef stb_prof
815 #define stb_prof(x)  0
816 #endif
817
818 #if defined(STB_VORBIS_NO_PUSHDATA_API)
819    #define IS_PUSH_MODE(f)   FALSE
820 #elif defined(STB_VORBIS_NO_PULLDATA_API)
821    #define IS_PUSH_MODE(f)   TRUE
822 #else
823    #define IS_PUSH_MODE(f)   ((f)->push_mode)
824 #endif
825
826 typedef struct stb_vorbis vorb;
827
828 static int error(vorb *f, enum STBVorbisError e)
829 {
830    f->error = e;
831    if (!f->eof && e != VORBIS_need_more_data) {
832       f->error=e; // breakpoint for debugging
833    }
834    return 0;
835 }
836
837
838 // these functions are used for allocating temporary memory
839 // while decoding. if you can afford the stack space, use
840 // alloca(); otherwise, provide a temp buffer and it will
841 // allocate out of those.
842
843 #define array_size_required(count,size)  (count*(sizeof(void *)+(size)))
844
845 #define temp_alloc(f,size)              (f->alloc.alloc_buffer ? setup_temp_malloc(f,size) : alloca(size))
846 #ifdef dealloca
847 #define temp_free(f,p)                  (f->alloc.alloc_buffer ? 0 : dealloca(size))
848 #else
849 #define temp_free(f,p)                  0
850 #endif
851 #define temp_alloc_save(f)              ((f)->temp_offset)
852 #define temp_alloc_restore(f,p)         ((f)->temp_offset = (p))
853
854 #define temp_block_array(f,count,size)  make_block_array(temp_alloc(f,array_size_required(count,size)), count, size)
855
856 // given a sufficiently large block of memory, make an array of pointers to subblocks of it
857 static void *make_block_array(void *mem, int count, int size)
858 {
859    int i;
860    void ** p = (void **) mem;
861    char *q = (char *) (p + count);
862    for (i=0; i < count; ++i) {
863       p[i] = q;
864       q += size;
865    }
866    return p;
867 }
868
869 static void *setup_malloc(vorb *f, int sz)
870 {
871    sz = (sz+3) & ~3;
872    f->setup_memory_required += sz;
873    if (f->alloc.alloc_buffer) {
874       void *p = (char *) f->alloc.alloc_buffer + f->setup_offset;
875       if (f->setup_offset + sz > f->temp_offset) return NULL;
876       f->setup_offset += sz;
877       return p;
878    }
879    return sz ? malloc(sz) : NULL;
880 }
881
882 static void setup_free(vorb *f, void *p)
883 {
884    if (f->alloc.alloc_buffer) return; // do nothing; setup mem is not a stack
885    free(p);
886 }
887
888 static void *setup_temp_malloc(vorb *f, int sz)
889 {
890    sz = (sz+3) & ~3;
891    if (f->alloc.alloc_buffer) {
892       if (f->temp_offset - sz < f->setup_offset) return NULL;
893       f->temp_offset -= sz;
894       return (char *) f->alloc.alloc_buffer + f->temp_offset;
895    }
896    return malloc(sz);
897 }
898
899 static void setup_temp_free(vorb *f, void *p, size_t sz)
900 {
901    if (f->alloc.alloc_buffer) {
902       f->temp_offset += (sz+3)&~3;
903       return;
904    }
905    free(p);
906 }
907
908 #define CRC32_POLY    0x04c11db7   // from spec
909
910 static uint32 crc_table[256];
911 static void crc32_init(void)
912 {
913    int i,j;
914    uint32 s;
915    for(i=0; i < 256; i++) {
916       for (s=i<<24, j=0; j < 8; ++j)
917          s = (s << 1) ^ (s >= (1<<31) ? CRC32_POLY : 0);
918       crc_table[i] = s;
919    }
920 }
921
922 static __forceinline uint32 crc32_update(uint32 crc, uint8 byte)
923 {
924    return (crc << 8) ^ crc_table[byte ^ (crc >> 24)];
925 }
926
927
928 // used in setup, and for huffman that doesn't go fast path
929 static unsigned int bit_reverse(unsigned int n)
930 {
931   n = ((n & 0xAAAAAAAA) >>  1) | ((n & 0x55555555) << 1);
932   n = ((n & 0xCCCCCCCC) >>  2) | ((n & 0x33333333) << 2);
933   n = ((n & 0xF0F0F0F0) >>  4) | ((n & 0x0F0F0F0F) << 4);
934   n = ((n & 0xFF00FF00) >>  8) | ((n & 0x00FF00FF) << 8);
935   return (n >> 16) | (n << 16);
936 }
937
938 static float square(float x)
939 {
940    return x*x;
941 }
942
943 // this is a weird definition of log2() for which log2(1) = 1, log2(2) = 2, log2(4) = 3
944 // as required by the specification. fast(?) implementation from stb.h
945 // @OPTIMIZE: called multiple times per-packet with "constants"; move to setup
946 static int ilog(int32 n)
947 {
948    static signed char log2_4[16] = { 0,1,2,2,3,3,3,3,4,4,4,4,4,4,4,4 };
949
950    // 2 compares if n < 16, 3 compares otherwise (4 if signed or n > 1<<29)
951    if (n < (1U << 14))
952         if (n < (1U <<  4))        return     0 + log2_4[n      ];
953         else if (n < (1U <<  9))      return  5 + log2_4[n >>  5];
954              else                     return 10 + log2_4[n >> 10];
955    else if (n < (1U << 24))
956              if (n < (1U << 19))      return 15 + log2_4[n >> 15];
957              else                     return 20 + log2_4[n >> 20];
958         else if (n < (1U << 29))      return 25 + log2_4[n >> 25];
959              else if (n < (1U << 31)) return 30 + log2_4[n >> 30];
960                   else                return 0; // signed n returns 0
961 }
962
963 #ifndef M_PI
964   #define M_PI  3.14159265358979323846264f  // from CRC
965 #endif
966
967 // code length assigned to a value with no huffman encoding
968 #define NO_CODE   255
969
970 /////////////////////// LEAF SETUP FUNCTIONS //////////////////////////
971 //
972 // these functions are only called at setup, and only a few times
973 // per file
974
975 static float float32_unpack(uint32 x)
976 {
977    // from the specification
978    uint32 mantissa = x & 0x1fffff;
979    uint32 sign = x & 0x80000000;
980    uint32 exp = (x & 0x7fe00000) >> 21;
981    double res = sign ? -(double)mantissa : (double)mantissa;
982    return (float) ldexp((float)res, exp-788);
983 }
984
985
986 // zlib & jpeg huffman tables assume that the output symbols
987 // can either be arbitrarily arranged, or have monotonically
988 // increasing frequencies--they rely on the lengths being sorted;
989 // this makes for a very simple generation algorithm.
990 // vorbis allows a huffman table with non-sorted lengths. This
991 // requires a more sophisticated construction, since symbols in
992 // order do not map to huffman codes "in order".
993 static void add_entry(Codebook *c, uint32 huff_code, int symbol, int count, int len, uint32 *values)
994 {
995    if (!c->sparse) {
996       c->codewords      [symbol] = huff_code;
997    } else {
998       c->codewords       [count] = huff_code;
999       c->codeword_lengths[count] = len;
1000       values             [count] = symbol;
1001    }
1002 }
1003
1004 static int compute_codewords(Codebook *c, uint8 *len, int n, uint32 *values)
1005 {
1006    int i,k,m=0;
1007    uint32 available[32];
1008
1009    memset(available, 0, sizeof(available));
1010    // find the first entry
1011    for (k=0; k < n; ++k) if (len[k] < NO_CODE) break;
1012    if (k == n) { assert(c->sorted_entries == 0); return TRUE; }
1013    // add to the list
1014    add_entry(c, 0, k, m++, len[k], values);
1015    // add all available leaves
1016    for (i=1; i <= len[k]; ++i)
1017       available[i] = 1 << (32-i);
1018    // note that the above code treats the first case specially,
1019    // but it's really the same as the following code, so they
1020    // could probably be combined (except the initial code is 0,
1021    // and I use 0 in available[] to mean 'empty')
1022    for (i=k+1; i < n; ++i) {
1023       uint32 res;
1024       int z = len[i], y;
1025       if (z == NO_CODE) continue;
1026       // find lowest available leaf (should always be earliest,
1027       // which is what the specification calls for)
1028       // note that this property, and the fact we can never have
1029       // more than one free leaf at a given level, isn't totally
1030       // trivial to prove, but it seems true and the assert never
1031       // fires, so!
1032       while (z > 0 && !available[z]) --z;
1033       if (z == 0) { assert(0); return FALSE; }
1034       res = available[z];
1035       available[z] = 0;
1036       add_entry(c, bit_reverse(res), i, m++, len[i], values);
1037       // propogate availability up the tree
1038       if (z != len[i]) {
1039          for (y=len[i]; y > z; --y) {
1040             assert(available[y] == 0);
1041             available[y] = res + (1 << (32-y));
1042          }
1043       }
1044    }
1045    return TRUE;
1046 }
1047
1048 // accelerated huffman table allows fast O(1) match of all symbols
1049 // of length <= STB_VORBIS_FAST_HUFFMAN_LENGTH
1050 static void compute_accelerated_huffman(Codebook *c)
1051 {
1052    int i, len;
1053    for (i=0; i < FAST_HUFFMAN_TABLE_SIZE; ++i)
1054       c->fast_huffman[i] = -1;
1055
1056    len = c->sparse ? c->sorted_entries : c->entries;
1057    #ifdef STB_VORBIS_FAST_HUFFMAN_SHORT
1058    if (len > 32767) len = 32767; // largest possible value we can encode!
1059    #endif
1060    for (i=0; i < len; ++i) {
1061       if (c->codeword_lengths[i] <= STB_VORBIS_FAST_HUFFMAN_LENGTH) {
1062          uint32 z = c->sparse ? bit_reverse(c->sorted_codewords[i]) : c->codewords[i];
1063          // set table entries for all bit combinations in the higher bits
1064          while (z < FAST_HUFFMAN_TABLE_SIZE) {
1065              c->fast_huffman[z] = i;
1066              z += 1 << c->codeword_lengths[i];
1067          }
1068       }
1069    }
1070 }
1071
1072 static int uint32_compare(const void *p, const void *q)
1073 {
1074    uint32 x = * (uint32 *) p;
1075    uint32 y = * (uint32 *) q;
1076    return x < y ? -1 : x > y;
1077 }
1078
1079 static int include_in_sort(Codebook *c, uint8 len)
1080 {
1081    if (c->sparse) { assert(len != NO_CODE); return TRUE; }
1082    if (len == NO_CODE) return FALSE;
1083    if (len > STB_VORBIS_FAST_HUFFMAN_LENGTH) return TRUE;
1084    return FALSE;
1085 }
1086
1087 // if the fast table above doesn't work, we want to binary
1088 // search them... need to reverse the bits
1089 static void compute_sorted_huffman(Codebook *c, uint8 *lengths, uint32 *values)
1090 {
1091    int i, len;
1092    // build a list of all the entries
1093    // OPTIMIZATION: don't include the short ones, since they'll be caught by FAST_HUFFMAN.
1094    // this is kind of a frivolous optimization--I don't see any performance improvement,
1095    // but it's like 4 extra lines of code, so.
1096    if (!c->sparse) {
1097       int k = 0;
1098       for (i=0; i < c->entries; ++i)
1099          if (include_in_sort(c, lengths[i])) 
1100             c->sorted_codewords[k++] = bit_reverse(c->codewords[i]);
1101       assert(k == c->sorted_entries);
1102    } else {
1103       for (i=0; i < c->sorted_entries; ++i)
1104          c->sorted_codewords[i] = bit_reverse(c->codewords[i]);
1105    }
1106
1107    qsort(c->sorted_codewords, c->sorted_entries, sizeof(c->sorted_codewords[0]), uint32_compare);
1108    c->sorted_codewords[c->sorted_entries] = 0xffffffff;
1109
1110    len = c->sparse ? c->sorted_entries : c->entries;
1111    // now we need to indicate how they correspond; we could either
1112    //   #1: sort a different data structure that says who they correspond to
1113    //   #2: for each sorted entry, search the original list to find who corresponds
1114    //   #3: for each original entry, find the sorted entry
1115    // #1 requires extra storage, #2 is slow, #3 can use binary search!
1116    for (i=0; i < len; ++i) {
1117       int huff_len = c->sparse ? lengths[values[i]] : lengths[i];
1118       if (include_in_sort(c,huff_len)) {
1119          uint32 code = bit_reverse(c->codewords[i]);
1120          int x=0, n=c->sorted_entries;
1121          while (n > 1) {
1122             // invariant: sc[x] <= code < sc[x+n]
1123             int m = x + (n >> 1);
1124             if (c->sorted_codewords[m] <= code) {
1125                x = m;
1126                n -= (n>>1);
1127             } else {
1128                n >>= 1;
1129             }
1130          }
1131          assert(c->sorted_codewords[x] == code);
1132          if (c->sparse) {
1133             c->sorted_values[x] = values[i];
1134             c->codeword_lengths[x] = huff_len;
1135          } else {
1136             c->sorted_values[x] = i;
1137          }
1138       }
1139    }
1140 }
1141
1142 // only run while parsing the header (3 times)
1143 static int vorbis_validate(uint8 *data)
1144 {
1145    static uint8 vorbis[6] = { 'v', 'o', 'r', 'b', 'i', 's' };
1146    return memcmp(data, vorbis, 6) == 0;
1147 }
1148
1149 // called from setup only, once per code book
1150 // (formula implied by specification)
1151 static int lookup1_values(int entries, int dim)
1152 {
1153    int r = (int) floor(exp((float) log((float) entries) / dim));
1154    if ((int) floor(pow((float) r+1, dim)) <= entries)   // (int) cast for MinGW warning;
1155       ++r;                                              // floor() to avoid _ftol() when non-CRT
1156    assert(pow((float) r+1, dim) > entries);
1157    assert((int) floor(pow((float) r, dim)) <= entries); // (int),floor() as above
1158    return r;
1159 }
1160
1161 // called twice per file
1162 static void compute_twiddle_factors(int n, float *A, float *B, float *C)
1163 {
1164    int n4 = n >> 2, n8 = n >> 3;
1165    int k,k2;
1166
1167    for (k=k2=0; k < n4; ++k,k2+=2) {
1168       A[k2  ] = (float)  cos(4*k*M_PI/n);
1169       A[k2+1] = (float) -sin(4*k*M_PI/n);
1170       B[k2  ] = (float)  cos((k2+1)*M_PI/n/2) * 0.5f;
1171       B[k2+1] = (float)  sin((k2+1)*M_PI/n/2) * 0.5f;
1172    }
1173    for (k=k2=0; k < n8; ++k,k2+=2) {
1174       C[k2  ] = (float)  cos(2*(k2+1)*M_PI/n);
1175       C[k2+1] = (float) -sin(2*(k2+1)*M_PI/n);
1176    }
1177 }
1178
1179 static void compute_window(int n, float *window)
1180 {
1181    int n2 = n >> 1, i;
1182    for (i=0; i < n2; ++i)
1183       window[i] = (float) sin(0.5 * M_PI * square((float) sin((i - 0 + 0.5) / n2 * 0.5 * M_PI)));
1184 }
1185
1186 static void compute_bitreverse(int n, uint16 *rev)
1187 {
1188    int ld = ilog(n) - 1; // ilog is off-by-one from normal definitions
1189    int i, n8 = n >> 3;
1190    for (i=0; i < n8; ++i)
1191       rev[i] = (bit_reverse(i) >> (32-ld+3)) << 2;
1192 }
1193
1194 static int init_blocksize(vorb *f, int b, int n)
1195 {
1196    int n2 = n >> 1, n4 = n >> 2, n8 = n >> 3;
1197    f->A[b] = (float *) setup_malloc(f, sizeof(float) * n2);
1198    f->B[b] = (float *) setup_malloc(f, sizeof(float) * n2);
1199    f->C[b] = (float *) setup_malloc(f, sizeof(float) * n4);
1200    if (!f->A[b] || !f->B[b] || !f->C[b]) return error(f, VORBIS_outofmem);
1201    compute_twiddle_factors(n, f->A[b], f->B[b], f->C[b]);
1202    f->window[b] = (float *) setup_malloc(f, sizeof(float) * n2);
1203    if (!f->window[b]) return error(f, VORBIS_outofmem);
1204    compute_window(n, f->window[b]);
1205    f->bit_reverse[b] = (uint16 *) setup_malloc(f, sizeof(uint16) * n8);
1206    if (!f->bit_reverse[b]) return error(f, VORBIS_outofmem);
1207    compute_bitreverse(n, f->bit_reverse[b]);
1208    return TRUE;
1209 }
1210
1211 static void neighbors(uint16 *x, int n, int *plow, int *phigh)
1212 {
1213    int low = -1;
1214    int high = 65536;
1215    int i;
1216    for (i=0; i < n; ++i) {
1217       if (x[i] > low  && x[i] < x[n]) { *plow  = i; low = x[i]; }
1218       if (x[i] < high && x[i] > x[n]) { *phigh = i; high = x[i]; }
1219    }
1220 }
1221
1222 // this has been repurposed so y is now the original index instead of y
1223 typedef struct
1224 {
1225    uint16 x,y;
1226 } Point;
1227
1228 int point_compare(const void *p, const void *q)
1229 {
1230    Point *a = (Point *) p;
1231    Point *b = (Point *) q;
1232    return a->x < b->x ? -1 : a->x > b->x;
1233 }
1234
1235 //
1236 /////////////////////// END LEAF SETUP FUNCTIONS //////////////////////////
1237
1238
1239 #if defined(STB_VORBIS_NO_STDIO)
1240    #define USE_MEMORY(z)    TRUE
1241 #else
1242    #define USE_MEMORY(z)    ((z)->stream)
1243 #endif
1244
1245 static uint8 get8(vorb *z)
1246 {
1247    if (USE_MEMORY(z)) {
1248       if (z->stream >= z->stream_end) { z->eof = TRUE; return 0; }
1249       return *z->stream++;
1250    }
1251
1252    #ifndef STB_VORBIS_NO_STDIO
1253    {
1254    int c = fgetc(z->f);
1255    if (c == EOF) { z->eof = TRUE; return 0; }
1256    return c;
1257    }
1258    #endif
1259 }
1260
1261 static uint32 get32(vorb *f)
1262 {
1263    uint32 x;
1264    x = get8(f);
1265    x += get8(f) << 8;
1266    x += get8(f) << 16;
1267    x += get8(f) << 24;
1268    return x;
1269 }
1270
1271 static int getn(vorb *z, uint8 *data, int n)
1272 {
1273    if (USE_MEMORY(z)) {
1274       if (z->stream+n > z->stream_end) { z->eof = 1; return 0; }
1275       memcpy(data, z->stream, n);
1276       z->stream += n;
1277       return 1;
1278    }
1279
1280    #ifndef STB_VORBIS_NO_STDIO   
1281    if (fread(data, n, 1, z->f) == 1)
1282       return 1;
1283    else {
1284       z->eof = 1;
1285       return 0;
1286    }
1287    #endif
1288 }
1289
1290 static void skip(vorb *z, int n)
1291 {
1292    if (USE_MEMORY(z)) {
1293       z->stream += n;
1294       if (z->stream >= z->stream_end) z->eof = 1;
1295       return;
1296    }
1297    #ifndef STB_VORBIS_NO_STDIO
1298    {
1299       long x = ftell(z->f);
1300       fseek(z->f, x+n, SEEK_SET);
1301    }
1302    #endif
1303 }
1304
1305 static int set_file_offset(stb_vorbis *f, unsigned int loc)
1306 {
1307    #ifndef STB_VORBIS_NO_PUSHDATA_API
1308    if (f->push_mode) return 0;
1309    #endif
1310    f->eof = 0;
1311    if (USE_MEMORY(f)) {
1312       if (f->stream_start + loc >= f->stream_end || f->stream_start + loc < f->stream_start) {
1313          f->stream = f->stream_end;
1314          f->eof = 1;
1315          return 0;
1316       } else {
1317          f->stream = f->stream_start + loc;
1318          return 1;
1319       }
1320    }
1321    #ifndef STB_VORBIS_NO_STDIO
1322    if (loc + f->f_start < loc || loc >= 0x80000000) {
1323       loc = 0x7fffffff;
1324       f->eof = 1;
1325    } else {
1326       loc += f->f_start;
1327    }
1328    if (!fseek(f->f, loc, SEEK_SET))
1329       return 1;
1330    f->eof = 1;
1331    fseek(f->f, f->f_start, SEEK_END);
1332    return 0;
1333    #endif
1334 }
1335
1336
1337 static uint8 ogg_page_header[4] = { 0x4f, 0x67, 0x67, 0x53 };
1338
1339 static int capture_pattern(vorb *f)
1340 {
1341    if (0x4f != get8(f)) return FALSE;
1342    if (0x67 != get8(f)) return FALSE;
1343    if (0x67 != get8(f)) return FALSE;
1344    if (0x53 != get8(f)) return FALSE;
1345    return TRUE;
1346 }
1347
1348 #define PAGEFLAG_continued_packet   1
1349 #define PAGEFLAG_first_page         2
1350 #define PAGEFLAG_last_page          4
1351
1352 static int start_page_no_capturepattern(vorb *f)
1353 {
1354    uint32 loc0,loc1,n,i;
1355    // stream structure version
1356    if (0 != get8(f)) return error(f, VORBIS_invalid_stream_structure_version);
1357    // header flag
1358    f->page_flag = get8(f);
1359    // absolute granule position
1360    loc0 = get32(f); 
1361    loc1 = get32(f);
1362    // @TODO: validate loc0,loc1 as valid positions?
1363    // stream serial number -- vorbis doesn't interleave, so discard
1364    get32(f);
1365    //if (f->serial != get32(f)) return error(f, VORBIS_incorrect_stream_serial_number);
1366    // page sequence number
1367    n = get32(f);
1368    f->last_page = n;
1369    // CRC32
1370    get32(f);
1371    // page_segments
1372    f->segment_count = get8(f);
1373    if (!getn(f, f->segments, f->segment_count))
1374       return error(f, VORBIS_unexpected_eof);
1375    // assume we _don't_ know any the sample position of any segments
1376    f->end_seg_with_known_loc = -2;
1377    if (loc0 != ~0 || loc1 != ~0) {
1378       // determine which packet is the last one that will complete
1379       for (i=f->segment_count-1; i >= 0; --i)
1380          if (f->segments[i] < 255)
1381             break;
1382       // 'i' is now the index of the _last_ segment of a packet that ends
1383       if (i >= 0) {
1384          f->end_seg_with_known_loc = i;
1385          f->known_loc_for_packet   = loc0;
1386       }
1387    }
1388    if (f->first_decode) {
1389       int i,len;
1390       ProbedPage p;
1391       len = 0;
1392       for (i=0; i < f->segment_count; ++i)
1393          len += f->segments[i];
1394       len += 27 + f->segment_count;
1395       p.page_start = f->first_audio_page_offset;
1396       p.page_end = p.page_start + len;
1397       p.after_previous_page_start = p.page_start;
1398       p.first_decoded_sample = 0;
1399       p.last_decoded_sample = loc0;
1400       f->p_first = p;
1401    }
1402    f->next_seg = 0;
1403    return TRUE;
1404 }
1405
1406 static int start_page(vorb *f)
1407 {
1408    if (!capture_pattern(f)) return error(f, VORBIS_missing_capture_pattern);
1409    return start_page_no_capturepattern(f);
1410 }
1411
1412 static int start_packet(vorb *f)
1413 {
1414    while (f->next_seg == -1) {
1415       if (!start_page(f)) return FALSE;
1416       if (f->page_flag & PAGEFLAG_continued_packet)
1417          return error(f, VORBIS_continued_packet_flag_invalid);
1418    }
1419    f->last_seg = FALSE;
1420    f->valid_bits = 0;
1421    f->packet_bytes = 0;
1422    f->bytes_in_seg = 0;
1423    // f->next_seg is now valid
1424    return TRUE;
1425 }
1426
1427 static int maybe_start_packet(vorb *f)
1428 {
1429    if (f->next_seg == -1) {
1430       int x = get8(f);
1431       if (f->eof) return FALSE; // EOF at page boundary is not an error!
1432       if (0x4f != x      ) return error(f, VORBIS_missing_capture_pattern);
1433       if (0x67 != get8(f)) return error(f, VORBIS_missing_capture_pattern);
1434       if (0x67 != get8(f)) return error(f, VORBIS_missing_capture_pattern);
1435       if (0x53 != get8(f)) return error(f, VORBIS_missing_capture_pattern);
1436       if (!start_page_no_capturepattern(f)) return FALSE;
1437       if (f->page_flag & PAGEFLAG_continued_packet) {
1438          // set up enough state that we can read this packet if we want,
1439          // e.g. during recovery
1440          f->last_seg = FALSE;
1441          f->bytes_in_seg = 0;
1442          return error(f, VORBIS_continued_packet_flag_invalid);
1443       }
1444    }
1445    return start_packet(f);
1446 }
1447
1448 static int next_segment(vorb *f)
1449 {
1450    int len;
1451    if (f->last_seg) return 0;
1452    if (f->next_seg == -1) {
1453       f->last_seg_which = f->segment_count-1; // in case start_page fails
1454       if (!start_page(f)) { f->last_seg = 1; return 0; }
1455       if (!(f->page_flag & PAGEFLAG_continued_packet)) return error(f, VORBIS_continued_packet_flag_invalid);
1456    }
1457    len = f->segments[f->next_seg++];
1458    if (len < 255) {
1459       f->last_seg = TRUE;
1460       f->last_seg_which = f->next_seg-1;
1461    }
1462    if (f->next_seg >= f->segment_count)
1463       f->next_seg = -1;
1464    assert(f->bytes_in_seg == 0);
1465    f->bytes_in_seg = len;
1466    return len;
1467 }
1468
1469 #define EOP    (-1)
1470 #define INVALID_BITS  (-1)
1471
1472 static int get8_packet_raw(vorb *f)
1473 {
1474    if (!f->bytes_in_seg)
1475       if (f->last_seg) return EOP;
1476       else if (!next_segment(f)) return EOP;
1477    assert(f->bytes_in_seg > 0);
1478    --f->bytes_in_seg;
1479    ++f->packet_bytes;
1480    return get8(f);
1481 }
1482
1483 static int get8_packet(vorb *f)
1484 {
1485    int x = get8_packet_raw(f);
1486    f->valid_bits = 0;
1487    return x;
1488 }
1489
1490 static void flush_packet(vorb *f)
1491 {
1492    while (get8_packet_raw(f) != EOP);
1493 }
1494
1495 // @OPTIMIZE: this is the secondary bit decoder, so it's probably not as important
1496 // as the huffman decoder?
1497 static uint32 get_bits(vorb *f, int n)
1498 {
1499    uint32 z;
1500
1501    if (f->valid_bits < 0) return 0;
1502    if (f->valid_bits < n) {
1503       if (n > 24) {
1504          // the accumulator technique below would not work correctly in this case
1505          z = get_bits(f, 24);
1506          z += get_bits(f, n-24) << 24;
1507          return z;
1508       }
1509       if (f->valid_bits == 0) f->acc = 0;
1510       while (f->valid_bits < n) {
1511          int z = get8_packet_raw(f);
1512          if (z == EOP) {
1513             f->valid_bits = INVALID_BITS;
1514             return 0;
1515          }
1516          f->acc += z << f->valid_bits;
1517          f->valid_bits += 8;
1518       }
1519    }
1520    if (f->valid_bits < 0) return 0;
1521    z = f->acc & ((1 << n)-1);
1522    f->acc >>= n;
1523    f->valid_bits -= n;
1524    return z;
1525 }
1526
1527 static int32 get_bits_signed(vorb *f, int n)
1528 {
1529    uint32 z = get_bits(f, n);
1530    if (z & (1 << (n-1)))
1531       z += ~((1 << n) - 1);
1532    return (int32) z;
1533 }
1534
1535 // @OPTIMIZE: primary accumulator for huffman
1536 // expand the buffer to as many bits as possible without reading off end of packet
1537 // it might be nice to allow f->valid_bits and f->acc to be stored in registers,
1538 // e.g. cache them locally and decode locally
1539 static __forceinline void prep_huffman(vorb *f)
1540 {
1541    if (f->valid_bits <= 24) {
1542       if (f->valid_bits == 0) f->acc = 0;
1543       do {
1544          int z;
1545          if (f->last_seg && !f->bytes_in_seg) return;
1546          z = get8_packet_raw(f);
1547          if (z == EOP) return;
1548          f->acc += z << f->valid_bits;
1549          f->valid_bits += 8;
1550       } while (f->valid_bits <= 24);
1551    }
1552 }
1553
1554 enum
1555 {
1556    VORBIS_packet_id = 1,
1557    VORBIS_packet_comment = 3,
1558    VORBIS_packet_setup = 5,
1559 };
1560
1561 static int codebook_decode_scalar_raw(vorb *f, Codebook *c)
1562 {
1563    int i;
1564    prep_huffman(f);
1565
1566    assert(c->sorted_codewords || c->codewords);
1567    // cases to use binary search: sorted_codewords && !c->codewords
1568    //                             sorted_codewords && c->entries > 8
1569    if (c->entries > 8 ? c->sorted_codewords!=NULL : !c->codewords) {
1570       // binary search
1571       uint32 code = bit_reverse(f->acc);
1572       int x=0, n=c->sorted_entries, len;
1573
1574       while (n > 1) {
1575          // invariant: sc[x] <= code < sc[x+n]
1576          int m = x + (n >> 1);
1577          if (c->sorted_codewords[m] <= code) {
1578             x = m;
1579             n -= (n>>1);
1580          } else {
1581             n >>= 1;
1582          }
1583       }
1584       // x is now the sorted index
1585       if (!c->sparse) x = c->sorted_values[x];
1586       // x is now sorted index if sparse, or symbol otherwise
1587       len = c->codeword_lengths[x];
1588       if (f->valid_bits >= len) {
1589          f->acc >>= len;
1590          f->valid_bits -= len;
1591          return x;
1592       }
1593
1594       f->valid_bits = 0;
1595       return -1;
1596    }
1597
1598    // if small, linear search
1599    assert(!c->sparse);
1600    for (i=0; i < c->entries; ++i) {
1601       if (c->codeword_lengths[i] == NO_CODE) continue;
1602       if (c->codewords[i] == (f->acc & ((1 << c->codeword_lengths[i])-1))) {
1603          if (f->valid_bits >= c->codeword_lengths[i]) {
1604             f->acc >>= c->codeword_lengths[i];
1605             f->valid_bits -= c->codeword_lengths[i];
1606             return i;
1607          }
1608          f->valid_bits = 0;
1609          return -1;
1610       }
1611    }
1612
1613    error(f, VORBIS_invalid_stream);
1614    f->valid_bits = 0;
1615    return -1;
1616 }
1617
1618 static int codebook_decode_scalar(vorb *f, Codebook *c)
1619 {
1620    int i;
1621    if (f->valid_bits < STB_VORBIS_FAST_HUFFMAN_LENGTH)
1622       prep_huffman(f);
1623    // fast huffman table lookup
1624    i = f->acc & FAST_HUFFMAN_TABLE_MASK;
1625    i = c->fast_huffman[i];
1626    if (i >= 0) {
1627       f->acc >>= c->codeword_lengths[i];
1628       f->valid_bits -= c->codeword_lengths[i];
1629       if (f->valid_bits < 0) { f->valid_bits = 0; return -1; }
1630       return i;
1631    }
1632    return codebook_decode_scalar_raw(f,c);
1633 }
1634
1635 #ifndef STB_VORBIS_NO_INLINE_DECODE
1636
1637 #define DECODE_RAW(var, f,c)                                  \
1638    if (f->valid_bits < STB_VORBIS_FAST_HUFFMAN_LENGTH)        \
1639       prep_huffman(f);                                        \
1640    var = f->acc & FAST_HUFFMAN_TABLE_MASK;                    \
1641    var = c->fast_huffman[var];                                \
1642    if (var >= 0) {                                            \
1643       int n = c->codeword_lengths[var];                       \
1644       f->acc >>= n;                                           \
1645       f->valid_bits -= n;                                     \
1646       if (f->valid_bits < 0) { f->valid_bits = 0; var = -1; } \
1647    } else {                                                   \
1648       var = codebook_decode_scalar_raw(f,c);                  \
1649    }
1650
1651 #else
1652
1653 #define DECODE_RAW(var,f,c)    var = codebook_decode_scalar(f,c);
1654
1655 #endif
1656
1657 #define DECODE(var,f,c)                                       \
1658    DECODE_RAW(var,f,c)                                        \
1659    if (c->sparse) var = c->sorted_values[var];
1660
1661 #ifndef STB_VORBIS_DIVIDES_IN_CODEBOOK
1662   #define DECODE_VQ(var,f,c)   DECODE_RAW(var,f,c)
1663 #else
1664   #define DECODE_VQ(var,f,c)   DECODE(var,f,c)
1665 #endif
1666
1667
1668
1669
1670
1671
1672 // CODEBOOK_ELEMENT_FAST is an optimization for the CODEBOOK_FLOATS case
1673 // where we avoid one addition
1674 #ifndef STB_VORBIS_CODEBOOK_FLOATS
1675    #define CODEBOOK_ELEMENT(c,off)          (c->multiplicands[off] * c->delta_value + c->minimum_value)
1676    #define CODEBOOK_ELEMENT_FAST(c,off)     (c->multiplicands[off] * c->delta_value)
1677    #define CODEBOOK_ELEMENT_BASE(c)         (c->minimum_value)
1678 #else
1679    #define CODEBOOK_ELEMENT(c,off)          (c->multiplicands[off])
1680    #define CODEBOOK_ELEMENT_FAST(c,off)     (c->multiplicands[off])
1681    #define CODEBOOK_ELEMENT_BASE(c)         (0)
1682 #endif
1683
1684 static int codebook_decode_start(vorb *f, Codebook *c, int len)
1685 {
1686    int z = -1;
1687
1688    // type 0 is only legal in a scalar context
1689    if (c->lookup_type == 0)
1690       error(f, VORBIS_invalid_stream);
1691    else {
1692       DECODE_VQ(z,f,c);
1693       if (c->sparse) assert(z < c->sorted_entries);
1694       if (z < 0) {  // check for EOP
1695          if (!f->bytes_in_seg)
1696             if (f->last_seg)
1697                return z;
1698          error(f, VORBIS_invalid_stream);
1699       }
1700    }
1701    return z;
1702 }
1703
1704 static int codebook_decode(vorb *f, Codebook *c, float *output, int len)
1705 {
1706    int i,z = codebook_decode_start(f,c,len);
1707    if (z < 0) return FALSE;
1708    if (len > c->dimensions) len = c->dimensions;
1709
1710 #ifdef STB_VORBIS_DIVIDES_IN_CODEBOOK
1711    if (c->lookup_type == 1) {
1712       float last = CODEBOOK_ELEMENT_BASE(c);
1713       int div = 1;
1714       for (i=0; i < len; ++i) {
1715          int off = (z / div) % c->lookup_values;
1716          float val = CODEBOOK_ELEMENT_FAST(c,off) + last;
1717          output[i] += val;
1718          if (c->sequence_p) last = val + c->minimum_value;
1719          div *= c->lookup_values;
1720       }
1721       return TRUE;
1722    }
1723 #endif
1724
1725    z *= c->dimensions;
1726    if (c->sequence_p) {
1727       float last = CODEBOOK_ELEMENT_BASE(c);
1728       for (i=0; i < len; ++i) {
1729          float val = CODEBOOK_ELEMENT_FAST(c,z+i) + last;
1730          output[i] += val;
1731          last = val + c->minimum_value;
1732       }
1733    } else {
1734       float last = CODEBOOK_ELEMENT_BASE(c);
1735       for (i=0; i < len; ++i) {
1736          output[i] += CODEBOOK_ELEMENT_FAST(c,z+i) + last;
1737       }
1738    }
1739
1740    return TRUE;
1741 }
1742
1743 static int codebook_decode_step(vorb *f, Codebook *c, float *output, int len, int step)
1744 {
1745    int i,z = codebook_decode_start(f,c,len);
1746    float last = CODEBOOK_ELEMENT_BASE(c);
1747    if (z < 0) return FALSE;
1748    if (len > c->dimensions) len = c->dimensions;
1749
1750 #ifdef STB_VORBIS_DIVIDES_IN_CODEBOOK
1751    if (c->lookup_type == 1) {
1752       int div = 1;
1753       for (i=0; i < len; ++i) {
1754          int off = (z / div) % c->lookup_values;
1755          float val = CODEBOOK_ELEMENT_FAST(c,off) + last;
1756          output[i*step] += val;
1757          if (c->sequence_p) last = val;
1758          div *= c->lookup_values;
1759       }
1760       return TRUE;
1761    }
1762 #endif
1763
1764    z *= c->dimensions;
1765    for (i=0; i < len; ++i) {
1766       float val = CODEBOOK_ELEMENT_FAST(c,z+i) + last;
1767       output[i*step] += val;
1768       if (c->sequence_p) last = val;
1769    }
1770
1771    return TRUE;
1772 }
1773
1774 static int codebook_decode_deinterleave_repeat(vorb *f, Codebook *c, float **outputs, int ch, int *c_inter_p, int *p_inter_p, int len, int total_decode)
1775 {
1776    int c_inter = *c_inter_p;
1777    int p_inter = *p_inter_p;
1778    int i,z, effective = c->dimensions;
1779
1780    // type 0 is only legal in a scalar context
1781    if (c->lookup_type == 0)   return error(f, VORBIS_invalid_stream);
1782
1783    while (total_decode > 0) {
1784       float last = CODEBOOK_ELEMENT_BASE(c);
1785       DECODE_VQ(z,f,c);
1786       #ifndef STB_VORBIS_DIVIDES_IN_CODEBOOK
1787       assert(!c->sparse || z < c->sorted_entries);
1788       #endif
1789       if (z < 0) {
1790          if (!f->bytes_in_seg)
1791             if (f->last_seg) return FALSE;
1792          return error(f, VORBIS_invalid_stream);
1793       }
1794
1795       // if this will take us off the end of the buffers, stop short!
1796       // we check by computing the length of the virtual interleaved
1797       // buffer (len*ch), our current offset within it (p_inter*ch)+(c_inter),
1798       // and the length we'll be using (effective)
1799       if (c_inter + p_inter*ch + effective > len * ch) {
1800          effective = len*ch - (p_inter*ch - c_inter);
1801       }
1802
1803    #ifdef STB_VORBIS_DIVIDES_IN_CODEBOOK
1804       if (c->lookup_type == 1) {
1805          int div = 1;
1806          for (i=0; i < effective; ++i) {
1807             int off = (z / div) % c->lookup_values;
1808             float val = CODEBOOK_ELEMENT_FAST(c,off) + last;
1809             outputs[c_inter][p_inter] += val;
1810             if (++c_inter == ch) { c_inter = 0; ++p_inter; }
1811             if (c->sequence_p) last = val;
1812             div *= c->lookup_values;
1813          }
1814       } else
1815    #endif
1816       {
1817          z *= c->dimensions;
1818          if (c->sequence_p) {
1819             for (i=0; i < effective; ++i) {
1820                float val = CODEBOOK_ELEMENT_FAST(c,z+i) + last;
1821                outputs[c_inter][p_inter] += val;
1822                if (++c_inter == ch) { c_inter = 0; ++p_inter; }
1823                last = val;
1824             }
1825          } else {
1826             for (i=0; i < effective; ++i) {
1827                float val = CODEBOOK_ELEMENT_FAST(c,z+i) + last;
1828                outputs[c_inter][p_inter] += val;
1829                if (++c_inter == ch) { c_inter = 0; ++p_inter; }
1830             }
1831          }
1832       }
1833
1834       total_decode -= effective;
1835    }
1836    *c_inter_p = c_inter;
1837    *p_inter_p = p_inter;
1838    return TRUE;
1839 }
1840
1841 #ifndef STB_VORBIS_DIVIDES_IN_CODEBOOK
1842 static int codebook_decode_deinterleave_repeat_2(vorb *f, Codebook *c, float **outputs, int *c_inter_p, int *p_inter_p, int len, int total_decode)
1843 {
1844    int c_inter = *c_inter_p;
1845    int p_inter = *p_inter_p;
1846    int i,z, effective = c->dimensions;
1847
1848    // type 0 is only legal in a scalar context
1849    if (c->lookup_type == 0)   return error(f, VORBIS_invalid_stream);
1850
1851    while (total_decode > 0) {
1852       float last = CODEBOOK_ELEMENT_BASE(c);
1853       DECODE_VQ(z,f,c);
1854
1855       if (z < 0) {
1856          if (!f->bytes_in_seg)
1857             if (f->last_seg) return FALSE;
1858          return error(f, VORBIS_invalid_stream);
1859       }
1860
1861       // if this will take us off the end of the buffers, stop short!
1862       // we check by computing the length of the virtual interleaved
1863       // buffer (len*ch), our current offset within it (p_inter*ch)+(c_inter),
1864       // and the length we'll be using (effective)
1865       if (c_inter + p_inter*2 + effective > len * 2) {
1866          effective = len*2 - (p_inter*2 - c_inter);
1867       }
1868
1869       {
1870          z *= c->dimensions;
1871          stb_prof(11);
1872          if (c->sequence_p) {
1873             // haven't optimized this case because I don't have any examples
1874             for (i=0; i < effective; ++i) {
1875                float val = CODEBOOK_ELEMENT_FAST(c,z+i) + last;
1876                outputs[c_inter][p_inter] += val;
1877                if (++c_inter == 2) { c_inter = 0; ++p_inter; }
1878                last = val;
1879             }
1880          } else {
1881             i=0;
1882             if (c_inter == 1) {
1883                float val = CODEBOOK_ELEMENT_FAST(c,z+i) + last;
1884                outputs[c_inter][p_inter] += val;
1885                c_inter = 0; ++p_inter;
1886                ++i;
1887             }
1888             {
1889                float *z0 = outputs[0];
1890                float *z1 = outputs[1];
1891                for (; i+1 < effective;) {
1892                   z0[p_inter] += CODEBOOK_ELEMENT_FAST(c,z+i) + last;
1893                   z1[p_inter] += CODEBOOK_ELEMENT_FAST(c,z+i+1) + last;
1894                   ++p_inter;
1895                   i += 2;
1896                }
1897             }
1898             if (i < effective) {
1899                float val = CODEBOOK_ELEMENT_FAST(c,z+i) + last;
1900                outputs[c_inter][p_inter] += val;
1901                if (++c_inter == 2) { c_inter = 0; ++p_inter; }
1902             }
1903          }
1904       }
1905
1906       total_decode -= effective;
1907    }
1908    *c_inter_p = c_inter;
1909    *p_inter_p = p_inter;
1910    return TRUE;
1911 }
1912 #endif
1913
1914 static int predict_point(int x, int x0, int x1, int y0, int y1)
1915 {
1916    int dy = y1 - y0;
1917    int adx = x1 - x0;
1918    // @OPTIMIZE: force int division to round in the right direction... is this necessary on x86?
1919    int err = abs(dy) * (x - x0);
1920    int off = err / adx;
1921    return dy < 0 ? y0 - off : y0 + off;
1922 }
1923
1924 // the following table is block-copied from the specification
1925 static float inverse_db_table[256] =
1926 {
1927   1.0649863e-07f, 1.1341951e-07f, 1.2079015e-07f, 1.2863978e-07f, 
1928   1.3699951e-07f, 1.4590251e-07f, 1.5538408e-07f, 1.6548181e-07f, 
1929   1.7623575e-07f, 1.8768855e-07f, 1.9988561e-07f, 2.1287530e-07f, 
1930   2.2670913e-07f, 2.4144197e-07f, 2.5713223e-07f, 2.7384213e-07f, 
1931   2.9163793e-07f, 3.1059021e-07f, 3.3077411e-07f, 3.5226968e-07f, 
1932   3.7516214e-07f, 3.9954229e-07f, 4.2550680e-07f, 4.5315863e-07f, 
1933   4.8260743e-07f, 5.1396998e-07f, 5.4737065e-07f, 5.8294187e-07f, 
1934   6.2082472e-07f, 6.6116941e-07f, 7.0413592e-07f, 7.4989464e-07f, 
1935   7.9862701e-07f, 8.5052630e-07f, 9.0579828e-07f, 9.6466216e-07f, 
1936   1.0273513e-06f, 1.0941144e-06f, 1.1652161e-06f, 1.2409384e-06f, 
1937   1.3215816e-06f, 1.4074654e-06f, 1.4989305e-06f, 1.5963394e-06f, 
1938   1.7000785e-06f, 1.8105592e-06f, 1.9282195e-06f, 2.0535261e-06f, 
1939   2.1869758e-06f, 2.3290978e-06f, 2.4804557e-06f, 2.6416497e-06f, 
1940   2.8133190e-06f, 2.9961443e-06f, 3.1908506e-06f, 3.3982101e-06f, 
1941   3.6190449e-06f, 3.8542308e-06f, 4.1047004e-06f, 4.3714470e-06f, 
1942   4.6555282e-06f, 4.9580707e-06f, 5.2802740e-06f, 5.6234160e-06f, 
1943   5.9888572e-06f, 6.3780469e-06f, 6.7925283e-06f, 7.2339451e-06f, 
1944   7.7040476e-06f, 8.2047000e-06f, 8.7378876e-06f, 9.3057248e-06f, 
1945   9.9104632e-06f, 1.0554501e-05f, 1.1240392e-05f, 1.1970856e-05f, 
1946   1.2748789e-05f, 1.3577278e-05f, 1.4459606e-05f, 1.5399272e-05f, 
1947   1.6400004e-05f, 1.7465768e-05f, 1.8600792e-05f, 1.9809576e-05f, 
1948   2.1096914e-05f, 2.2467911e-05f, 2.3928002e-05f, 2.5482978e-05f, 
1949   2.7139006e-05f, 2.8902651e-05f, 3.0780908e-05f, 3.2781225e-05f, 
1950   3.4911534e-05f, 3.7180282e-05f, 3.9596466e-05f, 4.2169667e-05f, 
1951   4.4910090e-05f, 4.7828601e-05f, 5.0936773e-05f, 5.4246931e-05f, 
1952   5.7772202e-05f, 6.1526565e-05f, 6.5524908e-05f, 6.9783085e-05f, 
1953   7.4317983e-05f, 7.9147585e-05f, 8.4291040e-05f, 8.9768747e-05f, 
1954   9.5602426e-05f, 0.00010181521f, 0.00010843174f, 0.00011547824f, 
1955   0.00012298267f, 0.00013097477f, 0.00013948625f, 0.00014855085f, 
1956   0.00015820453f, 0.00016848555f, 0.00017943469f, 0.00019109536f, 
1957   0.00020351382f, 0.00021673929f, 0.00023082423f, 0.00024582449f, 
1958   0.00026179955f, 0.00027881276f, 0.00029693158f, 0.00031622787f, 
1959   0.00033677814f, 0.00035866388f, 0.00038197188f, 0.00040679456f, 
1960   0.00043323036f, 0.00046138411f, 0.00049136745f, 0.00052329927f, 
1961   0.00055730621f, 0.00059352311f, 0.00063209358f, 0.00067317058f, 
1962   0.00071691700f, 0.00076350630f, 0.00081312324f, 0.00086596457f, 
1963   0.00092223983f, 0.00098217216f, 0.0010459992f,  0.0011139742f, 
1964   0.0011863665f,  0.0012634633f,  0.0013455702f,  0.0014330129f, 
1965   0.0015261382f,  0.0016253153f,  0.0017309374f,  0.0018434235f, 
1966   0.0019632195f,  0.0020908006f,  0.0022266726f,  0.0023713743f, 
1967   0.0025254795f,  0.0026895994f,  0.0028643847f,  0.0030505286f, 
1968   0.0032487691f,  0.0034598925f,  0.0036847358f,  0.0039241906f, 
1969   0.0041792066f,  0.0044507950f,  0.0047400328f,  0.0050480668f, 
1970   0.0053761186f,  0.0057254891f,  0.0060975636f,  0.0064938176f, 
1971   0.0069158225f,  0.0073652516f,  0.0078438871f,  0.0083536271f, 
1972   0.0088964928f,  0.009474637f,   0.010090352f,   0.010746080f, 
1973   0.011444421f,   0.012188144f,   0.012980198f,   0.013823725f, 
1974   0.014722068f,   0.015678791f,   0.016697687f,   0.017782797f, 
1975   0.018938423f,   0.020169149f,   0.021479854f,   0.022875735f, 
1976   0.024362330f,   0.025945531f,   0.027631618f,   0.029427276f, 
1977   0.031339626f,   0.033376252f,   0.035545228f,   0.037855157f, 
1978   0.040315199f,   0.042935108f,   0.045725273f,   0.048696758f, 
1979   0.051861348f,   0.055231591f,   0.058820850f,   0.062643361f, 
1980   0.066714279f,   0.071049749f,   0.075666962f,   0.080584227f, 
1981   0.085821044f,   0.091398179f,   0.097337747f,   0.10366330f, 
1982   0.11039993f,    0.11757434f,    0.12521498f,    0.13335215f, 
1983   0.14201813f,    0.15124727f,    0.16107617f,    0.17154380f, 
1984   0.18269168f,    0.19456402f,    0.20720788f,    0.22067342f, 
1985   0.23501402f,    0.25028656f,    0.26655159f,    0.28387361f, 
1986   0.30232132f,    0.32196786f,    0.34289114f,    0.36517414f, 
1987   0.38890521f,    0.41417847f,    0.44109412f,    0.46975890f, 
1988   0.50028648f,    0.53279791f,    0.56742212f,    0.60429640f, 
1989   0.64356699f,    0.68538959f,    0.72993007f,    0.77736504f, 
1990   0.82788260f,    0.88168307f,    0.9389798f,     1.0f
1991 };
1992
1993
1994 // @OPTIMIZE: if you want to replace this bresenham line-drawing routine,
1995 // note that you must produce bit-identical output to decode correctly;
1996 // this specific sequence of operations is specified in the spec (it's
1997 // drawing integer-quantized frequency-space lines that the encoder
1998 // expects to be exactly the same)
1999 //     ... also, isn't the whole point of Bresenham's algorithm to NOT
2000 // have to divide in the setup? sigh.
2001 #ifndef STB_VORBIS_NO_DEFER_FLOOR
2002 #define LINE_OP(a,b)   a *= b
2003 #else
2004 #define LINE_OP(a,b)   a = b
2005 #endif
2006
2007 #ifdef STB_VORBIS_DIVIDE_TABLE
2008 #define DIVTAB_NUMER   32
2009 #define DIVTAB_DENOM   64
2010 int8 integer_divide_table[DIVTAB_NUMER][DIVTAB_DENOM]; // 2KB
2011 #endif
2012
2013 static __forceinline void draw_line(float *output, int x0, int y0, int x1, int y1, int n)
2014 {
2015    int dy = y1 - y0;
2016    int adx = x1 - x0;
2017    int ady = abs(dy);
2018    int base;
2019    int x=x0,y=y0;
2020    int err = 0;
2021    int sy;
2022
2023 #ifdef STB_VORBIS_DIVIDE_TABLE
2024    if (adx < DIVTAB_DENOM && ady < DIVTAB_NUMER) {
2025       if (dy < 0) {
2026          base = -integer_divide_table[ady][adx];
2027          sy = base-1;
2028       } else {
2029          base =  integer_divide_table[ady][adx];
2030          sy = base+1;
2031       }
2032    } else {
2033       base = dy / adx;
2034       if (dy < 0)
2035          sy = base - 1;
2036       else
2037          sy = base+1;
2038    }
2039 #else
2040    base = dy / adx;
2041    if (dy < 0)
2042       sy = base - 1;
2043    else
2044       sy = base+1;
2045 #endif
2046    ady -= abs(base) * adx;
2047    if (x1 > n) x1 = n;
2048    LINE_OP(output[x], inverse_db_table[y]);
2049    for (++x; x < x1; ++x) {
2050       err += ady;
2051       if (err >= adx) {
2052          err -= adx;
2053          y += sy;
2054       } else
2055          y += base;
2056       LINE_OP(output[x], inverse_db_table[y]);
2057    }
2058 }
2059
2060 static int residue_decode(vorb *f, Codebook *book, float *target, int offset, int n, int rtype)
2061 {
2062    int k;
2063    if (rtype == 0) {
2064       int step = n / book->dimensions;
2065       for (k=0; k < step; ++k)
2066          if (!codebook_decode_step(f, book, target+offset+k, n-offset-k, step))
2067             return FALSE;
2068    } else {
2069       for (k=0; k < n; ) {
2070          if (!codebook_decode(f, book, target+offset, n-k))
2071             return FALSE;
2072          k += book->dimensions;
2073          offset += book->dimensions;
2074       }
2075    }
2076    return TRUE;
2077 }
2078
2079 static void decode_residue(vorb *f, float *residue_buffers[], int ch, int n, int rn, uint8 *do_not_decode)
2080 {
2081    int i,j,pass;
2082    Residue *r = f->residue_config + rn;
2083    int rtype = f->residue_types[rn];
2084    int c = r->classbook;
2085    int classwords = f->codebooks[c].dimensions;
2086    int n_read = r->end - r->begin;
2087    int part_read = n_read / r->part_size;
2088    int temp_alloc_point = temp_alloc_save(f);
2089    #ifndef STB_VORBIS_DIVIDES_IN_RESIDUE
2090    uint8 ***part_classdata = (uint8 ***) temp_block_array(f,f->channels, part_read * sizeof(**part_classdata));
2091    #else
2092    int **classifications = (int **) temp_block_array(f,f->channels, part_read * sizeof(**classifications));
2093    #endif
2094
2095    stb_prof(2);
2096    for (i=0; i < ch; ++i)
2097       if (!do_not_decode[i])
2098          memset(residue_buffers[i], 0, sizeof(float) * n);
2099
2100    if (rtype == 2 && ch != 1) {
2101       int len = ch * n;
2102       for (j=0; j < ch; ++j)
2103          if (!do_not_decode[j])
2104             break;
2105       if (j == ch)
2106          goto done;
2107
2108       stb_prof(3);
2109       for (pass=0; pass < 8; ++pass) {
2110          int pcount = 0, class_set = 0;
2111          if (ch == 2) {
2112             stb_prof(13);
2113             while (pcount < part_read) {
2114                int z = r->begin + pcount*r->part_size;
2115                int c_inter = (z & 1), p_inter = z>>1;
2116                if (pass == 0) {
2117                   Codebook *c = f->codebooks+r->classbook;
2118                   int q;
2119                   DECODE(q,f,c);
2120                   if (q == EOP) goto done;
2121                   #ifndef STB_VORBIS_DIVIDES_IN_RESIDUE
2122                   part_classdata[0][class_set] = r->classdata[q];
2123                   #else
2124                   for (i=classwords-1; i >= 0; --i) {
2125                      classifications[0][i+pcount] = q % r->classifications;
2126                      q /= r->classifications;
2127                   }
2128                   #endif
2129                }
2130                stb_prof(5);
2131                for (i=0; i < classwords && pcount < part_read; ++i, ++pcount) {
2132                   int z = r->begin + pcount*r->part_size;
2133                   #ifndef STB_VORBIS_DIVIDES_IN_RESIDUE
2134                   int c = part_classdata[0][class_set][i];
2135                   #else
2136                   int c = classifications[0][pcount];
2137                   #endif
2138                   int b = r->residue_books[c][pass];
2139                   if (b >= 0) {
2140                      Codebook *book = f->codebooks + b;
2141                      stb_prof(20);  // accounts for X time
2142                      #ifdef STB_VORBIS_DIVIDES_IN_CODEBOOK
2143                      if (!codebook_decode_deinterleave_repeat(f, book, residue_buffers, ch, &c_inter, &p_inter, n, r->part_size))
2144                         goto done;
2145                      #else
2146                      // saves 1%
2147                      if (!codebook_decode_deinterleave_repeat_2(f, book, residue_buffers, &c_inter, &p_inter, n, r->part_size))
2148                         goto done;
2149                      #endif
2150                      stb_prof(7);
2151                   } else {
2152                      z += r->part_size;
2153                      c_inter = z & 1;
2154                      p_inter = z >> 1;
2155                   }
2156                }
2157                stb_prof(8);
2158                #ifndef STB_VORBIS_DIVIDES_IN_RESIDUE
2159                ++class_set;
2160                #endif
2161             }
2162          } else if (ch == 1) {
2163             while (pcount < part_read) {
2164                int z = r->begin + pcount*r->part_size;
2165                int c_inter = 0, p_inter = z;
2166                if (pass == 0) {
2167                   Codebook *c = f->codebooks+r->classbook;
2168                   int q;
2169                   DECODE(q,f,c);
2170                   if (q == EOP) goto done;
2171                   #ifndef STB_VORBIS_DIVIDES_IN_RESIDUE
2172                   part_classdata[0][class_set] = r->classdata[q];
2173                   #else
2174                   for (i=classwords-1; i >= 0; --i) {
2175                      classifications[0][i+pcount] = q % r->classifications;
2176                      q /= r->classifications;
2177                   }
2178                   #endif
2179                }
2180                for (i=0; i < classwords && pcount < part_read; ++i, ++pcount) {
2181                   int z = r->begin + pcount*r->part_size;
2182                   #ifndef STB_VORBIS_DIVIDES_IN_RESIDUE
2183                   int c = part_classdata[0][class_set][i];
2184                   #else
2185                   int c = classifications[0][pcount];
2186                   #endif
2187                   int b = r->residue_books[c][pass];
2188                   if (b >= 0) {
2189                      Codebook *book = f->codebooks + b;
2190                      stb_prof(22);
2191                      if (!codebook_decode_deinterleave_repeat(f, book, residue_buffers, ch, &c_inter, &p_inter, n, r->part_size))
2192                         goto done;
2193                      stb_prof(3);
2194                   } else {
2195                      z += r->part_size;
2196                      c_inter = 0;
2197                      p_inter = z;
2198                   }
2199                }
2200                #ifndef STB_VORBIS_DIVIDES_IN_RESIDUE
2201                ++class_set;
2202                #endif
2203             }
2204          } else {
2205             while (pcount < part_read) {
2206                int z = r->begin + pcount*r->part_size;
2207                int c_inter = z % ch, p_inter = z/ch;
2208                if (pass == 0) {
2209                   Codebook *c = f->codebooks+r->classbook;
2210                   int q;
2211                   DECODE(q,f,c);
2212                   if (q == EOP) goto done;
2213                   #ifndef STB_VORBIS_DIVIDES_IN_RESIDUE
2214                   part_classdata[0][class_set] = r->classdata[q];
2215                   #else
2216                   for (i=classwords-1; i >= 0; --i) {
2217                      classifications[0][i+pcount] = q % r->classifications;
2218                      q /= r->classifications;
2219                   }
2220                   #endif
2221                }
2222                for (i=0; i < classwords && pcount < part_read; ++i, ++pcount) {
2223                   int z = r->begin + pcount*r->part_size;
2224                   #ifndef STB_VORBIS_DIVIDES_IN_RESIDUE
2225                   int c = part_classdata[0][class_set][i];
2226                   #else
2227                   int c = classifications[0][pcount];
2228                   #endif
2229                   int b = r->residue_books[c][pass];
2230                   if (b >= 0) {
2231                      Codebook *book = f->codebooks + b;
2232                      stb_prof(22);
2233                      if (!codebook_decode_deinterleave_repeat(f, book, residue_buffers, ch, &c_inter, &p_inter, n, r->part_size))
2234                         goto done;
2235                      stb_prof(3);
2236                   } else {
2237                      z += r->part_size;
2238                      c_inter = z % ch;
2239                      p_inter = z / ch;
2240                   }
2241                }
2242                #ifndef STB_VORBIS_DIVIDES_IN_RESIDUE
2243                ++class_set;
2244                #endif
2245             }
2246          }
2247       }
2248       goto done;
2249    }
2250    stb_prof(9);
2251
2252    for (pass=0; pass < 8; ++pass) {
2253       int pcount = 0, class_set=0;
2254       while (pcount < part_read) {
2255          if (pass == 0) {
2256             for (j=0; j < ch; ++j) {
2257                if (!do_not_decode[j]) {
2258                   Codebook *c = f->codebooks+r->classbook;
2259                   int temp;
2260                   DECODE(temp,f,c);
2261                   if (temp == EOP) goto done;
2262                   #ifndef STB_VORBIS_DIVIDES_IN_RESIDUE
2263                   part_classdata[j][class_set] = r->classdata[temp];
2264                   #else
2265                   for (i=classwords-1; i >= 0; --i) {
2266                      classifications[j][i+pcount] = temp % r->classifications;
2267                      temp /= r->classifications;
2268                   }
2269                   #endif
2270                }
2271             }
2272          }
2273          for (i=0; i < classwords && pcount < part_read; ++i, ++pcount) {
2274             for (j=0; j < ch; ++j) {
2275                if (!do_not_decode[j]) {
2276                   #ifndef STB_VORBIS_DIVIDES_IN_RESIDUE
2277                   int c = part_classdata[j][class_set][i];
2278                   #else
2279                   int c = classifications[j][pcount];
2280                   #endif
2281                   int b = r->residue_books[c][pass];
2282                   if (b >= 0) {
2283                      float *target = residue_buffers[j];
2284                      int offset = r->begin + pcount * r->part_size;
2285                      int n = r->part_size;
2286                      Codebook *book = f->codebooks + b;
2287                      if (!residue_decode(f, book, target, offset, n, rtype))
2288                         goto done;
2289                   }
2290                }
2291             }
2292          }
2293          #ifndef STB_VORBIS_DIVIDES_IN_RESIDUE
2294          ++class_set;
2295          #endif
2296       }
2297    }
2298   done:
2299    stb_prof(0);
2300    temp_alloc_restore(f,temp_alloc_point);
2301 }
2302
2303
2304 #if 0
2305 // slow way for debugging
2306 void inverse_mdct_slow(float *buffer, int n)
2307 {
2308    int i,j;
2309    int n2 = n >> 1;
2310    float *x = (float *) malloc(sizeof(*x) * n2);
2311    memcpy(x, buffer, sizeof(*x) * n2);
2312    for (i=0; i < n; ++i) {
2313       float acc = 0;
2314       for (j=0; j < n2; ++j)
2315          // formula from paper:
2316          //acc += n/4.0f * x[j] * (float) cos(M_PI / 2 / n * (2 * i + 1 + n/2.0)*(2*j+1));
2317          // formula from wikipedia
2318          //acc += 2.0f / n2 * x[j] * (float) cos(M_PI/n2 * (i + 0.5 + n2/2)*(j + 0.5));
2319          // these are equivalent, except the formula from the paper inverts the multiplier!
2320          // however, what actually works is NO MULTIPLIER!?!
2321          //acc += 64 * 2.0f / n2 * x[j] * (float) cos(M_PI/n2 * (i + 0.5 + n2/2)*(j + 0.5));
2322          acc += x[j] * (float) cos(M_PI / 2 / n * (2 * i + 1 + n/2.0)*(2*j+1));
2323       buffer[i] = acc;
2324    }
2325    free(x);
2326 }
2327 #elif 0
2328 // same as above, but just barely able to run in real time on modern machines
2329 void inverse_mdct_slow(float *buffer, int n, vorb *f, int blocktype)
2330 {
2331    float mcos[16384];
2332    int i,j;
2333    int n2 = n >> 1, nmask = (n << 2) -1;
2334    float *x = (float *) malloc(sizeof(*x) * n2);
2335    memcpy(x, buffer, sizeof(*x) * n2);
2336    for (i=0; i < 4*n; ++i)
2337       mcos[i] = (float) cos(M_PI / 2 * i / n);
2338
2339    for (i=0; i < n; ++i) {
2340       float acc = 0;
2341       for (j=0; j < n2; ++j)
2342          acc += x[j] * mcos[(2 * i + 1 + n2)*(2*j+1) & nmask];
2343       buffer[i] = acc;
2344    }
2345    free(x);
2346 }
2347 #else
2348 // transform to use a slow dct-iv; this is STILL basically trivial,
2349 // but only requires half as many ops
2350 void dct_iv_slow(float *buffer, int n)
2351 {
2352    float mcos[16384];
2353    float x[2048];
2354    int i,j;
2355    int n2 = n >> 1, nmask = (n << 3) - 1;
2356    memcpy(x, buffer, sizeof(*x) * n);
2357    for (i=0; i < 8*n; ++i)
2358       mcos[i] = (float) cos(M_PI / 4 * i / n);
2359    for (i=0; i < n; ++i) {
2360       float acc = 0;
2361       for (j=0; j < n; ++j)
2362          acc += x[j] * mcos[((2 * i + 1)*(2*j+1)) & nmask];
2363          //acc += x[j] * cos(M_PI / n * (i + 0.5) * (j + 0.5));
2364       buffer[i] = acc;
2365    }
2366 }
2367
2368 void inverse_mdct_slow(float *buffer, int n, vorb *f, int blocktype)
2369 {
2370    int i, n4 = n >> 2, n2 = n >> 1, n3_4 = n - n4;
2371    float temp[4096];
2372
2373    memcpy(temp, buffer, n2 * sizeof(float));
2374    dct_iv_slow(temp, n2);  // returns -c'-d, a-b'
2375
2376    for (i=0; i < n4  ; ++i) buffer[i] = temp[i+n4];            // a-b'
2377    for (   ; i < n3_4; ++i) buffer[i] = -temp[n3_4 - i - 1];   // b-a', c+d'
2378    for (   ; i < n   ; ++i) buffer[i] = -temp[i - n3_4];       // c'+d
2379 }
2380 #endif
2381
2382 #ifndef LIBVORBIS_MDCT
2383 #define LIBVORBIS_MDCT 0
2384 #endif
2385
2386 #if LIBVORBIS_MDCT
2387 // directly call the vorbis MDCT using an interface documented
2388 // by Jeff Roberts... useful for performance comparison
2389 typedef struct 
2390 {
2391   int n;
2392   int log2n;
2393   
2394   float *trig;
2395   int   *bitrev;
2396
2397   float scale;
2398 } mdct_lookup;
2399
2400 extern void mdct_init(mdct_lookup *lookup, int n);
2401 extern void mdct_clear(mdct_lookup *l);
2402 extern void mdct_backward(mdct_lookup *init, float *in, float *out);
2403
2404 mdct_lookup M1,M2;
2405
2406 void inverse_mdct(float *buffer, int n, vorb *f, int blocktype)
2407 {
2408    mdct_lookup *M;
2409    if (M1.n == n) M = &M1;
2410    else if (M2.n == n) M = &M2;
2411    else if (M1.n == 0) { mdct_init(&M1, n); M = &M1; }
2412    else { 
2413       if (M2.n) __asm int 3;
2414       mdct_init(&M2, n);
2415       M = &M2;
2416    }
2417
2418    mdct_backward(M, buffer, buffer);
2419 }
2420 #endif
2421
2422
2423 // the following were split out into separate functions while optimizing;
2424 // they could be pushed back up but eh. __forceinline showed no change;
2425 // they're probably already being inlined.
2426 static void imdct_step3_iter0_loop(int n, float *e, int i_off, int k_off, float *A)
2427 {
2428    float *ee0 = e + i_off;
2429    float *ee2 = ee0 + k_off;
2430    int i;
2431
2432    assert((n & 3) == 0);
2433    for (i=(n>>2); i > 0; --i) {
2434       float k00_20, k01_21;
2435       k00_20  = ee0[ 0] - ee2[ 0];
2436       k01_21  = ee0[-1] - ee2[-1];
2437       ee0[ 0] += ee2[ 0];//ee0[ 0] = ee0[ 0] + ee2[ 0];
2438       ee0[-1] += ee2[-1];//ee0[-1] = ee0[-1] + ee2[-1];
2439       ee2[ 0] = k00_20 * A[0] - k01_21 * A[1];
2440       ee2[-1] = k01_21 * A[0] + k00_20 * A[1];
2441       A += 8;
2442
2443       k00_20  = ee0[-2] - ee2[-2];
2444       k01_21  = ee0[-3] - ee2[-3];
2445       ee0[-2] += ee2[-2];//ee0[-2] = ee0[-2] + ee2[-2];
2446       ee0[-3] += ee2[-3];//ee0[-3] = ee0[-3] + ee2[-3];
2447       ee2[-2] = k00_20 * A[0] - k01_21 * A[1];
2448       ee2[-3] = k01_21 * A[0] + k00_20 * A[1];
2449       A += 8;
2450
2451       k00_20  = ee0[-4] - ee2[-4];
2452       k01_21  = ee0[-5] - ee2[-5];
2453       ee0[-4] += ee2[-4];//ee0[-4] = ee0[-4] + ee2[-4];
2454       ee0[-5] += ee2[-5];//ee0[-5] = ee0[-5] + ee2[-5];
2455       ee2[-4] = k00_20 * A[0] - k01_21 * A[1];
2456       ee2[-5] = k01_21 * A[0] + k00_20 * A[1];
2457       A += 8;
2458
2459       k00_20  = ee0[-6] - ee2[-6];
2460       k01_21  = ee0[-7] - ee2[-7];
2461       ee0[-6] += ee2[-6];//ee0[-6] = ee0[-6] + ee2[-6];
2462       ee0[-7] += ee2[-7];//ee0[-7] = ee0[-7] + ee2[-7];
2463       ee2[-6] = k00_20 * A[0] - k01_21 * A[1];
2464       ee2[-7] = k01_21 * A[0] + k00_20 * A[1];
2465       A += 8;
2466       ee0 -= 8;
2467       ee2 -= 8;
2468    }
2469 }
2470
2471 static void imdct_step3_inner_r_loop(int lim, float *e, int d0, int k_off, float *A, int k1)
2472 {
2473    int i;
2474    float k00_20, k01_21;
2475
2476    float *e0 = e + d0;
2477    float *e2 = e0 + k_off;
2478
2479    for (i=lim >> 2; i > 0; --i) {
2480       k00_20 = e0[-0] - e2[-0];
2481       k01_21 = e0[-1] - e2[-1];
2482       e0[-0] += e2[-0];//e0[-0] = e0[-0] + e2[-0];
2483       e0[-1] += e2[-1];//e0[-1] = e0[-1] + e2[-1];
2484       e2[-0] = (k00_20)*A[0] - (k01_21) * A[1];
2485       e2[-1] = (k01_21)*A[0] + (k00_20) * A[1];
2486
2487       A += k1;
2488
2489       k00_20 = e0[-2] - e2[-2];
2490       k01_21 = e0[-3] - e2[-3];
2491       e0[-2] += e2[-2];//e0[-2] = e0[-2] + e2[-2];
2492       e0[-3] += e2[-3];//e0[-3] = e0[-3] + e2[-3];
2493       e2[-2] = (k00_20)*A[0] - (k01_21) * A[1];
2494       e2[-3] = (k01_21)*A[0] + (k00_20) * A[1];
2495
2496       A += k1;
2497
2498       k00_20 = e0[-4] - e2[-4];
2499       k01_21 = e0[-5] - e2[-5];
2500       e0[-4] += e2[-4];//e0[-4] = e0[-4] + e2[-4];
2501       e0[-5] += e2[-5];//e0[-5] = e0[-5] + e2[-5];
2502       e2[-4] = (k00_20)*A[0] - (k01_21) * A[1];
2503       e2[-5] = (k01_21)*A[0] + (k00_20) * A[1];
2504
2505       A += k1;
2506
2507       k00_20 = e0[-6] - e2[-6];
2508       k01_21 = e0[-7] - e2[-7];
2509       e0[-6] += e2[-6];//e0[-6] = e0[-6] + e2[-6];
2510       e0[-7] += e2[-7];//e0[-7] = e0[-7] + e2[-7];
2511       e2[-6] = (k00_20)*A[0] - (k01_21) * A[1];
2512       e2[-7] = (k01_21)*A[0] + (k00_20) * A[1];
2513
2514       e0 -= 8;
2515       e2 -= 8;
2516
2517       A += k1;
2518    }
2519 }
2520
2521 static void imdct_step3_inner_s_loop(int n, float *e, int i_off, int k_off, float *A, int a_off, int k0)
2522 {
2523    int i;
2524    float A0 = A[0];
2525    float A1 = A[0+1];
2526    float A2 = A[0+a_off];
2527    float A3 = A[0+a_off+1];
2528    float A4 = A[0+a_off*2+0];
2529    float A5 = A[0+a_off*2+1];
2530    float A6 = A[0+a_off*3+0];
2531    float A7 = A[0+a_off*3+1];
2532
2533    float k00,k11;
2534
2535    float *ee0 = e  +i_off;
2536    float *ee2 = ee0+k_off;
2537
2538    for (i=n; i > 0; --i) {
2539       k00     = ee0[ 0] - ee2[ 0];
2540       k11     = ee0[-1] - ee2[-1];
2541       ee0[ 0] =  ee0[ 0] + ee2[ 0];
2542       ee0[-1] =  ee0[-1] + ee2[-1];
2543       ee2[ 0] = (k00) * A0 - (k11) * A1;
2544       ee2[-1] = (k11) * A0 + (k00) * A1;
2545
2546       k00     = ee0[-2] - ee2[-2];
2547       k11     = ee0[-3] - ee2[-3];
2548       ee0[-2] =  ee0[-2] + ee2[-2];
2549       ee0[-3] =  ee0[-3] + ee2[-3];
2550       ee2[-2] = (k00) * A2 - (k11) * A3;
2551       ee2[-3] = (k11) * A2 + (k00) * A3;
2552
2553       k00     = ee0[-4] - ee2[-4];
2554       k11     = ee0[-5] - ee2[-5];
2555       ee0[-4] =  ee0[-4] + ee2[-4];
2556       ee0[-5] =  ee0[-5] + ee2[-5];
2557       ee2[-4] = (k00) * A4 - (k11) * A5;
2558       ee2[-5] = (k11) * A4 + (k00) * A5;
2559
2560       k00     = ee0[-6] - ee2[-6];
2561       k11     = ee0[-7] - ee2[-7];
2562       ee0[-6] =  ee0[-6] + ee2[-6];
2563       ee0[-7] =  ee0[-7] + ee2[-7];
2564       ee2[-6] = (k00) * A6 - (k11) * A7;
2565       ee2[-7] = (k11) * A6 + (k00) * A7;
2566
2567       ee0 -= k0;
2568       ee2 -= k0;
2569    }
2570 }
2571
2572 static __forceinline void iter_54(float *z)
2573 {
2574    float k00,k11,k22,k33;
2575    float y0,y1,y2,y3;
2576
2577    k00  = z[ 0] - z[-4];
2578    y0   = z[ 0] + z[-4];
2579    y2   = z[-2] + z[-6];
2580    k22  = z[-2] - z[-6];
2581
2582    z[-0] = y0 + y2;      // z0 + z4 + z2 + z6
2583    z[-2] = y0 - y2;      // z0 + z4 - z2 - z6
2584
2585    // done with y0,y2
2586
2587    k33  = z[-3] - z[-7];
2588
2589    z[-4] = k00 + k33;    // z0 - z4 + z3 - z7
2590    z[-6] = k00 - k33;    // z0 - z4 - z3 + z7
2591
2592    // done with k33
2593
2594    k11  = z[-1] - z[-5];
2595    y1   = z[-1] + z[-5];
2596    y3   = z[-3] + z[-7];
2597
2598    z[-1] = y1 + y3;      // z1 + z5 + z3 + z7
2599    z[-3] = y1 - y3;      // z1 + z5 - z3 - z7
2600    z[-5] = k11 - k22;    // z1 - z5 + z2 - z6
2601    z[-7] = k11 + k22;    // z1 - z5 - z2 + z6
2602 }
2603
2604 static void imdct_step3_inner_s_loop_ld654(int n, float *e, int i_off, float *A, int base_n)
2605 {
2606    int k_off = -8;
2607    int a_off = base_n >> 3;
2608    float A2 = A[0+a_off];
2609    float *z = e + i_off;
2610    float *base = z - 16 * n;
2611
2612    while (z > base) {
2613       float k00,k11;
2614
2615       k00   = z[-0] - z[-8];
2616       k11   = z[-1] - z[-9];
2617       z[-0] = z[-0] + z[-8];
2618       z[-1] = z[-1] + z[-9];
2619       z[-8] =  k00;
2620       z[-9] =  k11 ;
2621
2622       k00    = z[ -2] - z[-10];
2623       k11    = z[ -3] - z[-11];
2624       z[ -2] = z[ -2] + z[-10];
2625       z[ -3] = z[ -3] + z[-11];
2626       z[-10] = (k00+k11) * A2;
2627       z[-11] = (k11-k00) * A2;
2628
2629       k00    = z[-12] - z[ -4];  // reverse to avoid a unary negation
2630       k11    = z[ -5] - z[-13];
2631       z[ -4] = z[ -4] + z[-12];
2632       z[ -5] = z[ -5] + z[-13];
2633       z[-12] = k11;
2634       z[-13] = k00;
2635
2636       k00    = z[-14] - z[ -6];  // reverse to avoid a unary negation
2637       k11    = z[ -7] - z[-15];
2638       z[ -6] = z[ -6] + z[-14];
2639       z[ -7] = z[ -7] + z[-15];
2640       z[-14] = (k00+k11) * A2;
2641       z[-15] = (k00-k11) * A2;
2642
2643       iter_54(z);
2644       iter_54(z-8);
2645       z -= 16;
2646    }
2647 }
2648
2649 static void inverse_mdct(float *buffer, int n, vorb *f, int blocktype)
2650 {
2651    int n2 = n >> 1, n4 = n >> 2, n8 = n >> 3, l;
2652    int n3_4 = n - n4, ld;
2653    // @OPTIMIZE: reduce register pressure by using fewer variables?
2654    int save_point = temp_alloc_save(f);
2655    float *buf2 = (float *) temp_alloc(f, n2 * sizeof(*buf2));
2656    float *u=NULL,*v=NULL;
2657    // twiddle factors
2658    float *A = f->A[blocktype];
2659
2660    // IMDCT algorithm from "The use of multirate filter banks for coding of high quality digital audio"
2661    // See notes about bugs in that paper in less-optimal implementation 'inverse_mdct_old' after this function.
2662
2663    // kernel from paper
2664
2665
2666    // merged:
2667    //   copy and reflect spectral data
2668    //   step 0
2669
2670    // note that it turns out that the items added together during
2671    // this step are, in fact, being added to themselves (as reflected
2672    // by step 0). inexplicable inefficiency! this became obvious
2673    // once I combined the passes.
2674
2675    // so there's a missing 'times 2' here (for adding X to itself).
2676    // this propogates through linearly to the end, where the numbers
2677    // are 1/2 too small, and need to be compensated for.
2678
2679    {
2680       float *d,*e, *AA, *e_stop;
2681       d = &buf2[n2-2];
2682       AA = A;
2683       e = &buffer[0];
2684       e_stop = &buffer[n2];
2685       while (e != e_stop) {
2686          d[1] = (e[0] * AA[0] - e[2]*AA[1]);
2687          d[0] = (e[0] * AA[1] + e[2]*AA[0]);
2688          d -= 2;
2689          AA += 2;
2690          e += 4;
2691       }
2692
2693       e = &buffer[n2-3];
2694       while (d >= buf2) {
2695          d[1] = (-e[2] * AA[0] - -e[0]*AA[1]);
2696          d[0] = (-e[2] * AA[1] + -e[0]*AA[0]);
2697          d -= 2;
2698          AA += 2;
2699          e -= 4;
2700       }
2701    }
2702
2703    // now we use symbolic names for these, so that we can
2704    // possibly swap their meaning as we change which operations
2705    // are in place
2706
2707    u = buffer;
2708    v = buf2;
2709
2710    // step 2    (paper output is w, now u)
2711    // this could be in place, but the data ends up in the wrong
2712    // place... _somebody_'s got to swap it, so this is nominated
2713    {
2714       float *AA = &A[n2-8];
2715       float *d0,*d1, *e0, *e1;
2716
2717       e0 = &v[n4];
2718       e1 = &v[0];
2719
2720       d0 = &u[n4];
2721       d1 = &u[0];
2722
2723       while (AA >= A) {
2724          float v40_20, v41_21;
2725
2726          v41_21 = e0[1] - e1[1];
2727          v40_20 = e0[0] - e1[0];
2728          d0[1]  = e0[1] + e1[1];
2729          d0[0]  = e0[0] + e1[0];
2730          d1[1]  = v41_21*AA[4] - v40_20*AA[5];
2731          d1[0]  = v40_20*AA[4] + v41_21*AA[5];
2732
2733          v41_21 = e0[3] - e1[3];
2734          v40_20 = e0[2] - e1[2];
2735          d0[3]  = e0[3] + e1[3];
2736          d0[2]  = e0[2] + e1[2];
2737          d1[3]  = v41_21*AA[0] - v40_20*AA[1];
2738          d1[2]  = v40_20*AA[0] + v41_21*AA[1];
2739
2740          AA -= 8;
2741
2742          d0 += 4;
2743          d1 += 4;
2744          e0 += 4;
2745          e1 += 4;
2746       }
2747    }
2748
2749    // step 3
2750    ld = ilog(n) - 1; // ilog is off-by-one from normal definitions
2751
2752    // optimized step 3:
2753
2754    // the original step3 loop can be nested r inside s or s inside r;
2755    // it's written originally as s inside r, but this is dumb when r
2756    // iterates many times, and s few. So I have two copies of it and
2757    // switch between them halfway.
2758
2759    // this is iteration 0 of step 3
2760    imdct_step3_iter0_loop(n >> 4, u, n2-1-n4*0, -(n >> 3), A);
2761    imdct_step3_iter0_loop(n >> 4, u, n2-1-n4*1, -(n >> 3), A);
2762
2763    // this is iteration 1 of step 3
2764    imdct_step3_inner_r_loop(n >> 5, u, n2-1 - n8*0, -(n >> 4), A, 16);
2765    imdct_step3_inner_r_loop(n >> 5, u, n2-1 - n8*1, -(n >> 4), A, 16);
2766    imdct_step3_inner_r_loop(n >> 5, u, n2-1 - n8*2, -(n >> 4), A, 16);
2767    imdct_step3_inner_r_loop(n >> 5, u, n2-1 - n8*3, -(n >> 4), A, 16);
2768
2769    l=2;
2770    for (; l < (ld-3)>>1; ++l) {
2771       int k0 = n >> (l+2), k0_2 = k0>>1;
2772       int lim = 1 << (l+1);
2773       int i;
2774       for (i=0; i < lim; ++i)
2775          imdct_step3_inner_r_loop(n >> (l+4), u, n2-1 - k0*i, -k0_2, A, 1 << (l+3));
2776    }
2777
2778    for (; l < ld-6; ++l) {
2779       int k0 = n >> (l+2), k1 = 1 << (l+3), k0_2 = k0>>1;
2780       int rlim = n >> (l+6), r;
2781       int lim = 1 << (l+1);
2782       int i_off;
2783       float *A0 = A;
2784       i_off = n2-1;
2785       for (r=rlim; r > 0; --r) {
2786          imdct_step3_inner_s_loop(lim, u, i_off, -k0_2, A0, k1, k0);
2787          A0 += k1*4;
2788          i_off -= 8;
2789       }
2790    }
2791
2792    // iterations with count:
2793    //   ld-6,-5,-4 all interleaved together
2794    //       the big win comes from getting rid of needless flops
2795    //         due to the constants on pass 5 & 4 being all 1 and 0;
2796    //       combining them to be simultaneous to improve cache made little difference
2797    imdct_step3_inner_s_loop_ld654(n >> 5, u, n2-1, A, n);
2798
2799    // output is u
2800
2801    // step 4, 5, and 6
2802    // cannot be in-place because of step 5
2803    {
2804       uint16 *bitrev = f->bit_reverse[blocktype];
2805       // weirdly, I'd have thought reading sequentially and writing
2806       // erratically would have been better than vice-versa, but in
2807       // fact that's not what my testing showed. (That is, with
2808       // j = bitreverse(i), do you read i and write j, or read j and write i.)
2809
2810       float *d0 = &v[n4-4];
2811       float *d1 = &v[n2-4];
2812       while (d0 >= v) {
2813          int k4;
2814
2815          k4 = bitrev[0];
2816          d1[3] = u[k4+0];
2817          d1[2] = u[k4+1];
2818          d0[3] = u[k4+2];
2819          d0[2] = u[k4+3];
2820
2821          k4 = bitrev[1];
2822          d1[1] = u[k4+0];
2823          d1[0] = u[k4+1];
2824          d0[1] = u[k4+2];
2825          d0[0] = u[k4+3];
2826          
2827          d0 -= 4;
2828          d1 -= 4;
2829          bitrev += 2;
2830       }
2831    }
2832    // (paper output is u, now v)
2833
2834
2835    // data must be in buf2
2836    assert(v == buf2);
2837
2838    // step 7   (paper output is v, now v)
2839    // this is now in place
2840    {
2841       float *C = f->C[blocktype];
2842       float *d, *e;
2843
2844       d = v;
2845       e = v + n2 - 4;
2846
2847       while (d < e) {
2848          float a02,a11,b0,b1,b2,b3;
2849
2850          a02 = d[0] - e[2];
2851          a11 = d[1] + e[3];
2852
2853          b0 = C[1]*a02 + C[0]*a11;
2854          b1 = C[1]*a11 - C[0]*a02;
2855
2856          b2 = d[0] + e[ 2];
2857          b3 = d[1] - e[ 3];
2858
2859          d[0] = b2 + b0;
2860          d[1] = b3 + b1;
2861          e[2] = b2 - b0;
2862          e[3] = b1 - b3;
2863
2864          a02 = d[2] - e[0];
2865          a11 = d[3] + e[1];
2866
2867          b0 = C[3]*a02 + C[2]*a11;
2868          b1 = C[3]*a11 - C[2]*a02;
2869
2870          b2 = d[2] + e[ 0];
2871          b3 = d[3] - e[ 1];
2872
2873          d[2] = b2 + b0;
2874          d[3] = b3 + b1;
2875          e[0] = b2 - b0;
2876          e[1] = b1 - b3;
2877
2878          C += 4;
2879          d += 4;
2880          e -= 4;
2881       }
2882    }
2883
2884    // data must be in buf2
2885
2886
2887    // step 8+decode   (paper output is X, now buffer)
2888    // this generates pairs of data a la 8 and pushes them directly through
2889    // the decode kernel (pushing rather than pulling) to avoid having
2890    // to make another pass later
2891
2892    // this cannot POSSIBLY be in place, so we refer to the buffers directly
2893
2894    {
2895       float *d0,*d1,*d2,*d3;
2896
2897       float *B = f->B[blocktype] + n2 - 8;
2898       float *e = buf2 + n2 - 8;
2899       d0 = &buffer[0];
2900       d1 = &buffer[n2-4];
2901       d2 = &buffer[n2];
2902       d3 = &buffer[n-4];
2903       while (e >= v) {
2904          float p0,p1,p2,p3;
2905
2906          p3 =  e[6]*B[7] - e[7]*B[6];
2907          p2 = -e[6]*B[6] - e[7]*B[7]; 
2908
2909          d0[0] =   p3;
2910          d1[3] = - p3;
2911          d2[0] =   p2;
2912          d3[3] =   p2;
2913
2914          p1 =  e[4]*B[5] - e[5]*B[4];
2915          p0 = -e[4]*B[4] - e[5]*B[5]; 
2916
2917          d0[1] =   p1;
2918          d1[2] = - p1;
2919          d2[1] =   p0;
2920          d3[2] =   p0;
2921
2922          p3 =  e[2]*B[3] - e[3]*B[2];
2923          p2 = -e[2]*B[2] - e[3]*B[3]; 
2924
2925          d0[2] =   p3;
2926          d1[1] = - p3;
2927          d2[2] =   p2;
2928          d3[1] =   p2;
2929
2930          p1 =  e[0]*B[1] - e[1]*B[0];
2931          p0 = -e[0]*B[0] - e[1]*B[1]; 
2932
2933          d0[3] =   p1;
2934          d1[0] = - p1;
2935          d2[3] =   p0;
2936          d3[0] =   p0;
2937
2938          B -= 8;
2939          e -= 8;
2940          d0 += 4;
2941          d2 += 4;
2942          d1 -= 4;
2943          d3 -= 4;
2944       }
2945    }
2946
2947    temp_alloc_restore(f,save_point);
2948 }
2949
2950 #if 0
2951 // this is the original version of the above code, if you want to optimize it from scratch
2952 void inverse_mdct_naive(float *buffer, int n)
2953 {
2954    float s;
2955    float A[1 << 12], B[1 << 12], C[1 << 11];
2956    int i,k,k2,k4, n2 = n >> 1, n4 = n >> 2, n8 = n >> 3, l;
2957    int n3_4 = n - n4, ld;
2958    // how can they claim this only uses N words?!
2959    // oh, because they're only used sparsely, whoops
2960    float u[1 << 13], X[1 << 13], v[1 << 13], w[1 << 13];
2961    // set up twiddle factors
2962
2963    for (k=k2=0; k < n4; ++k,k2+=2) {
2964       A[k2  ] = (float)  cos(4*k*M_PI/n);
2965       A[k2+1] = (float) -sin(4*k*M_PI/n);
2966       B[k2  ] = (float)  cos((k2+1)*M_PI/n/2);
2967       B[k2+1] = (float)  sin((k2+1)*M_PI/n/2);
2968    }
2969    for (k=k2=0; k < n8; ++k,k2+=2) {
2970       C[k2  ] = (float)  cos(2*(k2+1)*M_PI/n);
2971       C[k2+1] = (float) -sin(2*(k2+1)*M_PI/n);
2972    }
2973
2974    // IMDCT algorithm from "The use of multirate filter banks for coding of high quality digital audio"
2975    // Note there are bugs in that pseudocode, presumably due to them attempting
2976    // to rename the arrays nicely rather than representing the way their actual
2977    // implementation bounces buffers back and forth. As a result, even in the
2978    // "some formulars corrected" version, a direct implementation fails. These
2979    // are noted below as "paper bug".
2980
2981    // copy and reflect spectral data
2982    for (k=0; k < n2; ++k) u[k] = buffer[k];
2983    for (   ; k < n ; ++k) u[k] = -buffer[n - k - 1];
2984    // kernel from paper
2985    // step 1
2986    for (k=k2=k4=0; k < n4; k+=1, k2+=2, k4+=4) {
2987       v[n-k4-1] = (u[k4] - u[n-k4-1]) * A[k2]   - (u[k4+2] - u[n-k4-3])*A[k2+1];
2988       v[n-k4-3] = (u[k4] - u[n-k4-1]) * A[k2+1] + (u[k4+2] - u[n-k4-3])*A[k2];
2989    }
2990    // step 2
2991    for (k=k4=0; k < n8; k+=1, k4+=4) {
2992       w[n2+3+k4] = v[n2+3+k4] + v[k4+3];
2993       w[n2+1+k4] = v[n2+1+k4] + v[k4+1];
2994       w[k4+3]    = (v[n2+3+k4] - v[k4+3])*A[n2-4-k4] - (v[n2+1+k4]-v[k4+1])*A[n2-3-k4];
2995       w[k4+1]    = (v[n2+1+k4] - v[k4+1])*A[n2-4-k4] + (v[n2+3+k4]-v[k4+3])*A[n2-3-k4];
2996    }
2997    // step 3
2998    ld = ilog(n) - 1; // ilog is off-by-one from normal definitions
2999    for (l=0; l < ld-3; ++l) {
3000       int k0 = n >> (l+2), k1 = 1 << (l+3);
3001       int rlim = n >> (l+4), r4, r;
3002       int s2lim = 1 << (l+2), s2;
3003       for (r=r4=0; r < rlim; r4+=4,++r) {
3004          for (s2=0; s2 < s2lim; s2+=2) {
3005             u[n-1-k0*s2-r4] = w[n-1-k0*s2-r4] + w[n-1-k0*(s2+1)-r4];
3006             u[n-3-k0*s2-r4] = w[n-3-k0*s2-r4] + w[n-3-k0*(s2+1)-r4];
3007             u[n-1-k0*(s2+1)-r4] = (w[n-1-k0*s2-r4] - w[n-1-k0*(s2+1)-r4]) * A[r*k1]
3008                                 - (w[n-3-k0*s2-r4] - w[n-3-k0*(s2+1)-r4]) * A[r*k1+1];
3009             u[n-3-k0*(s2+1)-r4] = (w[n-3-k0*s2-r4] - w[n-3-k0*(s2+1)-r4]) * A[r*k1]
3010                                 + (w[n-1-k0*s2-r4] - w[n-1-k0*(s2+1)-r4]) * A[r*k1+1];
3011          }
3012       }
3013       if (l+1 < ld-3) {
3014          // paper bug: ping-ponging of u&w here is omitted
3015          memcpy(w, u, sizeof(u));
3016       }
3017    }
3018
3019    // step 4
3020    for (i=0; i < n8; ++i) {
3021       int j = bit_reverse(i) >> (32-ld+3);
3022       assert(j < n8);
3023       if (i == j) {
3024          // paper bug: original code probably swapped in place; if copying,
3025          //            need to directly copy in this case
3026          int i8 = i << 3;
3027          v[i8+1] = u[i8+1];
3028          v[i8+3] = u[i8+3];
3029          v[i8+5] = u[i8+5];
3030          v[i8+7] = u[i8+7];
3031       } else if (i < j) {
3032          int i8 = i << 3, j8 = j << 3;
3033          v[j8+1] = u[i8+1], v[i8+1] = u[j8 + 1];
3034          v[j8+3] = u[i8+3], v[i8+3] = u[j8 + 3];
3035          v[j8+5] = u[i8+5], v[i8+5] = u[j8 + 5];
3036          v[j8+7] = u[i8+7], v[i8+7] = u[j8 + 7];
3037       }
3038    }
3039    // step 5
3040    for (k=0; k < n2; ++k) {
3041       w[k] = v[k*2+1];
3042    }
3043    // step 6
3044    for (k=k2=k4=0; k < n8; ++k, k2 += 2, k4 += 4) {
3045       u[n-1-k2] = w[k4];
3046       u[n-2-k2] = w[k4+1];
3047       u[n3_4 - 1 - k2] = w[k4+2];
3048       u[n3_4 - 2 - k2] = w[k4+3];
3049    }
3050    // step 7
3051    for (k=k2=0; k < n8; ++k, k2 += 2) {
3052       v[n2 + k2 ] = ( u[n2 + k2] + u[n-2-k2] + C[k2+1]*(u[n2+k2]-u[n-2-k2]) + C[k2]*(u[n2+k2+1]+u[n-2-k2+1]))/2;
3053       v[n-2 - k2] = ( u[n2 + k2] + u[n-2-k2] - C[k2+1]*(u[n2+k2]-u[n-2-k2]) - C[k2]*(u[n2+k2+1]+u[n-2-k2+1]))/2;
3054       v[n2+1+ k2] = ( u[n2+1+k2] - u[n-1-k2] + C[k2+1]*(u[n2+1+k2]+u[n-1-k2]) - C[k2]*(u[n2+k2]-u[n-2-k2]))/2;
3055       v[n-1 - k2] = (-u[n2+1+k2] + u[n-1-k2] + C[k2+1]*(u[n2+1+k2]+u[n-1-k2]) - C[k2]*(u[n2+k2]-u[n-2-k2]))/2;
3056    }
3057    // step 8
3058    for (k=k2=0; k < n4; ++k,k2 += 2) {
3059       X[k]      = v[k2+n2]*B[k2  ] + v[k2+1+n2]*B[k2+1];
3060       X[n2-1-k] = v[k2+n2]*B[k2+1] - v[k2+1+n2]*B[k2  ];
3061    }
3062
3063    // decode kernel to output
3064    // determined the following value experimentally
3065    // (by first figuring out what made inverse_mdct_slow work); then matching that here
3066    // (probably vorbis encoder premultiplies by n or n/2, to save it on the decoder?)
3067    s = 0.5; // theoretically would be n4
3068
3069    // [[[ note! the s value of 0.5 is compensated for by the B[] in the current code,
3070    //     so it needs to use the "old" B values to behave correctly, or else
3071    //     set s to 1.0 ]]]
3072    for (i=0; i < n4  ; ++i) buffer[i] = s * X[i+n4];
3073    for (   ; i < n3_4; ++i) buffer[i] = -s * X[n3_4 - i - 1];
3074    for (   ; i < n   ; ++i) buffer[i] = -s * X[i - n3_4];
3075 }
3076 #endif
3077
3078 static float *get_window(vorb *f, int len)
3079 {
3080    len <<= 1;
3081    if (len == f->blocksize_0) return f->window[0];
3082    if (len == f->blocksize_1) return f->window[1];
3083    assert(0);
3084    return NULL;
3085 }
3086
3087 #ifndef STB_VORBIS_NO_DEFER_FLOOR
3088 typedef int16 YTYPE;
3089 #else
3090 typedef int YTYPE;
3091 #endif
3092 static int do_floor(vorb *f, Mapping *map, int i, int n, float *target, YTYPE *finalY, uint8 *step2_flag)
3093 {
3094    int n2 = n >> 1;
3095    int s = map->chan[i].mux, floor;
3096    floor = map->submap_floor[s];
3097    if (f->floor_types[floor] == 0) {
3098       return error(f, VORBIS_invalid_stream);
3099    } else {
3100       Floor1 *g = &f->floor_config[floor].floor1;
3101       int j,q;
3102       int lx = 0, ly = finalY[0] * g->floor1_multiplier;
3103       for (q=1; q < g->values; ++q) {
3104          j = g->sorted_order[q];
3105          #ifndef STB_VORBIS_NO_DEFER_FLOOR
3106          if (finalY[j] >= 0)
3107          #else
3108          if (step2_flag[j])
3109          #endif
3110          {
3111             int hy = finalY[j] * g->floor1_multiplier;
3112             int hx = g->Xlist[j];
3113             draw_line(target, lx,ly, hx,hy, n2);
3114             lx = hx, ly = hy;
3115          }
3116       }
3117       if (lx < n2)
3118          // optimization of: draw_line(target, lx,ly, n,ly, n2);
3119          for (j=lx; j < n2; ++j)
3120             LINE_OP(target[j], inverse_db_table[ly]);
3121    }
3122    return TRUE;
3123 }
3124
3125 static int vorbis_decode_initial(vorb *f, int *p_left_start, int *p_left_end, int *p_right_start, int *p_right_end, int *mode)
3126 {
3127    Mode *m;
3128    int i, n, prev, next, window_center;
3129    f->channel_buffer_start = f->channel_buffer_end = 0;
3130
3131   retry:
3132    if (f->eof) return FALSE;
3133    if (!maybe_start_packet(f))
3134       return FALSE;
3135    // check packet type
3136    if (get_bits(f,1) != 0) {
3137       if (IS_PUSH_MODE(f))
3138          return error(f,VORBIS_bad_packet_type);
3139       while (EOP != get8_packet(f));
3140       goto retry;
3141    }
3142
3143    if (f->alloc.alloc_buffer)
3144       assert(f->alloc.alloc_buffer_length_in_bytes == f->temp_offset);
3145
3146    i = get_bits(f, ilog(f->mode_count-1));
3147    if (i == EOP) return FALSE;
3148    if (i >= f->mode_count) return FALSE;
3149    *mode = i;
3150    m = f->mode_config + i;
3151    if (m->blockflag) {
3152       n = f->blocksize_1;
3153       prev = get_bits(f,1);
3154       next = get_bits(f,1);
3155    } else {
3156       prev = next = 0;
3157       n = f->blocksize_0;
3158    }
3159
3160 // WINDOWING
3161
3162    window_center = n >> 1;
3163    if (m->blockflag && !prev) {