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 */