V2: Use a single definition of MIN and MAX in sources
[flac.git] / src / libFLAC / fixed.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 #include <string.h>
38 #include "private/bitmath.h"
39 #include "private/fixed.h"
40 #include "private/macros.h"
41 #include "FLAC/assert.h"
42
43 #ifndef M_LN2
44 /* math.h in VC++ doesn't seem to have this (how Microsoft is that?) */
45 #define M_LN2 0.69314718055994530942
46 #endif
47
48 #ifdef local_abs
49 #undef local_abs
50 #endif
51 #define local_abs(x) ((unsigned)((x)<0? -(x) : (x)))
52
53 #ifdef FLAC__INTEGER_ONLY_LIBRARY
54 /* rbps stands for residual bits per sample
55  *
56  *             (ln(2) * err)
57  * rbps = log  (-----------)
58  *           2 (     n     )
59  */
60 static FLAC__fixedpoint local__compute_rbps_integerized(FLAC__uint32 err, FLAC__uint32 n)
61 {
62         FLAC__uint32 rbps;
63         unsigned bits; /* the number of bits required to represent a number */
64         int fracbits; /* the number of bits of rbps that comprise the fractional part */
65
66         FLAC__ASSERT(sizeof(rbps) == sizeof(FLAC__fixedpoint));
67         FLAC__ASSERT(err > 0);
68         FLAC__ASSERT(n > 0);
69
70         FLAC__ASSERT(n <= FLAC__MAX_BLOCK_SIZE);
71         if(err <= n)
72                 return 0;
73         /*
74          * The above two things tell us 1) n fits in 16 bits; 2) err/n > 1.
75          * These allow us later to know we won't lose too much precision in the
76          * fixed-point division (err<<fracbits)/n.
77          */
78
79         fracbits = (8*sizeof(err)) - (FLAC__bitmath_ilog2(err)+1);
80
81         err <<= fracbits;
82         err /= n;
83         /* err now holds err/n with fracbits fractional bits */
84
85         /*
86          * Whittle err down to 16 bits max.  16 significant bits is enough for
87          * our purposes.
88          */
89         FLAC__ASSERT(err > 0);
90         bits = FLAC__bitmath_ilog2(err)+1;
91         if(bits > 16) {
92                 err >>= (bits-16);
93                 fracbits -= (bits-16);
94         }
95         rbps = (FLAC__uint32)err;
96
97         /* Multiply by fixed-point version of ln(2), with 16 fractional bits */
98         rbps *= FLAC__FP_LN2;
99         fracbits += 16;
100         FLAC__ASSERT(fracbits >= 0);
101
102         /* FLAC__fixedpoint_log2 requires fracbits%4 to be 0 */
103         {
104                 const int f = fracbits & 3;
105                 if(f) {
106                         rbps >>= f;
107                         fracbits -= f;
108                 }
109         }
110
111         rbps = FLAC__fixedpoint_log2(rbps, fracbits, (unsigned)(-1));
112
113         if(rbps == 0)
114                 return 0;
115
116         /*
117          * The return value must have 16 fractional bits.  Since the whole part
118          * of the base-2 log of a 32 bit number must fit in 5 bits, and fracbits
119          * must be >= -3, these assertion allows us to be able to shift rbps
120          * left if necessary to get 16 fracbits without losing any bits of the
121          * whole part of rbps.
122          *
123          * There is a slight chance due to accumulated error that the whole part
124          * will require 6 bits, so we use 6 in the assertion.  Really though as
125          * long as it fits in 13 bits (32 - (16 - (-3))) we are fine.
126          */
127         FLAC__ASSERT((int)FLAC__bitmath_ilog2(rbps)+1 <= fracbits + 6);
128         FLAC__ASSERT(fracbits >= -3);
129
130         /* now shift the decimal point into place */
131         if(fracbits < 16)
132                 return rbps << (16-fracbits);
133         else if(fracbits > 16)
134                 return rbps >> (fracbits-16);
135         else
136                 return rbps;
137 }
138
139 static FLAC__fixedpoint local__compute_rbps_wide_integerized(FLAC__uint64 err, FLAC__uint32 n)
140 {
141         FLAC__uint32 rbps;
142         unsigned bits; /* the number of bits required to represent a number */
143         int fracbits; /* the number of bits of rbps that comprise the fractional part */
144
145         FLAC__ASSERT(sizeof(rbps) == sizeof(FLAC__fixedpoint));
146         FLAC__ASSERT(err > 0);
147         FLAC__ASSERT(n > 0);
148
149         FLAC__ASSERT(n <= FLAC__MAX_BLOCK_SIZE);
150         if(err <= n)
151                 return 0;
152         /*
153          * The above two things tell us 1) n fits in 16 bits; 2) err/n > 1.
154          * These allow us later to know we won't lose too much precision in the
155          * fixed-point division (err<<fracbits)/n.
156          */
157
158         fracbits = (8*sizeof(err)) - (FLAC__bitmath_ilog2_wide(err)+1);
159
160         err <<= fracbits;
161         err /= n;
162         /* err now holds err/n with fracbits fractional bits */
163
164         /*
165          * Whittle err down to 16 bits max.  16 significant bits is enough for
166          * our purposes.
167          */
168         FLAC__ASSERT(err > 0);
169         bits = FLAC__bitmath_ilog2_wide(err)+1;
170         if(bits > 16) {
171                 err >>= (bits-16);
172                 fracbits -= (bits-16);
173         }
174         rbps = (FLAC__uint32)err;
175
176         /* Multiply by fixed-point version of ln(2), with 16 fractional bits */
177         rbps *= FLAC__FP_LN2;
178         fracbits += 16;
179         FLAC__ASSERT(fracbits >= 0);
180
181         /* FLAC__fixedpoint_log2 requires fracbits%4 to be 0 */
182         {
183                 const int f = fracbits & 3;
184                 if(f) {
185                         rbps >>= f;
186                         fracbits -= f;
187                 }
188         }
189
190         rbps = FLAC__fixedpoint_log2(rbps, fracbits, (unsigned)(-1));
191
192         if(rbps == 0)
193                 return 0;
194
195         /*
196          * The return value must have 16 fractional bits.  Since the whole part
197          * of the base-2 log of a 32 bit number must fit in 5 bits, and fracbits
198          * must be >= -3, these assertion allows us to be able to shift rbps
199          * left if necessary to get 16 fracbits without losing any bits of the
200          * whole part of rbps.
201          *
202          * There is a slight chance due to accumulated error that the whole part
203          * will require 6 bits, so we use 6 in the assertion.  Really though as
204          * long as it fits in 13 bits (32 - (16 - (-3))) we are fine.
205          */
206         FLAC__ASSERT((int)FLAC__bitmath_ilog2(rbps)+1 <= fracbits + 6);
207         FLAC__ASSERT(fracbits >= -3);
208
209         /* now shift the decimal point into place */
210         if(fracbits < 16)
211                 return rbps << (16-fracbits);
212         else if(fracbits > 16)
213                 return rbps >> (fracbits-16);
214         else
215                 return rbps;
216 }
217 #endif
218
219 #ifndef FLAC__INTEGER_ONLY_LIBRARY
220 unsigned FLAC__fixed_compute_best_predictor(const FLAC__int32 data[], unsigned data_len, FLAC__float residual_bits_per_sample[FLAC__MAX_FIXED_ORDER+1])
221 #else
222 unsigned FLAC__fixed_compute_best_predictor(const FLAC__int32 data[], unsigned data_len, FLAC__fixedpoint residual_bits_per_sample[FLAC__MAX_FIXED_ORDER+1])
223 #endif
224 {
225         FLAC__int32 last_error_0 = data[-1];
226         FLAC__int32 last_error_1 = data[-1] - data[-2];
227         FLAC__int32 last_error_2 = last_error_1 - (data[-2] - data[-3]);
228         FLAC__int32 last_error_3 = last_error_2 - (data[-2] - 2*data[-3] + data[-4]);
229         FLAC__int32 error, save;
230         FLAC__uint32 total_error_0 = 0, total_error_1 = 0, total_error_2 = 0, total_error_3 = 0, total_error_4 = 0;
231         unsigned i, order;
232
233         for(i = 0; i < data_len; i++) {
234                 error  = data[i]     ; total_error_0 += local_abs(error);                      save = error;
235                 error -= last_error_0; total_error_1 += local_abs(error); last_error_0 = save; save = error;
236                 error -= last_error_1; total_error_2 += local_abs(error); last_error_1 = save; save = error;
237                 error -= last_error_2; total_error_3 += local_abs(error); last_error_2 = save; save = error;
238                 error -= last_error_3; total_error_4 += local_abs(error); last_error_3 = save;
239         }
240
241         if(total_error_0 < flac_min(flac_min(flac_min(total_error_1, total_error_2), total_error_3), total_error_4))
242                 order = 0;
243         else if(total_error_1 < flac_min(flac_min(total_error_2, total_error_3), total_error_4))
244                 order = 1;
245         else if(total_error_2 < flac_min(total_error_3, total_error_4))
246                 order = 2;
247         else if(total_error_3 < total_error_4)
248                 order = 3;
249         else
250                 order = 4;
251
252         /* Estimate the expected number of bits per residual signal sample. */
253         /* 'total_error*' is linearly related to the variance of the residual */
254         /* signal, so we use it directly to compute E(|x|) */
255         FLAC__ASSERT(data_len > 0 || total_error_0 == 0);
256         FLAC__ASSERT(data_len > 0 || total_error_1 == 0);
257         FLAC__ASSERT(data_len > 0 || total_error_2 == 0);
258         FLAC__ASSERT(data_len > 0 || total_error_3 == 0);
259         FLAC__ASSERT(data_len > 0 || total_error_4 == 0);
260 #ifndef FLAC__INTEGER_ONLY_LIBRARY
261         residual_bits_per_sample[0] = (FLAC__float)((total_error_0 > 0) ? log(M_LN2 * (FLAC__double)total_error_0 / (FLAC__double)data_len) / M_LN2 : 0.0);
262         residual_bits_per_sample[1] = (FLAC__float)((total_error_1 > 0) ? log(M_LN2 * (FLAC__double)total_error_1 / (FLAC__double)data_len) / M_LN2 : 0.0);
263         residual_bits_per_sample[2] = (FLAC__float)((total_error_2 > 0) ? log(M_LN2 * (FLAC__double)total_error_2 / (FLAC__double)data_len) / M_LN2 : 0.0);
264         residual_bits_per_sample[3] = (FLAC__float)((total_error_3 > 0) ? log(M_LN2 * (FLAC__double)total_error_3 / (FLAC__double)data_len) / M_LN2 : 0.0);
265         residual_bits_per_sample[4] = (FLAC__float)((total_error_4 > 0) ? log(M_LN2 * (FLAC__double)total_error_4 / (FLAC__double)data_len) / M_LN2 : 0.0);
266 #else
267         residual_bits_per_sample[0] = (total_error_0 > 0) ? local__compute_rbps_integerized(total_error_0, data_len) : 0;
268         residual_bits_per_sample[1] = (total_error_1 > 0) ? local__compute_rbps_integerized(total_error_1, data_len) : 0;
269         residual_bits_per_sample[2] = (total_error_2 > 0) ? local__compute_rbps_integerized(total_error_2, data_len) : 0;
270         residual_bits_per_sample[3] = (total_error_3 > 0) ? local__compute_rbps_integerized(total_error_3, data_len) : 0;
271         residual_bits_per_sample[4] = (total_error_4 > 0) ? local__compute_rbps_integerized(total_error_4, data_len) : 0;
272 #endif
273
274         return order;
275 }
276
277 #ifndef FLAC__INTEGER_ONLY_LIBRARY
278 unsigned FLAC__fixed_compute_best_predictor_wide(const FLAC__int32 data[], unsigned data_len, FLAC__float residual_bits_per_sample[FLAC__MAX_FIXED_ORDER+1])
279 #else
280 unsigned FLAC__fixed_compute_best_predictor_wide(const FLAC__int32 data[], unsigned data_len, FLAC__fixedpoint residual_bits_per_sample[FLAC__MAX_FIXED_ORDER+1])
281 #endif
282 {
283         FLAC__int32 last_error_0 = data[-1];
284         FLAC__int32 last_error_1 = data[-1] - data[-2];
285         FLAC__int32 last_error_2 = last_error_1 - (data[-2] - data[-3]);
286         FLAC__int32 last_error_3 = last_error_2 - (data[-2] - 2*data[-3] + data[-4]);
287         FLAC__int32 error, save;
288         /* total_error_* are 64-bits to avoid overflow when encoding
289          * erratic signals when the bits-per-sample and blocksize are
290          * large.
291          */
292         FLAC__uint64 total_error_0 = 0, total_error_1 = 0, total_error_2 = 0, total_error_3 = 0, total_error_4 = 0;
293         unsigned i, order;
294
295         for(i = 0; i < data_len; i++) {
296                 error  = data[i]     ; total_error_0 += local_abs(error);                      save = error;
297                 error -= last_error_0; total_error_1 += local_abs(error); last_error_0 = save; save = error;
298                 error -= last_error_1; total_error_2 += local_abs(error); last_error_1 = save; save = error;
299                 error -= last_error_2; total_error_3 += local_abs(error); last_error_2 = save; save = error;
300                 error -= last_error_3; total_error_4 += local_abs(error); last_error_3 = save;
301         }
302
303         if(total_error_0 < flac_min(flac_min(flac_min(total_error_1, total_error_2), total_error_3), total_error_4))
304                 order = 0;
305         else if(total_error_1 < flac_min(flac_min(total_error_2, total_error_3), total_error_4))
306                 order = 1;
307         else if(total_error_2 < flac_min(total_error_3, total_error_4))
308                 order = 2;
309         else if(total_error_3 < total_error_4)
310                 order = 3;
311         else
312                 order = 4;
313
314         /* Estimate the expected number of bits per residual signal sample. */
315         /* 'total_error*' is linearly related to the variance of the residual */
316         /* signal, so we use it directly to compute E(|x|) */
317         FLAC__ASSERT(data_len > 0 || total_error_0 == 0);
318         FLAC__ASSERT(data_len > 0 || total_error_1 == 0);
319         FLAC__ASSERT(data_len > 0 || total_error_2 == 0);
320         FLAC__ASSERT(data_len > 0 || total_error_3 == 0);
321         FLAC__ASSERT(data_len > 0 || total_error_4 == 0);
322 #ifndef FLAC__INTEGER_ONLY_LIBRARY
323 #if defined _MSC_VER || defined __MINGW32__
324         /* with MSVC you have to spoon feed it the casting */
325         residual_bits_per_sample[0] = (FLAC__float)((total_error_0 > 0) ? log(M_LN2 * (FLAC__double)(FLAC__int64)total_error_0 / (FLAC__double)data_len) / M_LN2 : 0.0);
326         residual_bits_per_sample[1] = (FLAC__float)((total_error_1 > 0) ? log(M_LN2 * (FLAC__double)(FLAC__int64)total_error_1 / (FLAC__double)data_len) / M_LN2 : 0.0);
327         residual_bits_per_sample[2] = (FLAC__float)((total_error_2 > 0) ? log(M_LN2 * (FLAC__double)(FLAC__int64)total_error_2 / (FLAC__double)data_len) / M_LN2 : 0.0);
328         residual_bits_per_sample[3] = (FLAC__float)((total_error_3 > 0) ? log(M_LN2 * (FLAC__double)(FLAC__int64)total_error_3 / (FLAC__double)data_len) / M_LN2 : 0.0);
329         residual_bits_per_sample[4] = (FLAC__float)((total_error_4 > 0) ? log(M_LN2 * (FLAC__double)(FLAC__int64)total_error_4 / (FLAC__double)data_len) / M_LN2 : 0.0);
330 #else
331         residual_bits_per_sample[0] = (FLAC__float)((total_error_0 > 0) ? log(M_LN2 * (FLAC__double)total_error_0 / (FLAC__double)data_len) / M_LN2 : 0.0);
332         residual_bits_per_sample[1] = (FLAC__float)((total_error_1 > 0) ? log(M_LN2 * (FLAC__double)total_error_1 / (FLAC__double)data_len) / M_LN2 : 0.0);
333         residual_bits_per_sample[2] = (FLAC__float)((total_error_2 > 0) ? log(M_LN2 * (FLAC__double)total_error_2 / (FLAC__double)data_len) / M_LN2 : 0.0);
334         residual_bits_per_sample[3] = (FLAC__float)((total_error_3 > 0) ? log(M_LN2 * (FLAC__double)total_error_3 / (FLAC__double)data_len) / M_LN2 : 0.0);
335         residual_bits_per_sample[4] = (FLAC__float)((total_error_4 > 0) ? log(M_LN2 * (FLAC__double)total_error_4 / (FLAC__double)data_len) / M_LN2 : 0.0);
336 #endif
337 #else
338         residual_bits_per_sample[0] = (total_error_0 > 0) ? local__compute_rbps_wide_integerized(total_error_0, data_len) : 0;
339         residual_bits_per_sample[1] = (total_error_1 > 0) ? local__compute_rbps_wide_integerized(total_error_1, data_len) : 0;
340         residual_bits_per_sample[2] = (total_error_2 > 0) ? local__compute_rbps_wide_integerized(total_error_2, data_len) : 0;
341         residual_bits_per_sample[3] = (total_error_3 > 0) ? local__compute_rbps_wide_integerized(total_error_3, data_len) : 0;
342         residual_bits_per_sample[4] = (total_error_4 > 0) ? local__compute_rbps_wide_integerized(total_error_4, data_len) : 0;
343 #endif
344
345         return order;
346 }
347
348 void FLAC__fixed_compute_residual(const FLAC__int32 data[], unsigned data_len, unsigned order, FLAC__int32 residual[])
349 {
350         const int idata_len = (int)data_len;
351         int i;
352
353         switch(order) {
354                 case 0:
355                         FLAC__ASSERT(sizeof(residual[0]) == sizeof(data[0]));
356                         memcpy(residual, data, sizeof(residual[0])*data_len);
357                         break;
358                 case 1:
359                         for(i = 0; i < idata_len; i++)
360                                 residual[i] = data[i] - data[i-1];
361                         break;
362                 case 2:
363                         for(i = 0; i < idata_len; i++)
364 #if 1 /* OPT: may be faster with some compilers on some systems */
365                                 residual[i] = data[i] - (data[i-1] << 1) + data[i-2];
366 #else
367                                 residual[i] = data[i] - 2*data[i-1] + data[i-2];
368 #endif
369                         break;
370                 case 3:
371                         for(i = 0; i < idata_len; i++)
372 #if 1 /* OPT: may be faster with some compilers on some systems */
373                                 residual[i] = data[i] - (((data[i-1]-data[i-2])<<1) + (data[i-1]-data[i-2])) - data[i-3];
374 #else
375                                 residual[i] = data[i] - 3*data[i-1] + 3*data[i-2] - data[i-3];
376 #endif
377                         break;
378                 case 4:
379                         for(i = 0; i < idata_len; i++)
380 #if 1 /* OPT: may be faster with some compilers on some systems */
381                                 residual[i] = data[i] - ((data[i-1]+data[i-3])<<2) + ((data[i-2]<<2) + (data[i-2]<<1)) + data[i-4];
382 #else
383                                 residual[i] = data[i] - 4*data[i-1] + 6*data[i-2] - 4*data[i-3] + data[i-4];
384 #endif
385                         break;
386                 default:
387                         FLAC__ASSERT(0);
388         }
389 }
390
391 void FLAC__fixed_restore_signal(const FLAC__int32 residual[], unsigned data_len, unsigned order, FLAC__int32 data[])
392 {
393         int i, idata_len = (int)data_len;
394
395         switch(order) {
396                 case 0:
397                         FLAC__ASSERT(sizeof(residual[0]) == sizeof(data[0]));
398                         memcpy(data, residual, sizeof(residual[0])*data_len);
399                         break;
400                 case 1:
401                         for(i = 0; i < idata_len; i++)
402                                 data[i] = residual[i] + data[i-1];
403                         break;
404                 case 2:
405                         for(i = 0; i < idata_len; i++)
406 #if 1 /* OPT: may be faster with some compilers on some systems */
407                                 data[i] = residual[i] + (data[i-1]<<1) - data[i-2];
408 #else
409                                 data[i] = residual[i] + 2*data[i-1] - data[i-2];
410 #endif
411                         break;
412                 case 3:
413                         for(i = 0; i < idata_len; i++)
414 #if 1 /* OPT: may be faster with some compilers on some systems */
415                                 data[i] = residual[i] + (((data[i-1]-data[i-2])<<1) + (data[i-1]-data[i-2])) + data[i-3];
416 #else
417                                 data[i] = residual[i] + 3*data[i-1] - 3*data[i-2] + data[i-3];
418 #endif
419                         break;
420                 case 4:
421                         for(i = 0; i < idata_len; i++)
422 #if 1 /* OPT: may be faster with some compilers on some systems */
423                                 data[i] = residual[i] + ((data[i-1]+data[i-3])<<2) - ((data[i-2]<<2) + (data[i-2]<<1)) - data[i-4];
424 #else
425                                 data[i] = residual[i] + 4*data[i-1] - 6*data[i-2] + 4*data[i-3] - data[i-4];
426 #endif
427                         break;
428                 default:
429                         FLAC__ASSERT(0);
430         }
431 }