Initial commit
[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    free(x);
2367 }
2368
2369 void inverse_mdct_slow(float *buffer, int n, vorb *f, int blocktype)
2370 {
2371    int i, n4 = n >> 2, n2 = n >> 1, n3_4 = n - n4;
2372    float temp[4096];
2373
2374    memcpy(temp, buffer, n2 * sizeof(float));
2375    dct_iv_slow(temp, n2);  // returns -c'-d, a-b'
2376
2377    for (i=0; i < n4  ; ++i) buffer[i] = temp[i+n4];            // a-b'
2378    for (   ; i < n3_4; ++i) buffer[i] = -temp[n3_4 - i - 1];   // b-a', c+d'
2379    for (   ; i < n   ; ++i) buffer[i] = -temp[i - n3_4];       // c'+d
2380 }
2381 #endif
2382
2383 #ifndef LIBVORBIS_MDCT
2384 #define LIBVORBIS_MDCT 0
2385 #endif
2386
2387 #if LIBVORBIS_MDCT
2388 // directly call the vorbis MDCT using an interface documented
2389 // by Jeff Roberts... useful for performance comparison
2390 typedef struct 
2391 {
2392   int n;
2393   int log2n;
2394   
2395   float *trig;
2396   int   *bitrev;
2397
2398   float scale;
2399 } mdct_lookup;
2400
2401 extern void mdct_init(mdct_lookup *lookup, int n);
2402 extern void mdct_clear(mdct_lookup *l);
2403 extern void mdct_backward(mdct_lookup *init, float *in, float *out);
2404
2405 mdct_lookup M1,M2;
2406
2407 void inverse_mdct(float *buffer, int n, vorb *f, int blocktype)
2408 {
2409    mdct_lookup *M;
2410    if (M1.n == n) M = &M1;
2411    else if (M2.n == n) M = &M2;
2412    else if (M1.n == 0) { mdct_init(&M1, n); M = &M1; }
2413    else { 
2414       if (M2.n) __asm int 3;
2415       mdct_init(&M2, n);
2416       M = &M2;
2417    }
2418
2419    mdct_backward(M, buffer, buffer);
2420 }
2421 #endif
2422
2423
2424 // the following were split out into separate functions while optimizing;
2425 // they could be pushed back up but eh. __forceinline showed no change;
2426 // they're probably already being inlined.
2427 static void imdct_step3_iter0_loop(int n, float *e, int i_off, int k_off, float *A)
2428 {
2429    float *ee0 = e + i_off;
2430    float *ee2 = ee0 + k_off;
2431    int i;
2432
2433    assert((n & 3) == 0);
2434    for (i=(n>>2); i > 0; --i) {
2435       float k00_20, k01_21;
2436       k00_20  = ee0[ 0] - ee2[ 0];
2437       k01_21  = ee0[-1] - ee2[-1];
2438       ee0[ 0] += ee2[ 0];//ee0[ 0] = ee0[ 0] + ee2[ 0];
2439       ee0[-1] += ee2[-1];//ee0[-1] = ee0[-1] + ee2[-1];
2440       ee2[ 0] = k00_20 * A[0] - k01_21 * A[1];
2441       ee2[-1] = k01_21 * A[0] + k00_20 * A[1];
2442       A += 8;
2443
2444       k00_20  = ee0[-2] - ee2[-2];
2445       k01_21  = ee0[-3] - ee2[-3];
2446       ee0[-2] += ee2[-2];//ee0[-2] = ee0[-2] + ee2[-2];
2447       ee0[-3] += ee2[-3];//ee0[-3] = ee0[-3] + ee2[-3];
2448       ee2[-2] = k00_20 * A[0] - k01_21 * A[1];
2449       ee2[-3] = k01_21 * A[0] + k00_20 * A[1];
2450       A += 8;
2451
2452       k00_20  = ee0[-4] - ee2[-4];
2453       k01_21  = ee0[-5] - ee2[-5];
2454       ee0[-4] += ee2[-4];//ee0[-4] = ee0[-4] + ee2[-4];
2455       ee0[-5] += ee2[-5];//ee0[-5] = ee0[-5] + ee2[-5];
2456       ee2[-4] = k00_20 * A[0] - k01_21 * A[1];
2457       ee2[-5] = k01_21 * A[0] + k00_20 * A[1];
2458       A += 8;
2459
2460       k00_20  = ee0[-6] - ee2[-6];
2461       k01_21  = ee0[-7] - ee2[-7];
2462       ee0[-6] += ee2[-6];//ee0[-6] = ee0[-6] + ee2[-6];
2463       ee0[-7] += ee2[-7];//ee0[-7] = ee0[-7] + ee2[-7];
2464       ee2[-6] = k00_20 * A[0] - k01_21 * A[1];
2465       ee2[-7] = k01_21 * A[0] + k00_20 * A[1];
2466       A += 8;
2467       ee0 -= 8;
2468       ee2 -= 8;
2469    }
2470 }
2471
2472 static void imdct_step3_inner_r_loop(int lim, float *e, int d0, int k_off, float *A, int k1)
2473 {
2474    int i;
2475    float k00_20, k01_21;
2476
2477    float *e0 = e + d0;
2478    float *e2 = e0 + k_off;
2479
2480    for (i=lim >> 2; i > 0; --i) {
2481       k00_20 = e0[-0] - e2[-0];
2482       k01_21 = e0[-1] - e2[-1];
2483       e0[-0] += e2[-0];//e0[-0] = e0[-0] + e2[-0];
2484       e0[-1] += e2[-1];//e0[-1] = e0[-1] + e2[-1];
2485       e2[-0] = (k00_20)*A[0] - (k01_21) * A[1];
2486       e2[-1] = (k01_21)*A[0] + (k00_20) * A[1];
2487
2488       A += k1;
2489
2490       k00_20 = e0[-2] - e2[-2];
2491       k01_21 = e0[-3] - e2[-3];
2492       e0[-2] += e2[-2];//e0[-2] = e0[-2] + e2[-2];
2493       e0[-3] += e2[-3];//e0[-3] = e0[-3] + e2[-3];
2494       e2[-2] = (k00_20)*A[0] - (k01_21) * A[1];
2495       e2[-3] = (k01_21)*A[0] + (k00_20) * A[1];
2496
2497       A += k1;
2498
2499       k00_20 = e0[-4] - e2[-4];
2500       k01_21 = e0[-5] - e2[-5];
2501       e0[-4] += e2[-4];//e0[-4] = e0[-4] + e2[-4];
2502       e0[-5] += e2[-5];//e0[-5] = e0[-5] + e2[-5];
2503       e2[-4] = (k00_20)*A[0] - (k01_21) * A[1];
2504       e2[-5] = (k01_21)*A[0] + (k00_20) * A[1];
2505
2506       A += k1;
2507
2508       k00_20 = e0[-6] - e2[-6];
2509       k01_21 = e0[-7] - e2[-7];
2510       e0[-6] += e2[-6];//e0[-6] = e0[-6] + e2[-6];
2511       e0[-7] += e2[-7];//e0[-7] = e0[-7] + e2[-7];
2512       e2[-6] = (k00_20)*A[0] - (k01_21) * A[1];
2513       e2[-7] = (k01_21)*A[0] + (k00_20) * A[1];
2514
2515       e0 -= 8;
2516       e2 -= 8;
2517
2518       A += k1;
2519    }
2520 }
2521
2522 static void imdct_step3_inner_s_loop(int n, float *e, int i_off, int k_off, float *A, int a_off, int k0)
2523 {
2524    int i;
2525    float A0 = A[0];
2526    float A1 = A[0+1];
2527    float A2 = A[0+a_off];
2528    float A3 = A[0+a_off+1];
2529    float A4 = A[0+a_off*2+0];
2530    float A5 = A[0+a_off*2+1];
2531    float A6 = A[0+a_off*3+0];
2532    float A7 = A[0+a_off*3+1];
2533
2534    float k00,k11;
2535
2536    float *ee0 = e  +i_off;
2537    float *ee2 = ee0+k_off;
2538
2539    for (i=n; i > 0; --i) {
2540       k00     = ee0[ 0] - ee2[ 0];
2541       k11     = ee0[-1] - ee2[-1];
2542       ee0[ 0] =  ee0[ 0] + ee2[ 0];
2543       ee0[-1] =  ee0[-1] + ee2[-1];
2544       ee2[ 0] = (k00) * A0 - (k11) * A1;
2545       ee2[-1] = (k11) * A0 + (k00) * A1;
2546
2547       k00     = ee0[-2] - ee2[-2];
2548       k11     = ee0[-3] - ee2[-3];
2549       ee0[-2] =  ee0[-2] + ee2[-2];
2550       ee0[-3] =  ee0[-3] + ee2[-3];
2551       ee2[-2] = (k00) * A2 - (k11) * A3;
2552       ee2[-3] = (k11) * A2 + (k00) * A3;
2553
2554       k00     = ee0[-4] - ee2[-4];
2555       k11     = ee0[-5] - ee2[-5];
2556       ee0[-4] =  ee0[-4] + ee2[-4];
2557       ee0[-5] =  ee0[-5] + ee2[-5];
2558       ee2[-4] = (k00) * A4 - (k11) * A5;
2559       ee2[-5] = (k11) * A4 + (k00) * A5;
2560
2561       k00     = ee0[-6] - ee2[-6];
2562       k11     = ee0[-7] - ee2[-7];
2563       ee0[-6] =  ee0[-6] + ee2[-6];
2564       ee0[-7] =  ee0[-7] + ee2[-7];
2565       ee2[-6] = (k00) * A6 - (k11) * A7;
2566       ee2[-7] = (k11) * A6 + (k00) * A7;
2567
2568       ee0 -= k0;
2569       ee2 -= k0;
2570    }
2571 }
2572
2573 static __forceinline void iter_54(float *z)
2574 {
2575    float k00,k11,k22,k33;
2576    float y0,y1,y2,y3;
2577
2578    k00  = z[ 0] - z[-4];
2579    y0   = z[ 0] + z[-4];
2580    y2   = z[-2] + z[-6];
2581    k22  = z[-2] - z[-6];
2582
2583    z[-0] = y0 + y2;      // z0 + z4 + z2 + z6
2584    z[-2] = y0 - y2;      // z0 + z4 - z2 - z6
2585
2586    // done with y0,y2
2587
2588    k33  = z[-3] - z[-7];
2589
2590    z[-4] = k00 + k33;    // z0 - z4 + z3 - z7
2591    z[-6] = k00 - k33;    // z0 - z4 - z3 + z7
2592
2593    // done with k33
2594
2595    k11  = z[-1] - z[-5];
2596    y1   = z[-1] + z[-5];
2597    y3   = z[-3] + z[-7];
2598
2599    z[-1] = y1 + y3;      // z1 + z5 + z3 + z7
2600    z[-3] = y1 - y3;      // z1 + z5 - z3 - z7
2601    z[-5] = k11 - k22;    // z1 - z5 + z2 - z6
2602    z[-7] = k11 + k22;    // z1 - z5 - z2 + z6
2603 }
2604
2605 static void imdct_step3_inner_s_loop_ld654(int n, float *e, int i_off, float *A, int base_n)
2606 {
2607    int k_off = -8;
2608    int a_off = base_n >> 3;
2609    float A2 = A[0+a_off];
2610    float *z = e + i_off;
2611    float *base = z - 16 * n;
2612
2613    while (z > base) {
2614       float k00,k11;
2615
2616       k00   = z[-0] - z[-8];
2617       k11   = z[-1] - z[-9];
2618       z[-0] = z[-0] + z[-8];
2619       z[-1] = z[-1] + z[-9];
2620       z[-8] =  k00;
2621       z[-9] =  k11 ;
2622
2623       k00    = z[ -2] - z[-10];
2624       k11    = z[ -3] - z[-11];
2625       z[ -2] = z[ -2] + z[-10];
2626       z[ -3] = z[ -3] + z[-11];
2627       z[-10] = (k00+k11) * A2;
2628       z[-11] = (k11-k00) * A2;
2629
2630       k00    = z[-12] - z[ -4];  // reverse to avoid a unary negation
2631       k11    = z[ -5] - z[-13];
2632       z[ -4] = z[ -4] + z[-12];
2633       z[ -5] = z[ -5] + z[-13];
2634       z[-12] = k11;
2635       z[-13] = k00;
2636
2637       k00    = z[-14] - z[ -6];  // reverse to avoid a unary negation
2638       k11    = z[ -7] - z[-15];
2639       z[ -6] = z[ -6] + z[-14];
2640       z[ -7] = z[ -7] + z[-15];
2641       z[-14] = (k00+k11) * A2;
2642       z[-15] = (k00-k11) * A2;
2643
2644       iter_54(z);
2645       iter_54(z-8);
2646       z -= 16;
2647    }
2648 }
2649
2650 static void inverse_mdct(float *buffer, int n, vorb *f, int blocktype)
2651 {
2652    int n2 = n >> 1, n4 = n >> 2, n8 = n >> 3, l;
2653    int n3_4 = n - n4, ld;
2654    // @OPTIMIZE: reduce register pressure by using fewer variables?
2655    int save_point = temp_alloc_save(f);
2656    float *buf2 = (float *) temp_alloc(f, n2 * sizeof(*buf2));
2657    float *u=NULL,*v=NULL;
2658    // twiddle factors
2659    float *A = f->A[blocktype];
2660
2661    // IMDCT algorithm from "The use of multirate filter banks for coding of high quality digital audio"
2662    // See notes about bugs in that paper in less-optimal implementation 'inverse_mdct_old' after this function.
2663
2664    // kernel from paper
2665
2666
2667    // merged:
2668    //   copy and reflect spectral data
2669    //   step 0
2670
2671    // note that it turns out that the items added together during
2672    // this step are, in fact, being added to themselves (as reflected
2673    // by step 0). inexplicable inefficiency! this became obvious
2674    // once I combined the passes.
2675
2676    // so there's a missing 'times 2' here (for adding X to itself).
2677    // this propogates through linearly to the end, where the numbers
2678    // are 1/2 too small, and need to be compensated for.
2679
2680    {
2681       float *d,*e, *AA, *e_stop;
2682       d = &buf2[n2-2];
2683       AA = A;
2684       e = &buffer[0];
2685       e_stop = &buffer[n2];
2686       while (e != e_stop) {
2687          d[1] = (e[0] * AA[0] - e[2]*AA[1]);
2688          d[0] = (e[0] * AA[1] + e[2]*AA[0]);
2689          d -= 2;
2690          AA += 2;
2691          e += 4;
2692       }
2693
2694       e = &buffer[n2-3];
2695       while (d >= buf2) {
2696          d[1] = (-e[2] * AA[0] - -e[0]*AA[1]);
2697          d[0] = (-e[2] * AA[1] + -e[0]*AA[0]);
2698          d -= 2;
2699          AA += 2;
2700          e -= 4;
2701       }
2702    }
2703
2704    // now we use symbolic names for these, so that we can
2705    // possibly swap their meaning as we change which operations
2706    // are in place
2707
2708    u = buffer;
2709    v = buf2;
2710
2711    // step 2    (paper output is w, now u)
2712    // this could be in place, but the data ends up in the wrong
2713    // place... _somebody_'s got to swap it, so this is nominated
2714    {
2715       float *AA = &A[n2-8];
2716       float *d0,*d1, *e0, *e1;
2717
2718       e0 = &v[n4];
2719       e1 = &v[0];
2720
2721       d0 = &u[n4];
2722       d1 = &u[0];
2723
2724       while (AA >= A) {
2725          float v40_20, v41_21;
2726
2727          v41_21 = e0[1] - e1[1];
2728          v40_20 = e0[0] - e1[0];
2729          d0[1]  = e0[1] + e1[1];
2730          d0[0]  = e0[0] + e1[0];
2731          d1[1]  = v41_21*AA[4] - v40_20*AA[5];
2732          d1[0]  = v40_20*AA[4] + v41_21*AA[5];
2733
2734          v41_21 = e0[3] - e1[3];
2735          v40_20 = e0[2] - e1[2];
2736          d0[3]  = e0[3] + e1[3];
2737          d0[2]  = e0[2] + e1[2];
2738          d1[3]  = v41_21*AA[0] - v40_20*AA[1];
2739          d1[2]  = v40_20*AA[0] + v41_21*AA[1];
2740
2741          AA -= 8;
2742
2743          d0 += 4;
2744          d1 += 4;
2745          e0 += 4;
2746          e1 += 4;
2747       }
2748    }
2749
2750    // step 3
2751    ld = ilog(n) - 1; // ilog is off-by-one from normal definitions
2752
2753    // optimized step 3:
2754
2755    // the original step3 loop can be nested r inside s or s inside r;
2756    // it's written originally as s inside r, but this is dumb when r
2757    // iterates many times, and s few. So I have two copies of it and
2758    // switch between them halfway.
2759
2760    // this is iteration 0 of step 3
2761    imdct_step3_iter0_loop(n >> 4, u, n2-1-n4*0, -(n >> 3), A);
2762    imdct_step3_iter0_loop(n >> 4, u, n2-1-n4*1, -(n >> 3), A);
2763
2764    // this is iteration 1 of step 3
2765    imdct_step3_inner_r_loop(n >> 5, u, n2-1 - n8*0, -(n >> 4), A, 16);
2766    imdct_step3_inner_r_loop(n >> 5, u, n2-1 - n8*1, -(n >> 4), A, 16);
2767    imdct_step3_inner_r_loop(n >> 5, u, n2-1 - n8*2, -(n >> 4), A, 16);
2768    imdct_step3_inner_r_loop(n >> 5, u, n2-1 - n8*3, -(n >> 4), A, 16);
2769
2770    l=2;
2771    for (; l < (ld-3)>>1; ++l) {
2772       int k0 = n >> (l+2), k0_2 = k0>>1;
2773       int lim = 1 << (l+1);
2774       int i;
2775       for (i=0; i < lim; ++i)
2776          imdct_step3_inner_r_loop(n >> (l+4), u, n2-1 - k0*i, -k0_2, A, 1 << (l+3));
2777    }
2778
2779    for (; l < ld-6; ++l) {
2780       int k0 = n >> (l+2), k1 = 1 << (l+3), k0_2 = k0>>1;
2781       int rlim = n >> (l+6), r;
2782       int lim = 1 << (l+1);
2783       int i_off;
2784       float *A0 = A;
2785       i_off = n2-1;
2786       for (r=rlim; r > 0; --r) {
2787          imdct_step3_inner_s_loop(lim, u, i_off, -k0_2, A0, k1, k0);
2788          A0 += k1*4;
2789          i_off -= 8;
2790       }
2791    }
2792
2793    // iterations with count:
2794    //   ld-6,-5,-4 all interleaved together
2795    //       the big win comes from getting rid of needless flops
2796    //         due to the constants on pass 5 & 4 being all 1 and 0;
2797    //       combining them to be simultaneous to improve cache made little difference
2798    imdct_step3_inner_s_loop_ld654(n >> 5, u, n2-1, A, n);
2799
2800    // output is u
2801
2802    // step 4, 5, and 6
2803    // cannot be in-place because of step 5
2804    {
2805       uint16 *bitrev = f->bit_reverse[blocktype];
2806       // weirdly, I'd have thought reading sequentially and writing
2807       // erratically would have been better than vice-versa, but in
2808       // fact that's not what my testing showed. (That is, with
2809       // j = bitreverse(i), do you read i and write j, or read j and write i.)
2810
2811       float *d0 = &v[n4-4];
2812       float *d1 = &v[n2-4];
2813       while (d0 >= v) {
2814          int k4;
2815
2816          k4 = bitrev[0];
2817          d1[3] = u[k4+0];
2818          d1[2] = u[k4+1];
2819          d0[3] = u[k4+2];
2820          d0[2] = u[k4+3];
2821
2822          k4 = bitrev[1];
2823          d1[1] = u[k4+0];
2824          d1[0] = u[k4+1];
2825          d0[1] = u[k4+2];
2826          d0[0] = u[k4+3];
2827          
2828          d0 -= 4;
2829          d1 -= 4;
2830          bitrev += 2;
2831       }
2832    }
2833    // (paper output is u, now v)
2834
2835
2836    // data must be in buf2
2837    assert(v == buf2);
2838
2839    // step 7   (paper output is v, now v)
2840    // this is now in place
2841    {
2842       float *C = f->C[blocktype];
2843       float *d, *e;
2844
2845       d = v;
2846       e = v + n2 - 4;
2847
2848       while (d < e) {
2849          float a02,a11,b0,b1,b2,b3;
2850
2851          a02 = d[0] - e[2];
2852          a11 = d[1] + e[3];
2853
2854          b0 = C[1]*a02 + C[0]*a11;
2855          b1 = C[1]*a11 - C[0]*a02;
2856
2857          b2 = d[0] + e[ 2];
2858          b3 = d[1] - e[ 3];
2859
2860          d[0] = b2 + b0;
2861          d[1] = b3 + b1;
2862          e[2] = b2 - b0;
2863          e[3] = b1 - b3;
2864
2865          a02 = d[2] - e[0];
2866          a11 = d[3] + e[1];
2867
2868          b0 = C[3]*a02 + C[2]*a11;
2869          b1 = C[3]*a11 - C[2]*a02;
2870
2871          b2 = d[2] + e[ 0];
2872          b3 = d[3] - e[ 1];
2873
2874          d[2] = b2 + b0;
2875          d[3] = b3 + b1;
2876          e[0] = b2 - b0;
2877          e[1] = b1 - b3;
2878
2879          C += 4;
2880          d += 4;
2881          e -= 4;
2882       }
2883    }
2884
2885    // data must be in buf2
2886
2887
2888    // step 8+decode   (paper output is X, now buffer)
2889    // this generates pairs of data a la 8 and pushes them directly through
2890    // the decode kernel (pushing rather than pulling) to avoid having
2891    // to make another pass later
2892
2893    // this cannot POSSIBLY be in place, so we refer to the buffers directly
2894
2895    {
2896       float *d0,*d1,*d2,*d3;
2897
2898       float *B = f->B[blocktype] + n2 - 8;
2899       float *e = buf2 + n2 - 8;
2900       d0 = &buffer[0];
2901       d1 = &buffer[n2-4];
2902       d2 = &buffer[n2];
2903       d3 = &buffer[n-4];
2904       while (e >= v) {
2905          float p0,p1,p2,p3;
2906
2907          p3 =  e[6]*B[7] - e[7]*B[6];
2908          p2 = -e[6]*B[6] - e[7]*B[7]; 
2909
2910          d0[0] =   p3;
2911          d1[3] = - p3;
2912          d2[0] =   p2;
2913          d3[3] =   p2;
2914
2915          p1 =  e[4]*B[5] - e[5]*B[4];
2916          p0 = -e[4]*B[4] - e[5]*B[5]; 
2917
2918          d0[1] =   p1;
2919          d1[2] = - p1;
2920          d2[1] =   p0;
2921          d3[2] =   p0;
2922
2923          p3 =  e[2]*B[3] - e[3]*B[2];
2924          p2 = -e[2]*B[2] - e[3]*B[3]; 
2925
2926          d0[2] =   p3;
2927          d1[1] = - p3;
2928          d2[2] =   p2;
2929          d3[1] =   p2;
2930
2931          p1 =  e[0]*B[1] - e[1]*B[0];
2932          p0 = -e[0]*B[0] - e[1]*B[1]; 
2933
2934          d0[3] =   p1;
2935          d1[0] = - p1;
2936          d2[3] =   p0;
2937          d3[0] =   p0;
2938
2939          B -= 8;
2940          e -= 8;
2941          d0 += 4;
2942          d2 += 4;
2943          d1 -= 4;
2944          d3 -= 4;
2945       }
2946    }
2947
2948    temp_alloc_restore(f,save_point);
2949 }
2950
2951 #if 0
2952 // this is the original version of the above code, if you want to optimize it from scratch
2953 void inverse_mdct_naive(float *buffer, int n)
2954 {
2955    float s;
2956    float A[1 << 12], B[1 << 12], C[1 << 11];
2957    int i,k,k2,k4, n2 = n >> 1, n4 = n >> 2, n8 = n >> 3, l;
2958    int n3_4 = n - n4, ld;
2959    // how can they claim this only uses N words?!
2960    // oh, because they're only used sparsely, whoops
2961    float u[1 << 13], X[1 << 13], v[1 << 13], w[1 << 13];
2962    // set up twiddle factors
2963
2964    for (k=k2=0; k < n4; ++k,k2+=2) {
2965       A[k2  ] = (float)  cos(4*k*M_PI/n);
2966       A[k2+1] = (float) -sin(4*k*M_PI/n);
2967       B[k2  ] = (float)  cos((k2+1)*M_PI/n/2);
2968       B[k2+1] = (float)  sin((k2+1)*M_PI/n/2);
2969    }
2970    for (k=k2=0; k < n8; ++k,k2+=2) {
2971       C[k2  ] = (float)  cos(2*(k2+1)*M_PI/n);
2972       C[k2+1] = (float) -sin(2*(k2+1)*M_PI/n);
2973    }
2974
2975    // IMDCT algorithm from "The use of multirate filter banks for coding of high quality digital audio"
2976    // Note there are bugs in that pseudocode, presumably due to them attempting
2977    // to rename the arrays nicely rather than representing the way their actual
2978    // implementation bounces buffers back and forth. As a result, even in the
2979    // "some formulars corrected" version, a direct implementation fails. These
2980    // are noted below as "paper bug".
2981
2982    // copy and reflect spectral data
2983    for (k=0; k < n2; ++k) u[k] = buffer[k];
2984    for (   ; k < n ; ++k) u[k] = -buffer[n - k - 1];
2985    // kernel from paper
2986    // step 1
2987    for (k=k2=k4=0; k < n4; k+=1, k2+=2, k4+=4) {
2988       v[n-k4-1] = (u[k4] - u[n-k4-1]) * A[k2]   - (u[k4+2] - u[n-k4-3])*A[k2+1];
2989       v[n-k4-3] = (u[k4] - u[n-k4-1]) * A[k2+1] + (u[k4+2] - u[n-k4-3])*A[k2];
2990    }
2991    // step 2
2992    for (k=k4=0; k < n8; k+=1, k4+=4) {
2993       w[n2+3+k4] = v[n2+3+k4] + v[k4+3];
2994       w[n2+1+k4] = v[n2+1+k4] + v[k4+1];
2995       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];
2996       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];
2997    }
2998    // step 3
2999    ld = ilog(n) - 1; // ilog is off-by-one from normal definitions
3000    for (l=0; l < ld-3; ++l) {
3001       int k0 = n >> (l+2), k1 = 1 << (l+3);
3002       int rlim = n >> (l+4), r4, r;
3003       int s2lim = 1 << (l+2), s2;
3004       for (r=r4=0; r < rlim; r4+=4,++r) {
3005          for (s2=0; s2 < s2lim; s2+=2) {
3006             u[n-1-k0*s2-r4] = w[n-1-k0*s2-r4] + w[n-1-k0*(s2+1)-r4];
3007             u[n-3-k0*s2-r4] = w[n-3-k0*s2-r4] + w[n-3-k0*(s2+1)-r4];
3008             u[n-1-k0*(s2+1)-r4] = (w[n-1-k0*s2-r4] - w[n-1-k0*(s2+1)-r4]) * A[r*k1]
3009                                 - (w[n-3-k0*s2-r4] - w[n-3-k0*(s2+1)-r4]) * A[r*k1+1];
3010             u[n-3-k0*(s2+1)-r4] = (w[n-3-k0*s2-r4] - w[n-3-k0*(s2+1)-r4]) * A[r*k1]
3011                                 + (w[n-1-k0*s2-r4] - w[n-1-k0*(s2+1)-r4]) * A[r*k1+1];
3012          }
3013       }
3014       if (l+1 < ld-3) {
3015          // paper bug: ping-ponging of u&w here is omitted
3016          memcpy(w, u, sizeof(u));
3017       }
3018    }
3019
3020    // step 4
3021    for (i=0; i < n8; ++i) {
3022       int j = bit_reverse(i) >> (32-ld+3);
3023       assert(j < n8);
3024       if (i == j) {
3025          // paper bug: original code probably swapped in place; if copying,
3026          //            need to directly copy in this case
3027          int i8 = i << 3;
3028          v[i8+1] = u[i8+1];
3029          v[i8+3] = u[i8+3];
3030          v[i8+5] = u[i8+5];
3031          v[i8+7] = u[i8+7];
3032       } else if (i < j) {
3033          int i8 = i << 3, j8 = j << 3;
3034          v[j8+1] = u[i8+1], v[i8+1] = u[j8 + 1];
3035          v[j8+3] = u[i8+3], v[i8+3] = u[j8 + 3];
3036          v[j8+5] = u[i8+5], v[i8+5] = u[j8 + 5];
3037          v[j8+7] = u[i8+7], v[i8+7] = u[j8 + 7];
3038       }
3039    }
3040    // step 5
3041    for (k=0; k < n2; ++k) {
3042       w[k] = v[k*2+1];
3043    }
3044    // step 6
3045    for (k=k2=k4=0; k < n8; ++k, k2 += 2, k4 += 4) {
3046       u[n-1-k2] = w[k4];
3047       u[n-2-k2] = w[k4+1];
3048       u[n3_4 - 1 - k2] = w[k4+2];
3049       u[n3_4 - 2 - k2] = w[k4+3];
3050    }
3051    // step 7
3052    for (k=k2=0; k < n8; ++k, k2 += 2) {
3053       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;
3054       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;
3055       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;
3056       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;
3057    }
3058    // step 8
3059    for (k=k2=0; k < n4; ++k,k2 += 2) {
3060       X[k]      = v[k2+n2]*B[k2  ] + v[k2+1+n2]*B[k2+1];
3061       X[n2-1-k] = v[k2+n2]*B[k2+1] - v[k2+1+n2]*B[k2  ];
3062    }
3063
3064    // decode kernel to output
3065    // determined the following value experimentally
3066    // (by first figuring out what made inverse_mdct_slow work); then matching that here
3067    // (probably vorbis encoder premultiplies by n or n/2, to save it on the decoder?)
3068    s = 0.5; // theoretically would be n4
3069
3070    // [[[ note! the s value of 0.5 is compensated for by the B[] in the current code,
3071    //     so it needs to use the "old" B values to behave correctly, or else
3072    //     set s to 1.0 ]]]
3073    for (i=0; i < n4  ; ++i) buffer[i] = s * X[i+n4];
3074    for (   ; i < n3_4; ++i) buffer[i] = -s * X[n3_4 - i - 1];
3075    for (   ; i < n   ; ++i) buffer[i] = -s * X[i - n3_4];
3076 }
3077 #endif
3078
3079 static float *get_window(vorb *f, int len)
3080 {
3081    len <<= 1;
3082    if (len == f->blocksize_0) return f->window[0];
3083    if (len == f->blocksize_1) return f->window[1];
3084    assert(0);
3085    return NULL;
3086 }
3087
3088 #ifndef STB_VORBIS_NO_DEFER_FLOOR
3089 typedef int16 YTYPE;
3090 #else
3091 typedef int YTYPE;
3092 #endif
3093 static int do_floor(vorb *f, Mapping *map, int i, int n, float *target, YTYPE *finalY, uint8 *step2_flag)
3094 {
3095    int n2 = n >> 1;
3096    int s = map->chan[i].mux, floor;
3097    floor = map->submap_floor[s];
3098    if (f->floor_types[floor] == 0) {
3099       return error(f, VORBIS_invalid_stream);
3100    } else {
3101       Floor1 *g = &f->floor_config[floor].floor1;
3102       int j,q;
3103       int lx = 0, ly = finalY[0] * g->floor1_multiplier;
3104       for (q=1; q < g->values; ++q) {
3105          j = g->sorted_order[q];
3106          #ifndef STB_VORBIS_NO_DEFER_FLOOR
3107          if (finalY[j] >= 0)
3108          #else
3109          if (step2_flag[j])
3110          #endif
3111          {
3112             int hy = finalY[j] * g->floor1_multiplier;
3113             int hx = g->Xlist[j];
3114             draw_line(target, lx,ly, hx,hy, n2);
3115             lx = hx, ly = hy;
3116          }
3117       }
3118       if (lx < n2)
3119          // optimization of: draw_line(target, lx,ly, n,ly, n2);
3120          for (j=lx; j < n2; ++j)
3121             LINE_OP(target[j], inverse_db_table[ly]);
3122    }
3123    return TRUE;
3124 }
3125
3126 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)
3127 {
3128    Mode *m;
3129    int i, n, prev, next, window_center;
3130    f->channel_buffer_start = f->channel_buffer_end = 0;
3131
3132   retry:
3133    if (f->eof) return FALSE;
3134    if (!maybe_start_packet(f))
3135       return FALSE;
3136    // check packet type
3137    if (get_bits(f,1) != 0) {
3138       if (IS_PUSH_MODE(f))
3139          return error(f,VORBIS_bad_packet_type);
3140       while (EOP != get8_packet(f));
3141       goto retry;
3142    }
3143
3144    if (f->alloc.alloc_buffer)
3145       assert(f->alloc.alloc_buffer_length_in_bytes == f->temp_offset);
3146
3147    i = get_bits(f, ilog(f->mode_count-1));
3148    if (i == EOP) return FALSE;
3149    if (i >= f->mode_count) return FALSE;
3150    *mode = i;
3151    m = f->mode_config + i;
3152    if (m->blockflag) {
3153       n = f->blocksize_1;
3154       prev = get_bits(f,1);
3155       next = get_bits(f,1);
3156    } else {
3157       prev = next = 0;
3158       n = f->blocksize_0;
3159    }
3160
3161 // WINDOWING
3162
3163    window_center = n >> 1;