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