]> 4ch.mooo.com Git - 16.git/blob - src/lib/dl/ext/flac/lpc.c
meh did some cleanings and i will work on mapread to mm thingy sometime soon! oops...
[16.git] / src / lib / dl / ext / flac / lpc.c
1 /* libFLAC - Free Lossless Audio Codec library
2  * Copyright (C) 2000,2001,2002,2003,2004,2005,2006,2007  Josh Coalson
3  *
4  * Redistribution and use in source and binary forms, with or without
5  * modification, are permitted provided that the following conditions
6  * are met:
7  *
8  * - Redistributions of source code must retain the above copyright
9  * notice, this list of conditions and the following disclaimer.
10  *
11  * - Redistributions in binary form must reproduce the above copyright
12  * notice, this list of conditions and the following disclaimer in the
13  * documentation and/or other materials provided with the distribution.
14  *
15  * - Neither the name of the Xiph.org Foundation nor the names of its
16  * contributors may be used to endorse or promote products derived from
17  * this software without specific prior written permission.
18  *
19  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
20  * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
21  * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
22  * A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
23  * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
24  * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
25  * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
26  * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
27  * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
28  * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
29  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
30  */
31
32 #if HAVE_CONFIG_H
33 #  include "config.h"
34 #endif
35
36 #include <math.h>
37 #include "flac/assert.h"
38 #include "flac/format.h"
39 #include "private/bitmath.h"
40 #include "private/lpc.h"
41 #if defined DEBUG || defined FLAC__OVERFLOW_DETECT || defined FLAC__OVERFLOW_DETECT_VERBOSE
42 #include <stdio.h>
43 #endif
44
45 #ifndef FLAC__INTEGER_ONLY_LIBRARY
46
47 #ifndef M_LN2
48 /* math.h in VC++ doesn't seem to have this (how Microsoft is that?) */
49 #define M_LN2 0.69314718055994530942
50 #endif
51
52 /* OPT: #undef'ing this may improve the speed on some architectures */
53 #define FLAC__LPC_UNROLLED_FILTER_LOOPS
54
55
56 void FLAC__lpc_window_data(const FLAC__int32 in[], const FLAC__real window[], FLAC__real out[], unsigned data_len)
57 {
58         unsigned i;
59         for(i = 0; i < data_len; i++)
60                 out[i] = in[i] * window[i];
61 }
62
63 void FLAC__lpc_compute_autocorrelation(const FLAC__real data[], unsigned data_len, unsigned lag, FLAC__real autoc[])
64 {
65         /* a readable, but slower, version */
66 #if 0
67         FLAC__real d;
68         unsigned i;
69
70         FLAC__ASSERT(lag > 0);
71         FLAC__ASSERT(lag <= data_len);
72
73         /*
74          * Technically we should subtract the mean first like so:
75          *   for(i = 0; i < data_len; i++)
76          *     data[i] -= mean;
77          * but it appears not to make enough of a difference to matter, and
78          * most signals are already closely centered around zero
79          */
80         while(lag--) {
81                 for(i = lag, d = 0.0; i < data_len; i++)
82                         d += data[i] * data[i - lag];
83                 autoc[lag] = d;
84         }
85 #endif
86
87         /*
88          * this version tends to run faster because of better data locality
89          * ('data_len' is usually much larger than 'lag')
90          */
91         FLAC__real d;
92         unsigned sample, coeff;
93         const unsigned limit = data_len - lag;
94
95         FLAC__ASSERT(lag > 0);
96         FLAC__ASSERT(lag <= data_len);
97
98         for(coeff = 0; coeff < lag; coeff++)
99                 autoc[coeff] = 0.0;
100         for(sample = 0; sample <= limit; sample++) {
101                 d = data[sample];
102                 for(coeff = 0; coeff < lag; coeff++)
103                         autoc[coeff] += d * data[sample+coeff];
104         }
105         for(; sample < data_len; sample++) {
106                 d = data[sample];
107                 for(coeff = 0; coeff < data_len - sample; coeff++)
108                         autoc[coeff] += d * data[sample+coeff];
109         }
110 }
111
112 void FLAC__lpc_compute_lp_coefficients(const FLAC__real autoc[], unsigned *max_order, FLAC__real lp_coeff[][FLAC__MAX_LPC_ORDER], FLAC__double error[])
113 {
114         unsigned i, j;
115         FLAC__double r, err, ref[FLAC__MAX_LPC_ORDER], lpc[FLAC__MAX_LPC_ORDER];
116
117         FLAC__ASSERT(0 != max_order);
118         FLAC__ASSERT(0 < *max_order);
119         FLAC__ASSERT(*max_order <= FLAC__MAX_LPC_ORDER);
120         FLAC__ASSERT(autoc[0] != 0.0);
121
122         err = autoc[0];
123
124         for(i = 0; i < *max_order; i++) {
125                 /* Sum up this iteration's reflection coefficient. */
126                 r = -autoc[i+1];
127                 for(j = 0; j < i; j++)
128                         r -= lpc[j] * autoc[i-j];
129                 ref[i] = (r/=err);
130
131                 /* Update LPC coefficients and total error. */
132                 lpc[i]=r;
133                 for(j = 0; j < (i>>1); j++) {
134                         FLAC__double tmp = lpc[j];
135                         lpc[j] += r * lpc[i-1-j];
136                         lpc[i-1-j] += r * tmp;
137                 }
138                 if(i & 1)
139                         lpc[j] += lpc[j] * r;
140
141                 err *= (1.0 - r * r);
142
143                 /* save this order */
144                 for(j = 0; j <= i; j++)
145                         lp_coeff[i][j] = (FLAC__real)(-lpc[j]); /* negate FIR filter coeff to get predictor coeff */
146                 error[i] = err;
147
148                 /* see SF bug #1601812 http://sourceforge.net/tracker/index.php?func=detail&aid=1601812&group_id=13478&atid=113478 */
149                 if(err == 0.0) {
150                         *max_order = i+1;
151                         return;
152                 }
153         }
154 }
155
156 int FLAC__lpc_quantize_coefficients(const FLAC__real lp_coeff[], unsigned order, unsigned precision, FLAC__int32 qlp_coeff[], int *shift)
157 {
158         unsigned i;
159         FLAC__double cmax;
160         FLAC__int32 qmax, qmin;
161
162         FLAC__ASSERT(precision > 0);
163         FLAC__ASSERT(precision >= FLAC__MIN_QLP_COEFF_PRECISION);
164
165         /* drop one bit for the sign; from here on out we consider only |lp_coeff[i]| */
166         precision--;
167         qmax = 1 << precision;
168         qmin = -qmax;
169         qmax--;
170
171         /* calc cmax = max( |lp_coeff[i]| ) */
172         cmax = 0.0;
173         for(i = 0; i < order; i++) {
174                 const FLAC__double d = fabs(lp_coeff[i]);
175                 if(d > cmax)
176                         cmax = d;
177         }
178
179         if(cmax <= 0.0) {
180                 /* => coefficients are all 0, which means our constant-detect didn't work */
181                 return 2;
182         }
183         else {
184                 const int max_shiftlimit = (1 << (FLAC__SUBFRAME_LPC_QLP_SHIFT_LEN-1)) - 1;
185                 const int min_shiftlimit = -max_shiftlimit - 1;
186                 int log2cmax;
187
188                 (void)frexp(cmax, &log2cmax);
189                 log2cmax--;
190                 *shift = (int)precision - log2cmax - 1;
191
192                 if(*shift > max_shiftlimit)
193                         *shift = max_shiftlimit;
194                 else if(*shift < min_shiftlimit)
195                         return 1;
196         }
197
198         if(*shift >= 0) {
199                 FLAC__double error = 0.0;
200                 FLAC__int32 q;
201                 for(i = 0; i < order; i++) {
202                         error += lp_coeff[i] * (1 << *shift);
203 #if 1 /* unfortunately lround() is C99 */
204                         if(error >= 0.0)
205                                 q = (FLAC__int32)(error + 0.5);
206                         else
207                                 q = (FLAC__int32)(error - 0.5);
208 #else
209                         q = lround(error);
210 #endif
211 #ifdef FLAC__OVERFLOW_DETECT
212                         if(q > qmax+1) /* we expect q==qmax+1 occasionally due to rounding */
213                                 fprintf(stderr,"FLAC__lpc_quantize_coefficients: quantizer overflow: q>qmax %d>%d shift=%d cmax=%f precision=%u lpc[%u]=%f\n",q,qmax,*shift,cmax,precision+1,i,lp_coeff[i]);
214                         else if(q < qmin)
215                                 fprintf(stderr,"FLAC__lpc_quantize_coefficients: quantizer overflow: q<qmin %d<%d shift=%d cmax=%f precision=%u lpc[%u]=%f\n",q,qmin,*shift,cmax,precision+1,i,lp_coeff[i]);
216 #endif
217                         if(q > qmax)
218                                 q = qmax;
219                         else if(q < qmin)
220                                 q = qmin;
221                         error -= q;
222                         qlp_coeff[i] = q;
223                 }
224         }
225         /* negative shift is very rare but due to design flaw, negative shift is
226          * a NOP in the decoder, so it must be handled specially by scaling down
227          * coeffs
228          */
229         else {
230                 const int nshift = -(*shift);
231                 FLAC__double error = 0.0;
232                 FLAC__int32 q;
233 #ifdef DEBUG
234                 fprintf(stderr,"FLAC__lpc_quantize_coefficients: negative shift=%d order=%u cmax=%f\n", *shift, order, cmax);
235 #endif
236                 for(i = 0; i < order; i++) {
237                         error += lp_coeff[i] / (1 << nshift);
238 #if 1 /* unfortunately lround() is C99 */
239                         if(error >= 0.0)
240                                 q = (FLAC__int32)(error + 0.5);
241                         else
242                                 q = (FLAC__int32)(error - 0.5);
243 #else
244                         q = lround(error);
245 #endif
246 #ifdef FLAC__OVERFLOW_DETECT
247                         if(q > qmax+1) /* we expect q==qmax+1 occasionally due to rounding */
248                                 fprintf(stderr,"FLAC__lpc_quantize_coefficients: quantizer overflow: q>qmax %d>%d shift=%d cmax=%f precision=%u lpc[%u]=%f\n",q,qmax,*shift,cmax,precision+1,i,lp_coeff[i]);
249                         else if(q < qmin)
250                                 fprintf(stderr,"FLAC__lpc_quantize_coefficients: quantizer overflow: q<qmin %d<%d shift=%d cmax=%f precision=%u lpc[%u]=%f\n",q,qmin,*shift,cmax,precision+1,i,lp_coeff[i]);
251 #endif
252                         if(q > qmax)
253                                 q = qmax;
254                         else if(q < qmin)
255                                 q = qmin;
256                         error -= q;
257                         qlp_coeff[i] = q;
258                 }
259                 *shift = 0;
260         }
261
262         return 0;
263 }
264
265 void FLAC__lpc_compute_residual_from_qlp_coefficients(const FLAC__int32 *data, unsigned data_len, const FLAC__int32 qlp_coeff[], unsigned order, int lp_quantization, FLAC__int32 residual[])
266 #if defined(FLAC__OVERFLOW_DETECT) || !defined(FLAC__LPC_UNROLLED_FILTER_LOOPS)
267 {
268         FLAC__int64 sumo;
269         unsigned i, j;
270         FLAC__int32 sum;
271         const FLAC__int32 *history;
272
273 #ifdef FLAC__OVERFLOW_DETECT_VERBOSE
274         fprintf(stderr,"FLAC__lpc_compute_residual_from_qlp_coefficients: data_len=%d, order=%u, lpq=%d",data_len,order,lp_quantization);
275         for(i=0;i<order;i++)
276                 fprintf(stderr,", q[%u]=%d",i,qlp_coeff[i]);
277         fprintf(stderr,"\n");
278 #endif
279         FLAC__ASSERT(order > 0);
280
281         for(i = 0; i < data_len; i++) {
282                 sumo = 0;
283                 sum = 0;
284                 history = data;
285                 for(j = 0; j < order; j++) {
286                         sum += qlp_coeff[j] * (*(--history));
287                         sumo += (FLAC__int64)qlp_coeff[j] * (FLAC__int64)(*history);
288 #if defined _MSC_VER
289                         if(sumo > 2147483647I64 || sumo < -2147483648I64)
290                                 fprintf(stderr,"FLAC__lpc_compute_residual_from_qlp_coefficients: OVERFLOW, i=%u, j=%u, c=%d, d=%d, sumo=%I64d\n",i,j,qlp_coeff[j],*history,sumo);
291 #else
292                         if(sumo > 2147483647ll || sumo < -2147483648ll)
293                                 fprintf(stderr,"FLAC__lpc_compute_residual_from_qlp_coefficients: OVERFLOW, i=%u, j=%u, c=%d, d=%d, sumo=%lld\n",i,j,qlp_coeff[j],*history,(long long)sumo);
294 #endif
295                 }
296                 *(residual++) = *(data++) - (sum >> lp_quantization);
297         }
298
299         /* Here's a slower but clearer version:
300         for(i = 0; i < data_len; i++) {
301                 sum = 0;
302                 for(j = 0; j < order; j++)
303                         sum += qlp_coeff[j] * data[i-j-1];
304                 residual[i] = data[i] - (sum >> lp_quantization);
305         }
306         */
307 }
308 #else /* fully unrolled version for normal use */
309 {
310         int i;
311         FLAC__int32 sum;
312
313         FLAC__ASSERT(order > 0);
314         FLAC__ASSERT(order <= 32);
315
316         /*
317          * We do unique versions up to 12th order since that's the subset limit.
318          * Also they are roughly ordered to match frequency of occurrence to
319          * minimize branching.
320          */
321         if(order <= 12) {
322                 if(order > 8) {
323                         if(order > 10) {
324                                 if(order == 12) {
325                                         for(i = 0; i < (int)data_len; i++) {
326                                                 sum = 0;
327                                                 sum += qlp_coeff[11] * data[i-12];
328                                                 sum += qlp_coeff[10] * data[i-11];
329                                                 sum += qlp_coeff[9] * data[i-10];
330                                                 sum += qlp_coeff[8] * data[i-9];
331                                                 sum += qlp_coeff[7] * data[i-8];
332                                                 sum += qlp_coeff[6] * data[i-7];
333                                                 sum += qlp_coeff[5] * data[i-6];
334                                                 sum += qlp_coeff[4] * data[i-5];
335                                                 sum += qlp_coeff[3] * data[i-4];
336                                                 sum += qlp_coeff[2] * data[i-3];
337                                                 sum += qlp_coeff[1] * data[i-2];
338                                                 sum += qlp_coeff[0] * data[i-1];
339                                                 residual[i] = data[i] - (sum >> lp_quantization);
340                                         }
341                                 }
342                                 else { /* order == 11 */
343                                         for(i = 0; i < (int)data_len; i++) {
344                                                 sum = 0;
345                                                 sum += qlp_coeff[10] * data[i-11];
346                                                 sum += qlp_coeff[9] * data[i-10];
347                                                 sum += qlp_coeff[8] * data[i-9];
348                                                 sum += qlp_coeff[7] * data[i-8];
349                                                 sum += qlp_coeff[6] * data[i-7];
350                                                 sum += qlp_coeff[5] * data[i-6];
351                                                 sum += qlp_coeff[4] * data[i-5];
352                                                 sum += qlp_coeff[3] * data[i-4];
353                                                 sum += qlp_coeff[2] * data[i-3];
354                                                 sum += qlp_coeff[1] * data[i-2];
355                                                 sum += qlp_coeff[0] * data[i-1];
356                                                 residual[i] = data[i] - (sum >> lp_quantization);
357                                         }
358                                 }
359                         }
360                         else {
361                                 if(order == 10) {
362                                         for(i = 0; i < (int)data_len; i++) {
363                                                 sum = 0;
364                                                 sum += qlp_coeff[9] * data[i-10];
365                                                 sum += qlp_coeff[8] * data[i-9];
366                                                 sum += qlp_coeff[7] * data[i-8];
367                                                 sum += qlp_coeff[6] * data[i-7];
368                                                 sum += qlp_coeff[5] * data[i-6];
369                                                 sum += qlp_coeff[4] * data[i-5];
370                                                 sum += qlp_coeff[3] * data[i-4];
371                                                 sum += qlp_coeff[2] * data[i-3];
372                                                 sum += qlp_coeff[1] * data[i-2];
373                                                 sum += qlp_coeff[0] * data[i-1];
374                                                 residual[i] = data[i] - (sum >> lp_quantization);
375                                         }
376                                 }
377                                 else { /* order == 9 */
378                                         for(i = 0; i < (int)data_len; i++) {
379                                                 sum = 0;
380                                                 sum += qlp_coeff[8] * data[i-9];
381                                                 sum += qlp_coeff[7] * data[i-8];
382                                                 sum += qlp_coeff[6] * data[i-7];
383                                                 sum += qlp_coeff[5] * data[i-6];
384                                                 sum += qlp_coeff[4] * data[i-5];
385                                                 sum += qlp_coeff[3] * data[i-4];
386                                                 sum += qlp_coeff[2] * data[i-3];
387                                                 sum += qlp_coeff[1] * data[i-2];
388                                                 sum += qlp_coeff[0] * data[i-1];
389                                                 residual[i] = data[i] - (sum >> lp_quantization);
390                                         }
391                                 }
392                         }
393                 }
394                 else if(order > 4) {
395                         if(order > 6) {
396                                 if(order == 8) {
397                                         for(i = 0; i < (int)data_len; i++) {
398                                                 sum = 0;
399                                                 sum += qlp_coeff[7] * data[i-8];
400                                                 sum += qlp_coeff[6] * data[i-7];
401                                                 sum += qlp_coeff[5] * data[i-6];
402                                                 sum += qlp_coeff[4] * data[i-5];
403                                                 sum += qlp_coeff[3] * data[i-4];
404                                                 sum += qlp_coeff[2] * data[i-3];
405                                                 sum += qlp_coeff[1] * data[i-2];
406                                                 sum += qlp_coeff[0] * data[i-1];
407                                                 residual[i] = data[i] - (sum >> lp_quantization);
408                                         }
409                                 }
410                                 else { /* order == 7 */
411                                         for(i = 0; i < (int)data_len; i++) {
412                                                 sum = 0;
413                                                 sum += qlp_coeff[6] * data[i-7];
414                                                 sum += qlp_coeff[5] * data[i-6];
415                                                 sum += qlp_coeff[4] * data[i-5];
416                                                 sum += qlp_coeff[3] * data[i-4];
417                                                 sum += qlp_coeff[2] * data[i-3];
418                                                 sum += qlp_coeff[1] * data[i-2];
419                                                 sum += qlp_coeff[0] * data[i-1];
420                                                 residual[i] = data[i] - (sum >> lp_quantization);
421                                         }
422                                 }
423                         }
424                         else {
425                                 if(order == 6) {
426                                         for(i = 0; i < (int)data_len; i++) {
427                                                 sum = 0;
428                                                 sum += qlp_coeff[5] * data[i-6];
429                                                 sum += qlp_coeff[4] * data[i-5];
430                                                 sum += qlp_coeff[3] * data[i-4];
431                                                 sum += qlp_coeff[2] * data[i-3];
432                                                 sum += qlp_coeff[1] * data[i-2];
433                                                 sum += qlp_coeff[0] * data[i-1];
434                                                 residual[i] = data[i] - (sum >> lp_quantization);
435                                         }
436                                 }
437                                 else { /* order == 5 */
438                                         for(i = 0; i < (int)data_len; i++) {
439                                                 sum = 0;
440                                                 sum += qlp_coeff[4] * data[i-5];
441                                                 sum += qlp_coeff[3] * data[i-4];
442                                                 sum += qlp_coeff[2] * data[i-3];
443                                                 sum += qlp_coeff[1] * data[i-2];
444                                                 sum += qlp_coeff[0] * data[i-1];
445                                                 residual[i] = data[i] - (sum >> lp_quantization);
446                                         }
447                                 }
448                         }
449                 }
450                 else {
451                         if(order > 2) {
452                                 if(order == 4) {
453                                         for(i = 0; i < (int)data_len; i++) {
454                                                 sum = 0;
455                                                 sum += qlp_coeff[3] * data[i-4];
456                                                 sum += qlp_coeff[2] * data[i-3];
457                                                 sum += qlp_coeff[1] * data[i-2];
458                                                 sum += qlp_coeff[0] * data[i-1];
459                                                 residual[i] = data[i] - (sum >> lp_quantization);
460                                         }
461                                 }
462                                 else { /* order == 3 */
463                                         for(i = 0; i < (int)data_len; i++) {
464                                                 sum = 0;
465                                                 sum += qlp_coeff[2] * data[i-3];
466                                                 sum += qlp_coeff[1] * data[i-2];
467                                                 sum += qlp_coeff[0] * data[i-1];
468                                                 residual[i] = data[i] - (sum >> lp_quantization);
469                                         }
470                                 }
471                         }
472                         else {
473                                 if(order == 2) {
474                                         for(i = 0; i < (int)data_len; i++) {
475                                                 sum = 0;
476                                                 sum += qlp_coeff[1] * data[i-2];
477                                                 sum += qlp_coeff[0] * data[i-1];
478                                                 residual[i] = data[i] - (sum >> lp_quantization);
479                                         }
480                                 }
481                                 else { /* order == 1 */
482                                         for(i = 0; i < (int)data_len; i++)
483                                                 residual[i] = data[i] - ((qlp_coeff[0] * data[i-1]) >> lp_quantization);
484                                 }
485                         }
486                 }
487         }
488         else { /* order > 12 */
489                 for(i = 0; i < (int)data_len; i++) {
490                         sum = 0;
491                         switch(order) {
492                                 case 32: sum += qlp_coeff[31] * data[i-32];
493                                 case 31: sum += qlp_coeff[30] * data[i-31];
494                                 case 30: sum += qlp_coeff[29] * data[i-30];
495                                 case 29: sum += qlp_coeff[28] * data[i-29];
496                                 case 28: sum += qlp_coeff[27] * data[i-28];
497                                 case 27: sum += qlp_coeff[26] * data[i-27];
498                                 case 26: sum += qlp_coeff[25] * data[i-26];
499                                 case 25: sum += qlp_coeff[24] * data[i-25];
500                                 case 24: sum += qlp_coeff[23] * data[i-24];
501                                 case 23: sum += qlp_coeff[22] * data[i-23];
502                                 case 22: sum += qlp_coeff[21] * data[i-22];
503                                 case 21: sum += qlp_coeff[20] * data[i-21];
504                                 case 20: sum += qlp_coeff[19] * data[i-20];
505                                 case 19: sum += qlp_coeff[18] * data[i-19];
506                                 case 18: sum += qlp_coeff[17] * data[i-18];
507                                 case 17: sum += qlp_coeff[16] * data[i-17];
508                                 case 16: sum += qlp_coeff[15] * data[i-16];
509                                 case 15: sum += qlp_coeff[14] * data[i-15];
510                                 case 14: sum += qlp_coeff[13] * data[i-14];
511                                 case 13: sum += qlp_coeff[12] * data[i-13];
512                                          sum += qlp_coeff[11] * data[i-12];
513                                          sum += qlp_coeff[10] * data[i-11];
514                                          sum += qlp_coeff[ 9] * data[i-10];
515                                          sum += qlp_coeff[ 8] * data[i- 9];
516                                          sum += qlp_coeff[ 7] * data[i- 8];
517                                          sum += qlp_coeff[ 6] * data[i- 7];
518                                          sum += qlp_coeff[ 5] * data[i- 6];
519                                          sum += qlp_coeff[ 4] * data[i- 5];
520                                          sum += qlp_coeff[ 3] * data[i- 4];
521                                          sum += qlp_coeff[ 2] * data[i- 3];
522                                          sum += qlp_coeff[ 1] * data[i- 2];
523                                          sum += qlp_coeff[ 0] * data[i- 1];
524                         }
525                         residual[i] = data[i] - (sum >> lp_quantization);
526                 }
527         }
528 }
529 #endif
530
531 void FLAC__lpc_compute_residual_from_qlp_coefficients_wide(const FLAC__int32 *data, unsigned data_len, const FLAC__int32 qlp_coeff[], unsigned order, int lp_quantization, FLAC__int32 residual[])
532 #if defined(FLAC__OVERFLOW_DETECT) || !defined(FLAC__LPC_UNROLLED_FILTER_LOOPS)
533 {
534         unsigned i, j;
535         FLAC__int64 sum;
536         const FLAC__int32 *history;
537
538 #ifdef FLAC__OVERFLOW_DETECT_VERBOSE
539         fprintf(stderr,"FLAC__lpc_compute_residual_from_qlp_coefficients_wide: data_len=%d, order=%u, lpq=%d",data_len,order,lp_quantization);
540         for(i=0;i<order;i++)
541                 fprintf(stderr,", q[%u]=%d",i,qlp_coeff[i]);
542         fprintf(stderr,"\n");
543 #endif
544         FLAC__ASSERT(order > 0);
545
546         for(i = 0; i < data_len; i++) {
547                 sum = 0;
548                 history = data;
549                 for(j = 0; j < order; j++)
550                         sum += (FLAC__int64)qlp_coeff[j] * (FLAC__int64)(*(--history));
551                 if(FLAC__bitmath_silog2_wide(sum >> lp_quantization) > 32) {
552 #if defined _MSC_VER
553                         fprintf(stderr,"FLAC__lpc_compute_residual_from_qlp_coefficients_wide: OVERFLOW, i=%u, sum=%I64d\n", i, sum >> lp_quantization);
554 #else
555                         fprintf(stderr,"FLAC__lpc_compute_residual_from_qlp_coefficients_wide: OVERFLOW, i=%u, sum=%lld\n", i, (long long)(sum >> lp_quantization));
556 #endif
557                         break;
558                 }
559                 if(FLAC__bitmath_silog2_wide((FLAC__int64)(*data) - (sum >> lp_quantization)) > 32) {
560 #if defined _MSC_VER
561                         fprintf(stderr,"FLAC__lpc_compute_residual_from_qlp_coefficients_wide: OVERFLOW, i=%u, data=%d, sum=%I64d, residual=%I64d\n", i, *data, sum >> lp_quantization, (FLAC__int64)(*data) - (sum >> lp_quantization));
562 #else
563                         fprintf(stderr,"FLAC__lpc_compute_residual_from_qlp_coefficients_wide: OVERFLOW, i=%u, data=%d, sum=%lld, residual=%lld\n", i, *data, (long long)(sum >> lp_quantization), (long long)((FLAC__int64)(*data) - (sum >> lp_quantization)));
564 #endif
565                         break;
566                 }
567                 *(residual++) = *(data++) - (FLAC__int32)(sum >> lp_quantization);
568         }
569 }
570 #else /* fully unrolled version for normal use */
571 {
572         int i;
573         FLAC__int64 sum;
574
575         FLAC__ASSERT(order > 0);
576         FLAC__ASSERT(order <= 32);
577
578         /*
579          * We do unique versions up to 12th order since that's the subset limit.
580          * Also they are roughly ordered to match frequency of occurrence to
581          * minimize branching.
582          */
583         if(order <= 12) {
584                 if(order > 8) {
585                         if(order > 10) {
586                                 if(order == 12) {
587                                         for(i = 0; i < (int)data_len; i++) {
588                                                 sum = 0;
589                                                 sum += qlp_coeff[11] * (FLAC__int64)data[i-12];
590                                                 sum += qlp_coeff[10] * (FLAC__int64)data[i-11];
591                                                 sum += qlp_coeff[9] * (FLAC__int64)data[i-10];
592                                                 sum += qlp_coeff[8] * (FLAC__int64)data[i-9];
593                                                 sum += qlp_coeff[7] * (FLAC__int64)data[i-8];
594                                                 sum += qlp_coeff[6] * (FLAC__int64)data[i-7];
595                                                 sum += qlp_coeff[5] * (FLAC__int64)data[i-6];
596                                                 sum += qlp_coeff[4] * (FLAC__int64)data[i-5];
597                                                 sum += qlp_coeff[3] * (FLAC__int64)data[i-4];
598                                                 sum += qlp_coeff[2] * (FLAC__int64)data[i-3];
599                                                 sum += qlp_coeff[1] * (FLAC__int64)data[i-2];
600                                                 sum += qlp_coeff[0] * (FLAC__int64)data[i-1];
601                                                 residual[i] = data[i] - (FLAC__int32)(sum >> lp_quantization);
602                                         }
603                                 }
604                                 else { /* order == 11 */
605                                         for(i = 0; i < (int)data_len; i++) {
606                                                 sum = 0;
607                                                 sum += qlp_coeff[10] * (FLAC__int64)data[i-11];
608                                                 sum += qlp_coeff[9] * (FLAC__int64)data[i-10];
609                                                 sum += qlp_coeff[8] * (FLAC__int64)data[i-9];
610                                                 sum += qlp_coeff[7] * (FLAC__int64)data[i-8];
611                                                 sum += qlp_coeff[6] * (FLAC__int64)data[i-7];
612                                                 sum += qlp_coeff[5] * (FLAC__int64)data[i-6];
613                                                 sum += qlp_coeff[4] * (FLAC__int64)data[i-5];
614                                                 sum += qlp_coeff[3] * (FLAC__int64)data[i-4];
615                                                 sum += qlp_coeff[2] * (FLAC__int64)data[i-3];
616                                                 sum += qlp_coeff[1] * (FLAC__int64)data[i-2];
617                                                 sum += qlp_coeff[0] * (FLAC__int64)data[i-1];
618                                                 residual[i] = data[i] - (FLAC__int32)(sum >> lp_quantization);
619                                         }
620                                 }
621                         }
622                         else {
623                                 if(order == 10) {
624                                         for(i = 0; i < (int)data_len; i++) {
625                                                 sum = 0;
626                                                 sum += qlp_coeff[9] * (FLAC__int64)data[i-10];
627                                                 sum += qlp_coeff[8] * (FLAC__int64)data[i-9];
628                                                 sum += qlp_coeff[7] * (FLAC__int64)data[i-8];
629                                                 sum += qlp_coeff[6] * (FLAC__int64)data[i-7];
630                                                 sum += qlp_coeff[5] * (FLAC__int64)data[i-6];
631                                                 sum += qlp_coeff[4] * (FLAC__int64)data[i-5];
632                                                 sum += qlp_coeff[3] * (FLAC__int64)data[i-4];
633                                                 sum += qlp_coeff[2] * (FLAC__int64)data[i-3];
634                                                 sum += qlp_coeff[1] * (FLAC__int64)data[i-2];
635                                                 sum += qlp_coeff[0] * (FLAC__int64)data[i-1];
636                                                 residual[i] = data[i] - (FLAC__int32)(sum >> lp_quantization);
637                                         }
638                                 }
639                                 else { /* order == 9 */
640                                         for(i = 0; i < (int)data_len; i++) {
641                                                 sum = 0;
642                                                 sum += qlp_coeff[8] * (FLAC__int64)data[i-9];
643                                                 sum += qlp_coeff[7] * (FLAC__int64)data[i-8];
644                                                 sum += qlp_coeff[6] * (FLAC__int64)data[i-7];
645                                                 sum += qlp_coeff[5] * (FLAC__int64)data[i-6];
646                                                 sum += qlp_coeff[4] * (FLAC__int64)data[i-5];
647                                                 sum += qlp_coeff[3] * (FLAC__int64)data[i-4];
648                                                 sum += qlp_coeff[2] * (FLAC__int64)data[i-3];
649                                                 sum += qlp_coeff[1] * (FLAC__int64)data[i-2];
650                                                 sum += qlp_coeff[0] * (FLAC__int64)data[i-1];
651                                                 residual[i] = data[i] - (FLAC__int32)(sum >> lp_quantization);
652                                         }
653                                 }
654                         }
655                 }
656                 else if(order > 4) {
657                         if(order > 6) {
658                                 if(order == 8) {
659                                         for(i = 0; i < (int)data_len; i++) {
660                                                 sum = 0;
661                                                 sum += qlp_coeff[7] * (FLAC__int64)data[i-8];
662                                                 sum += qlp_coeff[6] * (FLAC__int64)data[i-7];
663                                                 sum += qlp_coeff[5] * (FLAC__int64)data[i-6];
664                                                 sum += qlp_coeff[4] * (FLAC__int64)data[i-5];
665                                                 sum += qlp_coeff[3] * (FLAC__int64)data[i-4];
666                                                 sum += qlp_coeff[2] * (FLAC__int64)data[i-3];
667                                                 sum += qlp_coeff[1] * (FLAC__int64)data[i-2];
668                                                 sum += qlp_coeff[0] * (FLAC__int64)data[i-1];
669                                                 residual[i] = data[i] - (FLAC__int32)(sum >> lp_quantization);
670                                         }
671                                 }
672                                 else { /* order == 7 */
673                                         for(i = 0; i < (int)data_len; i++) {
674                                                 sum = 0;
675                                                 sum += qlp_coeff[6] * (FLAC__int64)data[i-7];
676                                                 sum += qlp_coeff[5] * (FLAC__int64)data[i-6];
677                                                 sum += qlp_coeff[4] * (FLAC__int64)data[i-5];
678                                                 sum += qlp_coeff[3] * (FLAC__int64)data[i-4];
679                                                 sum += qlp_coeff[2] * (FLAC__int64)data[i-3];
680                                                 sum += qlp_coeff[1] * (FLAC__int64)data[i-2];
681                                                 sum += qlp_coeff[0] * (FLAC__int64)data[i-1];
682                                                 residual[i] = data[i] - (FLAC__int32)(sum >> lp_quantization);
683                                         }
684                                 }
685                         }
686                         else {
687                                 if(order == 6) {
688                                         for(i = 0; i < (int)data_len; i++) {
689                                                 sum = 0;
690                                                 sum += qlp_coeff[5] * (FLAC__int64)data[i-6];
691                                                 sum += qlp_coeff[4] * (FLAC__int64)data[i-5];
692                                                 sum += qlp_coeff[3] * (FLAC__int64)data[i-4];
693                                                 sum += qlp_coeff[2] * (FLAC__int64)data[i-3];
694                                                 sum += qlp_coeff[1] * (FLAC__int64)data[i-2];
695                                                 sum += qlp_coeff[0] * (FLAC__int64)data[i-1];
696                                                 residual[i] = data[i] - (FLAC__int32)(sum >> lp_quantization);
697                                         }
698                                 }
699                                 else { /* order == 5 */
700                                         for(i = 0; i < (int)data_len; i++) {
701                                                 sum = 0;
702                                                 sum += qlp_coeff[4] * (FLAC__int64)data[i-5];
703                                                 sum += qlp_coeff[3] * (FLAC__int64)data[i-4];
704                                                 sum += qlp_coeff[2] * (FLAC__int64)data[i-3];
705                                                 sum += qlp_coeff[1] * (FLAC__int64)data[i-2];
706                                                 sum += qlp_coeff[0] * (FLAC__int64)data[i-1];
707                                                 residual[i] = data[i] - (FLAC__int32)(sum >> lp_quantization);
708                                         }
709                                 }
710                         }
711                 }
712                 else {
713                         if(order > 2) {
714                                 if(order == 4) {
715                                         for(i = 0; i < (int)data_len; i++) {
716                                                 sum = 0;
717                                                 sum += qlp_coeff[3] * (FLAC__int64)data[i-4];
718                                                 sum += qlp_coeff[2] * (FLAC__int64)data[i-3];
719                                                 sum += qlp_coeff[1] * (FLAC__int64)data[i-2];
720                                                 sum += qlp_coeff[0] * (FLAC__int64)data[i-1];
721                                                 residual[i] = data[i] - (FLAC__int32)(sum >> lp_quantization);
722                                         }
723                                 }
724                                 else { /* order == 3 */
725                                         for(i = 0; i < (int)data_len; i++) {
726                                                 sum = 0;
727                                                 sum += qlp_coeff[2] * (FLAC__int64)data[i-3];
728                                                 sum += qlp_coeff[1] * (FLAC__int64)data[i-2];
729                                                 sum += qlp_coeff[0] * (FLAC__int64)data[i-1];
730                                                 residual[i] = data[i] - (FLAC__int32)(sum >> lp_quantization);
731                                         }
732                                 }
733                         }
734                         else {
735                                 if(order == 2) {
736                                         for(i = 0; i < (int)data_len; i++) {
737                                                 sum = 0;
738                                                 sum += qlp_coeff[1] * (FLAC__int64)data[i-2];
739                                                 sum += qlp_coeff[0] * (FLAC__int64)data[i-1];
740                                                 residual[i] = data[i] - (FLAC__int32)(sum >> lp_quantization);
741                                         }
742                                 }
743                                 else { /* order == 1 */
744                                         for(i = 0; i < (int)data_len; i++)
745                                                 residual[i] = data[i] - (FLAC__int32)((qlp_coeff[0] * (FLAC__int64)data[i-1]) >> lp_quantization);
746                                 }
747                         }
748                 }
749         }
750         else { /* order > 12 */
751                 for(i = 0; i < (int)data_len; i++) {
752                         sum = 0;
753                         switch(order) {
754                                 case 32: sum += qlp_coeff[31] * (FLAC__int64)data[i-32];
755                                 case 31: sum += qlp_coeff[30] * (FLAC__int64)data[i-31];
756                                 case 30: sum += qlp_coeff[29] * (FLAC__int64)data[i-30];
757                                 case 29: sum += qlp_coeff[28] * (FLAC__int64)data[i-29];
758                                 case 28: sum += qlp_coeff[27] * (FLAC__int64)data[i-28];
759                                 case 27: sum += qlp_coeff[26] * (FLAC__int64)data[i-27];
760                                 case 26: sum += qlp_coeff[25] * (FLAC__int64)data[i-26];
761                                 case 25: sum += qlp_coeff[24] * (FLAC__int64)data[i-25];
762                                 case 24: sum += qlp_coeff[23] * (FLAC__int64)data[i-24];
763                                 case 23: sum += qlp_coeff[22] * (FLAC__int64)data[i-23];
764                                 case 22: sum += qlp_coeff[21] * (FLAC__int64)data[i-22];
765                                 case 21: sum += qlp_coeff[20] * (FLAC__int64)data[i-21];
766                                 case 20: sum += qlp_coeff[19] * (FLAC__int64)data[i-20];
767                                 case 19: sum += qlp_coeff[18] * (FLAC__int64)data[i-19];
768                                 case 18: sum += qlp_coeff[17] * (FLAC__int64)data[i-18];
769                                 case 17: sum += qlp_coeff[16] * (FLAC__int64)data[i-17];
770                                 case 16: sum += qlp_coeff[15] * (FLAC__int64)data[i-16];
771                                 case 15: sum += qlp_coeff[14] * (FLAC__int64)data[i-15];
772                                 case 14: sum += qlp_coeff[13] * (FLAC__int64)data[i-14];
773                                 case 13: sum += qlp_coeff[12] * (FLAC__int64)data[i-13];
774                                          sum += qlp_coeff[11] * (FLAC__int64)data[i-12];
775                                          sum += qlp_coeff[10] * (FLAC__int64)data[i-11];
776                                          sum += qlp_coeff[ 9] * (FLAC__int64)data[i-10];
777                                          sum += qlp_coeff[ 8] * (FLAC__int64)data[i- 9];
778                                          sum += qlp_coeff[ 7] * (FLAC__int64)data[i- 8];
779                                          sum += qlp_coeff[ 6] * (FLAC__int64)data[i- 7];
780                                          sum += qlp_coeff[ 5] * (FLAC__int64)data[i- 6];
781                                          sum += qlp_coeff[ 4] * (FLAC__int64)data[i- 5];
782                                          sum += qlp_coeff[ 3] * (FLAC__int64)data[i- 4];
783                                          sum += qlp_coeff[ 2] * (FLAC__int64)data[i- 3];
784                                          sum += qlp_coeff[ 1] * (FLAC__int64)data[i- 2];
785                                          sum += qlp_coeff[ 0] * (FLAC__int64)data[i- 1];
786                         }
787                         residual[i] = data[i] - (FLAC__int32)(sum >> lp_quantization);
788                 }
789         }
790 }
791 #endif
792
793 #endif /* !defined FLAC__INTEGER_ONLY_LIBRARY */
794
795 void FLAC__lpc_restore_signal(const FLAC__int32 residual[], unsigned data_len, const FLAC__int32 qlp_coeff[], unsigned order, int lp_quantization, FLAC__int32 data[])
796 #if defined(FLAC__OVERFLOW_DETECT) || !defined(FLAC__LPC_UNROLLED_FILTER_LOOPS)
797 {
798         FLAC__int64 sumo;
799         unsigned i, j;
800         FLAC__int32 sum;
801         const FLAC__int32 *r = residual, *history;
802
803 #ifdef FLAC__OVERFLOW_DETECT_VERBOSE
804         fprintf(stderr,"FLAC__lpc_restore_signal: data_len=%d, order=%u, lpq=%d",data_len,order,lp_quantization);
805         for(i=0;i<order;i++)
806                 fprintf(stderr,", q[%u]=%d",i,qlp_coeff[i]);
807         fprintf(stderr,"\n");
808 #endif
809         FLAC__ASSERT(order > 0);
810
811         for(i = 0; i < data_len; i++) {
812                 sumo = 0;
813                 sum = 0;
814                 history = data;
815                 for(j = 0; j < order; j++) {
816                         sum += qlp_coeff[j] * (*(--history));
817                         sumo += (FLAC__int64)qlp_coeff[j] * (FLAC__int64)(*history);
818 #if defined _MSC_VER
819                         if(sumo > 2147483647I64 || sumo < -2147483648I64)
820                                 fprintf(stderr,"FLAC__lpc_restore_signal: OVERFLOW, i=%u, j=%u, c=%d, d=%d, sumo=%I64d\n",i,j,qlp_coeff[j],*history,sumo);
821 #else
822                         if(sumo > 2147483647ll || sumo < -2147483648ll)
823                                 fprintf(stderr,"FLAC__lpc_restore_signal: OVERFLOW, i=%u, j=%u, c=%d, d=%d, sumo=%lld\n",i,j,qlp_coeff[j],*history,(long long)sumo);
824 #endif
825                 }
826                 *(data++) = *(r++) + (sum >> lp_quantization);
827         }
828
829         /* Here's a slower but clearer version:
830         for(i = 0; i < data_len; i++) {
831                 sum = 0;
832                 for(j = 0; j < order; j++)
833                         sum += qlp_coeff[j] * data[i-j-1];
834                 data[i] = residual[i] + (sum >> lp_quantization);
835         }
836         */
837 }
838 #else /* fully unrolled version for normal use */
839 {
840         int i;
841         FLAC__int32 sum;
842
843         FLAC__ASSERT(order > 0);
844         FLAC__ASSERT(order <= 32);
845
846         /*
847          * We do unique versions up to 12th order since that's the subset limit.
848          * Also they are roughly ordered to match frequency of occurrence to
849          * minimize branching.
850          */
851         if(order <= 12) {
852                 if(order > 8) {
853                         if(order > 10) {
854                                 if(order == 12) {
855                                         for(i = 0; i < (int)data_len; i++) {
856                                                 sum = 0;
857                                                 sum += qlp_coeff[11] * data[i-12];
858                                                 sum += qlp_coeff[10] * data[i-11];
859                                                 sum += qlp_coeff[9] * data[i-10];
860                                                 sum += qlp_coeff[8] * data[i-9];
861                                                 sum += qlp_coeff[7] * data[i-8];
862                                                 sum += qlp_coeff[6] * data[i-7];
863                                                 sum += qlp_coeff[5] * data[i-6];
864                                                 sum += qlp_coeff[4] * data[i-5];
865                                                 sum += qlp_coeff[3] * data[i-4];
866                                                 sum += qlp_coeff[2] * data[i-3];
867                                                 sum += qlp_coeff[1] * data[i-2];
868                                                 sum += qlp_coeff[0] * data[i-1];
869                                                 data[i] = residual[i] + (sum >> lp_quantization);
870                                         }
871                                 }
872                                 else { /* order == 11 */
873                                         for(i = 0; i < (int)data_len; i++) {
874                                                 sum = 0;
875                                                 sum += qlp_coeff[10] * data[i-11];
876                                                 sum += qlp_coeff[9] * data[i-10];
877                                                 sum += qlp_coeff[8] * data[i-9];
878                                                 sum += qlp_coeff[7] * data[i-8];
879                                                 sum += qlp_coeff[6] * data[i-7];
880                                                 sum += qlp_coeff[5] * data[i-6];
881                                                 sum += qlp_coeff[4] * data[i-5];
882                                                 sum += qlp_coeff[3] * data[i-4];
883                                                 sum += qlp_coeff[2] * data[i-3];
884                                                 sum += qlp_coeff[1] * data[i-2];
885                                                 sum += qlp_coeff[0] * data[i-1];
886                                                 data[i] = residual[i] + (sum >> lp_quantization);
887                                         }
888                                 }
889                         }
890                         else {
891                                 if(order == 10) {
892                                         for(i = 0; i < (int)data_len; i++) {
893                                                 sum = 0;
894                                                 sum += qlp_coeff[9] * data[i-10];
895                                                 sum += qlp_coeff[8] * data[i-9];
896                                                 sum += qlp_coeff[7] * data[i-8];
897                                                 sum += qlp_coeff[6] * data[i-7];
898                                                 sum += qlp_coeff[5] * data[i-6];
899                                                 sum += qlp_coeff[4] * data[i-5];
900                                                 sum += qlp_coeff[3] * data[i-4];
901                                                 sum += qlp_coeff[2] * data[i-3];
902                                                 sum += qlp_coeff[1] * data[i-2];
903                                                 sum += qlp_coeff[0] * data[i-1];
904                                                 data[i] = residual[i] + (sum >> lp_quantization);
905                                         }
906                                 }
907                                 else { /* order == 9 */
908                                         for(i = 0; i < (int)data_len; i++) {
909                                                 sum = 0;
910                                                 sum += qlp_coeff[8] * data[i-9];
911                                                 sum += qlp_coeff[7] * data[i-8];
912                                                 sum += qlp_coeff[6] * data[i-7];
913                                                 sum += qlp_coeff[5] * data[i-6];
914                                                 sum += qlp_coeff[4] * data[i-5];
915                                                 sum += qlp_coeff[3] * data[i-4];
916                                                 sum += qlp_coeff[2] * data[i-3];
917                                                 sum += qlp_coeff[1] * data[i-2];
918                                                 sum += qlp_coeff[0] * data[i-1];
919                                                 data[i] = residual[i] + (sum >> lp_quantization);
920                                         }
921                                 }
922                         }
923                 }
924                 else if(order > 4) {
925                         if(order > 6) {
926                                 if(order == 8) {
927                                         for(i = 0; i < (int)data_len; i++) {
928                                                 sum = 0;
929                                                 sum += qlp_coeff[7] * data[i-8];
930                                                 sum += qlp_coeff[6] * data[i-7];
931                                                 sum += qlp_coeff[5] * data[i-6];
932                                                 sum += qlp_coeff[4] * data[i-5];
933                                                 sum += qlp_coeff[3] * data[i-4];
934                                                 sum += qlp_coeff[2] * data[i-3];
935                                                 sum += qlp_coeff[1] * data[i-2];
936                                                 sum += qlp_coeff[0] * data[i-1];
937                                                 data[i] = residual[i] + (sum >> lp_quantization);
938                                         }
939                                 }
940                                 else { /* order == 7 */
941                                         for(i = 0; i < (int)data_len; i++) {
942                                                 sum = 0;
943                                                 sum += qlp_coeff[6] * data[i-7];
944                                                 sum += qlp_coeff[5] * data[i-6];
945                                                 sum += qlp_coeff[4] * data[i-5];
946                                                 sum += qlp_coeff[3] * data[i-4];
947                                                 sum += qlp_coeff[2] * data[i-3];
948                                                 sum += qlp_coeff[1] * data[i-2];
949                                                 sum += qlp_coeff[0] * data[i-1];
950                                                 data[i] = residual[i] + (sum >> lp_quantization);
951                                         }
952                                 }
953                         }
954                         else {
955                                 if(order == 6) {
956                                         for(i = 0; i < (int)data_len; i++) {
957                                                 sum = 0;
958                                                 sum += qlp_coeff[5] * data[i-6];
959                                                 sum += qlp_coeff[4] * data[i-5];
960                                                 sum += qlp_coeff[3] * data[i-4];
961                                                 sum += qlp_coeff[2] * data[i-3];
962                                                 sum += qlp_coeff[1] * data[i-2];
963                                                 sum += qlp_coeff[0] * data[i-1];
964                                                 data[i] = residual[i] + (sum >> lp_quantization);
965                                         }
966                                 }
967                                 else { /* order == 5 */
968                                         for(i = 0; i < (int)data_len; i++) {
969                                                 sum = 0;
970                                                 sum += qlp_coeff[4] * data[i-5];
971                                                 sum += qlp_coeff[3] * data[i-4];
972                                                 sum += qlp_coeff[2] * data[i-3];
973                                                 sum += qlp_coeff[1] * data[i-2];
974                                                 sum += qlp_coeff[0] * data[i-1];
975                                                 data[i] = residual[i] + (sum >> lp_quantization);
976                                         }
977                                 }
978                         }
979                 }
980                 else {
981                         if(order > 2) {
982                                 if(order == 4) {
983                                         for(i = 0; i < (int)data_len; i++) {
984                                                 sum = 0;
985                                                 sum += qlp_coeff[3] * data[i-4];
986                                                 sum += qlp_coeff[2] * data[i-3];
987                                                 sum += qlp_coeff[1] * data[i-2];
988                                                 sum += qlp_coeff[0] * data[i-1];
989                                                 data[i] = residual[i] + (sum >> lp_quantization);
990                                         }
991                                 }
992                                 else { /* order == 3 */
993                                         for(i = 0; i < (int)data_len; i++) {
994                                                 sum = 0;
995                                                 sum += qlp_coeff[2] * data[i-3];
996                                                 sum += qlp_coeff[1] * data[i-2];
997                                                 sum += qlp_coeff[0] * data[i-1];
998                                                 data[i] = residual[i] + (sum >> lp_quantization);
999                                         }
1000                                 }
1001                         }
1002                         else {
1003                                 if(order == 2) {
1004                                         for(i = 0; i < (int)data_len; i++) {
1005                                                 sum = 0;
1006                                                 sum += qlp_coeff[1] * data[i-2];
1007                                                 sum += qlp_coeff[0] * data[i-1];
1008                                                 data[i] = residual[i] + (sum >> lp_quantization);
1009                                         }
1010                                 }
1011                                 else { /* order == 1 */
1012                                         for(i = 0; i < (int)data_len; i++)
1013                                                 data[i] = residual[i] + ((qlp_coeff[0] * data[i-1]) >> lp_quantization);
1014                                 }
1015                         }
1016                 }
1017         }
1018         else { /* order > 12 */
1019                 for(i = 0; i < (int)data_len; i++) {
1020                         sum = 0;
1021                         switch(order) {
1022                                 case 32: sum += qlp_coeff[31] * data[i-32];
1023                                 case 31: sum += qlp_coeff[30] * data[i-31];
1024                                 case 30: sum += qlp_coeff[29] * data[i-30];
1025                                 case 29: sum += qlp_coeff[28] * data[i-29];
1026                                 case 28: sum += qlp_coeff[27] * data[i-28];
1027                                 case 27: sum += qlp_coeff[26] * data[i-27];
1028                                 case 26: sum += qlp_coeff[25] * data[i-26];
1029                                 case 25: sum += qlp_coeff[24] * data[i-25];
1030                                 case 24: sum += qlp_coeff[23] * data[i-24];
1031                                 case 23: sum += qlp_coeff[22] * data[i-23];
1032                                 case 22: sum += qlp_coeff[21] * data[i-22];
1033                                 case 21: sum += qlp_coeff[20] * data[i-21];
1034                                 case 20: sum += qlp_coeff[19] * data[i-20];
1035                                 case 19: sum += qlp_coeff[18] * data[i-19];
1036                                 case 18: sum += qlp_coeff[17] * data[i-18];
1037                                 case 17: sum += qlp_coeff[16] * data[i-17];
1038                                 case 16: sum += qlp_coeff[15] * data[i-16];
1039                                 case 15: sum += qlp_coeff[14] * data[i-15];
1040                                 case 14: sum += qlp_coeff[13] * data[i-14];
1041                                 case 13: sum += qlp_coeff[12] * data[i-13];
1042                                          sum += qlp_coeff[11] * data[i-12];
1043                                          sum += qlp_coeff[10] * data[i-11];
1044                                          sum += qlp_coeff[ 9] * data[i-10];
1045                                          sum += qlp_coeff[ 8] * data[i- 9];
1046                                          sum += qlp_coeff[ 7] * data[i- 8];
1047                                          sum += qlp_coeff[ 6] * data[i- 7];
1048                                          sum += qlp_coeff[ 5] * data[i- 6];
1049                                          sum += qlp_coeff[ 4] * data[i- 5];
1050                                          sum += qlp_coeff[ 3] * data[i- 4];
1051                                          sum += qlp_coeff[ 2] * data[i- 3];
1052                                          sum += qlp_coeff[ 1] * data[i- 2];
1053                                          sum += qlp_coeff[ 0] * data[i- 1];
1054                         }
1055                         data[i] = residual[i] + (sum >> lp_quantization);
1056                 }
1057         }
1058 }
1059 #endif
1060
1061 void FLAC__lpc_restore_signal_wide(const FLAC__int32 residual[], unsigned data_len, const FLAC__int32 qlp_coeff[], unsigned order, int lp_quantization, FLAC__int32 data[])
1062 #if defined(FLAC__OVERFLOW_DETECT) || !defined(FLAC__LPC_UNROLLED_FILTER_LOOPS)
1063 {
1064         unsigned i, j;
1065         FLAC__int64 sum;
1066         const FLAC__int32 *r = residual, *history;
1067
1068 #ifdef FLAC__OVERFLOW_DETECT_VERBOSE
1069         fprintf(stderr,"FLAC__lpc_restore_signal_wide: data_len=%d, order=%u, lpq=%d",data_len,order,lp_quantization);
1070         for(i=0;i<order;i++)
1071                 fprintf(stderr,", q[%u]=%d",i,qlp_coeff[i]);
1072         fprintf(stderr,"\n");
1073 #endif
1074         FLAC__ASSERT(order > 0);
1075
1076         for(i = 0; i < data_len; i++) {
1077                 sum = 0;
1078                 history = data;
1079                 for(j = 0; j < order; j++)
1080                         sum += (FLAC__int64)qlp_coeff[j] * (FLAC__int64)(*(--history));
1081                 if(FLAC__bitmath_silog2_wide(sum >> lp_quantization) > 32) {
1082 #ifdef _MSC_VER
1083                         fprintf(stderr,"FLAC__lpc_restore_signal_wide: OVERFLOW, i=%u, sum=%I64d\n", i, sum >> lp_quantization);
1084 #else
1085                         fprintf(stderr,"FLAC__lpc_restore_signal_wide: OVERFLOW, i=%u, sum=%lld\n", i, (long long)(sum >> lp_quantization));
1086 #endif
1087                         break;
1088                 }
1089                 if(FLAC__bitmath_silog2_wide((FLAC__int64)(*r) + (sum >> lp_quantization)) > 32) {
1090 #ifdef _MSC_VER
1091                         fprintf(stderr,"FLAC__lpc_restore_signal_wide: OVERFLOW, i=%u, residual=%d, sum=%I64d, data=%I64d\n", i, *r, sum >> lp_quantization, (FLAC__int64)(*r) + (sum >> lp_quantization));
1092 #else
1093                         fprintf(stderr,"FLAC__lpc_restore_signal_wide: OVERFLOW, i=%u, residual=%d, sum=%lld, data=%lld\n", i, *r, (long long)(sum >> lp_quantization), (long long)((FLAC__int64)(*r) + (sum >> lp_quantization)));
1094 #endif
1095                         break;
1096                 }
1097                 *(data++) = *(r++) + (FLAC__int32)(sum >> lp_quantization);
1098         }
1099 }
1100 #else /* fully unrolled version for normal use */
1101 {
1102         int i;
1103         FLAC__int64 sum;
1104
1105         FLAC__ASSERT(order > 0);
1106         FLAC__ASSERT(order <= 32);
1107
1108         /*
1109          * We do unique versions up to 12th order since that's the subset limit.
1110          * Also they are roughly ordered to match frequency of occurrence to
1111          * minimize branching.
1112          */
1113         if(order <= 12) {
1114                 if(order > 8) {
1115                         if(order > 10) {
1116                                 if(order == 12) {
1117                                         for(i = 0; i < (int)data_len; i++) {
1118                                                 sum = 0;
1119                                                 sum += qlp_coeff[11] * (FLAC__int64)data[i-12];
1120                                                 sum += qlp_coeff[10] * (FLAC__int64)data[i-11];
1121                                                 sum += qlp_coeff[9] * (FLAC__int64)data[i-10];
1122                                                 sum += qlp_coeff[8] * (FLAC__int64)data[i-9];
1123                                                 sum += qlp_coeff[7] * (FLAC__int64)data[i-8];
1124                                                 sum += qlp_coeff[6] * (FLAC__int64)data[i-7];
1125                                                 sum += qlp_coeff[5] * (FLAC__int64)data[i-6];
1126                                                 sum += qlp_coeff[4] * (FLAC__int64)data[i-5];
1127                                                 sum += qlp_coeff[3] * (FLAC__int64)data[i-4];
1128                                                 sum += qlp_coeff[2] * (FLAC__int64)data[i-3];
1129                                                 sum += qlp_coeff[1] * (FLAC__int64)data[i-2];
1130                                                 sum += qlp_coeff[0] * (FLAC__int64)data[i-1];
1131                                                 data[i] = residual[i] + (FLAC__int32)(sum >> lp_quantization);
1132                                         }
1133                                 }
1134                                 else { /* order == 11 */
1135                                         for(i = 0; i < (int)data_len; i++) {
1136                                                 sum = 0;
1137                                                 sum += qlp_coeff[10] * (FLAC__int64)data[i-11];
1138                                                 sum += qlp_coeff[9] * (FLAC__int64)data[i-10];
1139                                                 sum += qlp_coeff[8] * (FLAC__int64)data[i-9];
1140                                                 sum += qlp_coeff[7] * (FLAC__int64)data[i-8];
1141                                                 sum += qlp_coeff[6] * (FLAC__int64)data[i-7];
1142                                                 sum += qlp_coeff[5] * (FLAC__int64)data[i-6];
1143                                                 sum += qlp_coeff[4] * (FLAC__int64)data[i-5];
1144                                                 sum += qlp_coeff[3] * (FLAC__int64)data[i-4];
1145                                                 sum += qlp_coeff[2] * (FLAC__int64)data[i-3];
1146                                                 sum += qlp_coeff[1] * (FLAC__int64)data[i-2];
1147                                                 sum += qlp_coeff[0] * (FLAC__int64)data[i-1];
1148                                                 data[i] = residual[i] + (FLAC__int32)(sum >> lp_quantization);
1149                                         }
1150                                 }
1151                         }
1152                         else {
1153                                 if(order == 10) {
1154                                         for(i = 0; i < (int)data_len; i++) {
1155                                                 sum = 0;
1156                                                 sum += qlp_coeff[9] * (FLAC__int64)data[i-10];
1157                                                 sum += qlp_coeff[8] * (FLAC__int64)data[i-9];
1158                                                 sum += qlp_coeff[7] * (FLAC__int64)data[i-8];
1159                                                 sum += qlp_coeff[6] * (FLAC__int64)data[i-7];
1160                                                 sum += qlp_coeff[5] * (FLAC__int64)data[i-6];
1161                                                 sum += qlp_coeff[4] * (FLAC__int64)data[i-5];
1162                                                 sum += qlp_coeff[3] * (FLAC__int64)data[i-4];
1163                                                 sum += qlp_coeff[2] * (FLAC__int64)data[i-3];
1164                                                 sum += qlp_coeff[1] * (FLAC__int64)data[i-2];
1165                                                 sum += qlp_coeff[0] * (FLAC__int64)data[i-1];
1166                                                 data[i] = residual[i] + (FLAC__int32)(sum >> lp_quantization);
1167                                         }
1168                                 }
1169                                 else { /* order == 9 */
1170                                         for(i = 0; i < (int)data_len; i++) {
1171                                                 sum = 0;
1172                                                 sum += qlp_coeff[8] * (FLAC__int64)data[i-9];
1173                                                 sum += qlp_coeff[7] * (FLAC__int64)data[i-8];
1174                                                 sum += qlp_coeff[6] * (FLAC__int64)data[i-7];
1175                                                 sum += qlp_coeff[5] * (FLAC__int64)data[i-6];
1176                                                 sum += qlp_coeff[4] * (FLAC__int64)data[i-5];
1177                                                 sum += qlp_coeff[3] * (FLAC__int64)data[i-4];
1178                                                 sum += qlp_coeff[2] * (FLAC__int64)data[i-3];
1179                                                 sum += qlp_coeff[1] * (FLAC__int64)data[i-2];
1180                                                 sum += qlp_coeff[0] * (FLAC__int64)data[i-1];
1181                                                 data[i] = residual[i] + (FLAC__int32)(sum >> lp_quantization);
1182                                         }
1183                                 }
1184                         }
1185                 }
1186                 else if(order > 4) {
1187                         if(order > 6) {
1188                                 if(order == 8) {
1189                                         for(i = 0; i < (int)data_len; i++) {
1190                                                 sum = 0;
1191                                                 sum += qlp_coeff[7] * (FLAC__int64)data[i-8];
1192                                                 sum += qlp_coeff[6] * (FLAC__int64)data[i-7];
1193                                                 sum += qlp_coeff[5] * (FLAC__int64)data[i-6];
1194                                                 sum += qlp_coeff[4] * (FLAC__int64)data[i-5];
1195                                                 sum += qlp_coeff[3] * (FLAC__int64)data[i-4];
1196                                                 sum += qlp_coeff[2] * (FLAC__int64)data[i-3];
1197                                                 sum += qlp_coeff[1] * (FLAC__int64)data[i-2];
1198                                                 sum += qlp_coeff[0] * (FLAC__int64)data[i-1];
1199                                                 data[i] = residual[i] + (FLAC__int32)(sum >> lp_quantization);
1200                                         }
1201                                 }
1202                                 else { /* order == 7 */
1203                                         for(i = 0; i < (int)data_len; i++) {
1204                                                 sum = 0;
1205                                                 sum += qlp_coeff[6] * (FLAC__int64)data[i-7];
1206                                                 sum += qlp_coeff[5] * (FLAC__int64)data[i-6];
1207                                                 sum += qlp_coeff[4] * (FLAC__int64)data[i-5];
1208                                                 sum += qlp_coeff[3] * (FLAC__int64)data[i-4];
1209                                                 sum += qlp_coeff[2] * (FLAC__int64)data[i-3];
1210                                                 sum += qlp_coeff[1] * (FLAC__int64)data[i-2];
1211                                                 sum += qlp_coeff[0] * (FLAC__int64)data[i-1];
1212                                                 data[i] = residual[i] + (FLAC__int32)(sum >> lp_quantization);
1213                                         }
1214                                 }
1215                         }
1216                         else {
1217                                 if(order == 6) {
1218                                         for(i = 0; i < (int)data_len; i++) {
1219                                                 sum = 0;
1220                                                 sum += qlp_coeff[5] * (FLAC__int64)data[i-6];
1221                                                 sum += qlp_coeff[4] * (FLAC__int64)data[i-5];
1222                                                 sum += qlp_coeff[3] * (FLAC__int64)data[i-4];
1223                                                 sum += qlp_coeff[2] * (FLAC__int64)data[i-3];
1224                                                 sum += qlp_coeff[1] * (FLAC__int64)data[i-2];
1225                                                 sum += qlp_coeff[0] * (FLAC__int64)data[i-1];
1226                                                 data[i] = residual[i] + (FLAC__int32)(sum >> lp_quantization);
1227                                         }
1228                                 }
1229                                 else { /* order == 5 */
1230                                         for(i = 0; i < (int)data_len; i++) {
1231                                                 sum = 0;
1232                                                 sum += qlp_coeff[4] * (FLAC__int64)data[i-5];
1233                                                 sum += qlp_coeff[3] * (FLAC__int64)data[i-4];
1234                                                 sum += qlp_coeff[2] * (FLAC__int64)data[i-3];
1235                                                 sum += qlp_coeff[1] * (FLAC__int64)data[i-2];
1236                                                 sum += qlp_coeff[0] * (FLAC__int64)data[i-1];
1237                                                 data[i] = residual[i] + (FLAC__int32)(sum >> lp_quantization);
1238                                         }
1239                                 }
1240                         }
1241                 }
1242                 else {
1243                         if(order > 2) {
1244                                 if(order == 4) {
1245                                         for(i = 0; i < (int)data_len; i++) {
1246                                                 sum = 0;
1247                                                 sum += qlp_coeff[3] * (FLAC__int64)data[i-4];
1248                                                 sum += qlp_coeff[2] * (FLAC__int64)data[i-3];
1249                                                 sum += qlp_coeff[1] * (FLAC__int64)data[i-2];
1250                                                 sum += qlp_coeff[0] * (FLAC__int64)data[i-1];
1251                                                 data[i] = residual[i] + (FLAC__int32)(sum >> lp_quantization);
1252                                         }
1253                                 }
1254                                 else { /* order == 3 */
1255                                         for(i = 0; i < (int)data_len; i++) {
1256                                                 sum = 0;
1257                                                 sum += qlp_coeff[2] * (FLAC__int64)data[i-3];
1258                                                 sum += qlp_coeff[1] * (FLAC__int64)data[i-2];
1259                                                 sum += qlp_coeff[0] * (FLAC__int64)data[i-1];
1260                                                 data[i] = residual[i] + (FLAC__int32)(sum >> lp_quantization);
1261                                         }
1262                                 }
1263                         }
1264                         else {
1265                                 if(order == 2) {
1266                                         for(i = 0; i < (int)data_len; i++) {
1267                                                 sum = 0;
1268                                                 sum += qlp_coeff[1] * (FLAC__int64)data[i-2];
1269                                                 sum += qlp_coeff[0] * (FLAC__int64)data[i-1];
1270                                                 data[i] = residual[i] + (FLAC__int32)(sum >> lp_quantization);
1271                                         }
1272                                 }
1273                                 else { /* order == 1 */
1274                                         for(i = 0; i < (int)data_len; i++)
1275                                                 data[i] = residual[i] + (FLAC__int32)((qlp_coeff[0] * (FLAC__int64)data[i-1]) >> lp_quantization);
1276                                 }
1277                         }
1278                 }
1279         }
1280         else { /* order > 12 */
1281                 for(i = 0; i < (int)data_len; i++) {
1282                         sum = 0;
1283                         switch(order) {
1284                                 case 32: sum += qlp_coeff[31] * (FLAC__int64)data[i-32];
1285                                 case 31: sum += qlp_coeff[30] * (FLAC__int64)data[i-31];
1286                                 case 30: sum += qlp_coeff[29] * (FLAC__int64)data[i-30];
1287                                 case 29: sum += qlp_coeff[28] * (FLAC__int64)data[i-29];
1288                                 case 28: sum += qlp_coeff[27] * (FLAC__int64)data[i-28];
1289                                 case 27: sum += qlp_coeff[26] * (FLAC__int64)data[i-27];
1290                                 case 26: sum += qlp_coeff[25] * (FLAC__int64)data[i-26];
1291                                 case 25: sum += qlp_coeff[24] * (FLAC__int64)data[i-25];
1292                                 case 24: sum += qlp_coeff[23] * (FLAC__int64)data[i-24];
1293                                 case 23: sum += qlp_coeff[22] * (FLAC__int64)data[i-23];
1294                                 case 22: sum += qlp_coeff[21] * (FLAC__int64)data[i-22];
1295                                 case 21: sum += qlp_coeff[20] * (FLAC__int64)data[i-21];
1296                                 case 20: sum += qlp_coeff[19] * (FLAC__int64)data[i-20];
1297                                 case 19: sum += qlp_coeff[18] * (FLAC__int64)data[i-19];
1298                                 case 18: sum += qlp_coeff[17] * (FLAC__int64)data[i-18];
1299                                 case 17: sum += qlp_coeff[16] * (FLAC__int64)data[i-17];
1300                                 case 16: sum += qlp_coeff[15] * (FLAC__int64)data[i-16];
1301                                 case 15: sum += qlp_coeff[14] * (FLAC__int64)data[i-15];
1302                                 case 14: sum += qlp_coeff[13] * (FLAC__int64)data[i-14];
1303                                 case 13: sum += qlp_coeff[12] * (FLAC__int64)data[i-13];
1304                                          sum += qlp_coeff[11] * (FLAC__int64)data[i-12];
1305                                          sum += qlp_coeff[10] * (FLAC__int64)data[i-11];
1306                                          sum += qlp_coeff[ 9] * (FLAC__int64)data[i-10];
1307                                          sum += qlp_coeff[ 8] * (FLAC__int64)data[i- 9];
1308                                          sum += qlp_coeff[ 7] * (FLAC__int64)data[i- 8];
1309                                          sum += qlp_coeff[ 6] * (FLAC__int64)data[i- 7];
1310                                          sum += qlp_coeff[ 5] * (FLAC__int64)data[i- 6];
1311                                          sum += qlp_coeff[ 4] * (FLAC__int64)data[i- 5];
1312                                          sum += qlp_coeff[ 3] * (FLAC__int64)data[i- 4];
1313                                          sum += qlp_coeff[ 2] * (FLAC__int64)data[i- 3];
1314                                          sum += qlp_coeff[ 1] * (FLAC__int64)data[i- 2];
1315                                          sum += qlp_coeff[ 0] * (FLAC__int64)data[i- 1];
1316                         }
1317                         data[i] = residual[i] + (FLAC__int32)(sum >> lp_quantization);
1318                 }
1319         }
1320 }
1321 #endif
1322
1323 #ifndef FLAC__INTEGER_ONLY_LIBRARY
1324
1325 FLAC__double FLAC__lpc_compute_expected_bits_per_residual_sample(FLAC__double lpc_error, unsigned total_samples)
1326 {
1327         FLAC__double error_scale;
1328
1329         FLAC__ASSERT(total_samples > 0);
1330
1331         error_scale = 0.5 * M_LN2 * M_LN2 / (FLAC__double)total_samples;
1332
1333         return FLAC__lpc_compute_expected_bits_per_residual_sample_with_error_scale(lpc_error, error_scale);
1334 }
1335
1336 FLAC__double FLAC__lpc_compute_expected_bits_per_residual_sample_with_error_scale(FLAC__double lpc_error, FLAC__double error_scale)
1337 {
1338         if(lpc_error > 0.0) {
1339                 FLAC__double bps = (FLAC__double)0.5 * log(error_scale * lpc_error) / M_LN2;
1340                 if(bps >= 0.0)
1341                         return bps;
1342                 else
1343                         return 0.0;
1344         }
1345         else if(lpc_error < 0.0) { /* error should not be negative but can happen due to inadequate floating-point resolution */
1346                 return 1e32;
1347         }
1348         else {
1349                 return 0.0;
1350         }
1351 }
1352
1353 unsigned FLAC__lpc_compute_best_order(const FLAC__double lpc_error[], unsigned max_order, unsigned total_samples, unsigned overhead_bits_per_order)
1354 {
1355         unsigned order, index, best_index; /* 'index' the index into lpc_error; index==order-1 since lpc_error[0] is for order==1, lpc_error[1] is for order==2, etc */
1356         FLAC__double bits, best_bits, error_scale;
1357
1358         FLAC__ASSERT(max_order > 0);
1359         FLAC__ASSERT(total_samples > 0);
1360
1361         error_scale = 0.5 * M_LN2 * M_LN2 / (FLAC__double)total_samples;
1362
1363         best_index = 0;
1364         best_bits = (unsigned)(-1);
1365
1366         for(index = 0, order = 1; index < max_order; index++, order++) {
1367                 bits = FLAC__lpc_compute_expected_bits_per_residual_sample_with_error_scale(lpc_error[index], error_scale) * (FLAC__double)(total_samples - order) + (FLAC__double)(order * overhead_bits_per_order);
1368                 if(bits < best_bits) {
1369                         best_index = index;
1370                         best_bits = bits;
1371                 }
1372         }
1373
1374         return best_index+1; /* +1 since index of lpc_error[] is order-1 */
1375 }
1376
1377 #endif /* !defined FLAC__INTEGER_ONLY_LIBRARY */