libFLAC/cpu.c: Get rid of OS_IS_ANDROID function
[flac.git] / src / libFLAC / fixed.c
1 /* libFLAC - Free Lossless Audio Codec library
2  * Copyright (C) 2000-2009  Josh Coalson
3  * Copyright (C) 2011-2016  Xiph.Org Foundation
4  *
5  * Redistribution and use in source and binary forms, with or without
6  * modification, are permitted provided that the following conditions
7  * are met:
8  *
9  * - Redistributions of source code must retain the above copyright
10  * notice, this list of conditions and the following disclaimer.
11  *
12  * - Redistributions in binary form must reproduce the above copyright
13  * notice, this list of conditions and the following disclaimer in the
14  * documentation and/or other materials provided with the distribution.
15  *
16  * - Neither the name of the Xiph.org Foundation nor the names of its
17  * contributors may be used to endorse or promote products derived from
18  * this software without specific prior written permission.
19  *
20  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21  * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22  * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
23  * A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
24  * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
25  * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
26  * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
27  * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
28  * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
29  * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
30  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
31  */
32
33 #ifdef HAVE_CONFIG_H
34 #  include <config.h>
35 #endif
36
37 #include <math.h>
38 #include <string.h>
39 #include "share/compat.h"
40 #include "private/bitmath.h"
41 #include "private/fixed.h"
42 #include "private/macros.h"
43 #include "FLAC/assert.h"
44
45 #ifdef local_abs
46 #undef local_abs
47 #endif
48 #define local_abs(x) ((unsigned)((x)<0? -(x) : (x)))
49
50 #ifdef FLAC__INTEGER_ONLY_LIBRARY
51 /* rbps stands for residual bits per sample
52  *
53  *             (ln(2) * err)
54  * rbps = log  (-----------)
55  *           2 (     n     )
56  */
57 static FLAC__fixedpoint local__compute_rbps_integerized(FLAC__uint32 err, FLAC__uint32 n)
58 {
59         FLAC__uint32 rbps;
60         unsigned bits; /* the number of bits required to represent a number */
61         int fracbits; /* the number of bits of rbps that comprise the fractional part */
62
63         FLAC__ASSERT(sizeof(rbps) == sizeof(FLAC__fixedpoint));
64         FLAC__ASSERT(err > 0);
65         FLAC__ASSERT(n > 0);
66
67         FLAC__ASSERT(n <= FLAC__MAX_BLOCK_SIZE);
68         if(err <= n)
69                 return 0;
70         /*
71          * The above two things tell us 1) n fits in 16 bits; 2) err/n > 1.
72          * These allow us later to know we won't lose too much precision in the
73          * fixed-point division (err<<fracbits)/n.
74          */
75
76         fracbits = (8*sizeof(err)) - (FLAC__bitmath_ilog2(err)+1);
77
78         err <<= fracbits;
79         err /= n;
80         /* err now holds err/n with fracbits fractional bits */
81
82         /*
83          * Whittle err down to 16 bits max.  16 significant bits is enough for
84          * our purposes.
85          */
86         FLAC__ASSERT(err > 0);
87         bits = FLAC__bitmath_ilog2(err)+1;
88         if(bits > 16) {
89                 err >>= (bits-16);
90                 fracbits -= (bits-16);
91         }
92         rbps = (FLAC__uint32)err;
93
94         /* Multiply by fixed-point version of ln(2), with 16 fractional bits */
95         rbps *= FLAC__FP_LN2;
96         fracbits += 16;
97         FLAC__ASSERT(fracbits >= 0);
98
99         /* FLAC__fixedpoint_log2 requires fracbits%4 to be 0 */
100         {
101                 const int f = fracbits & 3;
102                 if(f) {
103                         rbps >>= f;
104                         fracbits -= f;
105                 }
106         }
107
108         rbps = FLAC__fixedpoint_log2(rbps, fracbits, (unsigned)(-1));
109
110         if(rbps == 0)
111                 return 0;
112
113         /*
114          * The return value must have 16 fractional bits.  Since the whole part
115          * of the base-2 log of a 32 bit number must fit in 5 bits, and fracbits
116          * must be >= -3, these assertion allows us to be able to shift rbps
117          * left if necessary to get 16 fracbits without losing any bits of the
118          * whole part of rbps.
119          *
120          * There is a slight chance due to accumulated error that the whole part
121          * will require 6 bits, so we use 6 in the assertion.  Really though as
122          * long as it fits in 13 bits (32 - (16 - (-3))) we are fine.
123          */
124         FLAC__ASSERT((int)FLAC__bitmath_ilog2(rbps)+1 <= fracbits + 6);
125         FLAC__ASSERT(fracbits >= -3);
126
127         /* now shift the decimal point into place */
128         if(fracbits < 16)
129                 return rbps << (16-fracbits);
130         else if(fracbits > 16)
131                 return rbps >> (fracbits-16);
132         else
133                 return rbps;
134 }
135
136 static FLAC__fixedpoint local__compute_rbps_wide_integerized(FLAC__uint64 err, FLAC__uint32 n)
137 {
138         FLAC__uint32 rbps;
139         unsigned bits; /* the number of bits required to represent a number */
140         int fracbits; /* the number of bits of rbps that comprise the fractional part */
141
142         FLAC__ASSERT(sizeof(rbps) == sizeof(FLAC__fixedpoint));
143         FLAC__ASSERT(err > 0);
144         FLAC__ASSERT(n > 0);
145
146         FLAC__ASSERT(n <= FLAC__MAX_BLOCK_SIZE);
147         if(err <= n)
148                 return 0;
149         /*
150          * The above two things tell us 1) n fits in 16 bits; 2) err/n > 1.
151          * These allow us later to know we won't lose too much precision in the
152          * fixed-point division (err<<fracbits)/n.
153          */
154
155         fracbits = (8*sizeof(err)) - (FLAC__bitmath_ilog2_wide(err)+1);
156
157         err <<= fracbits;
158         err /= n;
159         /* err now holds err/n with fracbits fractional bits */
160
161         /*
162          * Whittle err down to 16 bits max.  16 significant bits is enough for
163          * our purposes.
164          */
165         FLAC__ASSERT(err > 0);
166         bits = FLAC__bitmath_ilog2_wide(err)+1;
167         if(bits > 16) {
168                 err >>= (bits-16);
169                 fracbits -= (bits-16);
170         }
171         rbps = (FLAC__uint32)err;
172
173         /* Multiply by fixed-point version of ln(2), with 16 fractional bits */
174         rbps *= FLAC__FP_LN2;
175         fracbits += 16;
176         FLAC__ASSERT(fracbits >= 0);
177
178         /* FLAC__fixedpoint_log2 requires fracbits%4 to be 0 */
179         {
180                 const int f = fracbits & 3;
181                 if(f) {
182                         rbps >>= f;
183                         fracbits -= f;
184                 }
185         }
186
187         rbps = FLAC__fixedpoint_log2(rbps, fracbits, (unsigned)(-1));
188
189         if(rbps == 0)
190                 return 0;
191
192         /*
193          * The return value must have 16 fractional bits.  Since the whole part
194          * of the base-2 log of a 32 bit number must fit in 5 bits, and fracbits
195          * must be >= -3, these assertion allows us to be able to shift rbps
196          * left if necessary to get 16 fracbits without losing any bits of the
197          * whole part of rbps.
198          *
199          * There is a slight chance due to accumulated error that the whole part
200          * will require 6 bits, so we use 6 in the assertion.  Really though as
201          * long as it fits in 13 bits (32 - (16 - (-3))) we are fine.
202          */
203         FLAC__ASSERT((int)FLAC__bitmath_ilog2(rbps)+1 <= fracbits + 6);
204         FLAC__ASSERT(fracbits >= -3);
205
206         /* now shift the decimal point into place */
207         if(fracbits < 16)
208                 return rbps << (16-fracbits);
209         else if(fracbits > 16)
210                 return rbps >> (fracbits-16);
211         else
212                 return rbps;
213 }
214 #endif
215
216 #ifndef FLAC__INTEGER_ONLY_LIBRARY
217 unsigned FLAC__fixed_compute_best_predictor(const FLAC__int32 data[], unsigned data_len, float residual_bits_per_sample[FLAC__MAX_FIXED_ORDER+1])
218 #else
219 unsigned FLAC__fixed_compute_best_predictor(const FLAC__int32 data[], unsigned data_len, FLAC__fixedpoint residual_bits_per_sample[FLAC__MAX_FIXED_ORDER+1])
220 #endif
221 {
222         FLAC__int32 last_error_0 = data[-1];
223         FLAC__int32 last_error_1 = data[-1] - data[-2];
224         FLAC__int32 last_error_2 = last_error_1 - (data[-2] - data[-3]);
225         FLAC__int32 last_error_3 = last_error_2 - (data[-2] - 2*data[-3] + data[-4]);
226         FLAC__int32 error, save;
227         FLAC__uint32 total_error_0 = 0, total_error_1 = 0, total_error_2 = 0, total_error_3 = 0, total_error_4 = 0;
228         unsigned i, order;
229
230         for(i = 0; i < data_len; i++) {
231                 error  = data[i]     ; total_error_0 += local_abs(error);                      save = error;
232                 error -= last_error_0; total_error_1 += local_abs(error); last_error_0 = save; save = error;
233                 error -= last_error_1; total_error_2 += local_abs(error); last_error_1 = save; save = error;
234                 error -= last_error_2; total_error_3 += local_abs(error); last_error_2 = save; save = error;
235                 error -= last_error_3; total_error_4 += local_abs(error); last_error_3 = save;
236         }
237
238         if(total_error_0 < flac_min(flac_min(flac_min(total_error_1, total_error_2), total_error_3), total_error_4))
239                 order = 0;
240         else if(total_error_1 < flac_min(flac_min(total_error_2, total_error_3), total_error_4))
241                 order = 1;
242         else if(total_error_2 < flac_min(total_error_3, total_error_4))
243                 order = 2;
244         else if(total_error_3 < total_error_4)
245                 order = 3;
246         else
247                 order = 4;
248
249         /* Estimate the expected number of bits per residual signal sample. */
250         /* 'total_error*' is linearly related to the variance of the residual */
251         /* signal, so we use it directly to compute E(|x|) */
252         FLAC__ASSERT(data_len > 0 || total_error_0 == 0);
253         FLAC__ASSERT(data_len > 0 || total_error_1 == 0);
254         FLAC__ASSERT(data_len > 0 || total_error_2 == 0);
255         FLAC__ASSERT(data_len > 0 || total_error_3 == 0);
256         FLAC__ASSERT(data_len > 0 || total_error_4 == 0);
257 #ifndef FLAC__INTEGER_ONLY_LIBRARY
258         residual_bits_per_sample[0] = (float)((total_error_0 > 0) ? log(M_LN2 * (double)total_error_0 / (double)data_len) / M_LN2 : 0.0);
259         residual_bits_per_sample[1] = (float)((total_error_1 > 0) ? log(M_LN2 * (double)total_error_1 / (double)data_len) / M_LN2 : 0.0);
260         residual_bits_per_sample[2] = (float)((total_error_2 > 0) ? log(M_LN2 * (double)total_error_2 / (double)data_len) / M_LN2 : 0.0);
261         residual_bits_per_sample[3] = (float)((total_error_3 > 0) ? log(M_LN2 * (double)total_error_3 / (double)data_len) / M_LN2 : 0.0);
262         residual_bits_per_sample[4] = (float)((total_error_4 > 0) ? log(M_LN2 * (double)total_error_4 / (double)data_len) / M_LN2 : 0.0);
263 #else
264         residual_bits_per_sample[0] = (total_error_0 > 0) ? local__compute_rbps_integerized(total_error_0, data_len) : 0;
265         residual_bits_per_sample[1] = (total_error_1 > 0) ? local__compute_rbps_integerized(total_error_1, data_len) : 0;
266         residual_bits_per_sample[2] = (total_error_2 > 0) ? local__compute_rbps_integerized(total_error_2, data_len) : 0;
267         residual_bits_per_sample[3] = (total_error_3 > 0) ? local__compute_rbps_integerized(total_error_3, data_len) : 0;
268         residual_bits_per_sample[4] = (total_error_4 > 0) ? local__compute_rbps_integerized(total_error_4, data_len) : 0;
269 #endif
270
271         return order;
272 }
273
274 #ifndef FLAC__INTEGER_ONLY_LIBRARY
275 unsigned FLAC__fixed_compute_best_predictor_wide(const FLAC__int32 data[], unsigned data_len, float residual_bits_per_sample[FLAC__MAX_FIXED_ORDER+1])
276 #else
277 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])
278 #endif
279 {
280         FLAC__int32 last_error_0 = data[-1];
281         FLAC__int32 last_error_1 = data[-1] - data[-2];
282         FLAC__int32 last_error_2 = last_error_1 - (data[-2] - data[-3]);
283         FLAC__int32 last_error_3 = last_error_2 - (data[-2] - 2*data[-3] + data[-4]);
284         FLAC__int32 error, save;
285         /* total_error_* are 64-bits to avoid overflow when encoding
286          * erratic signals when the bits-per-sample and blocksize are
287          * large.
288          */
289         FLAC__uint64 total_error_0 = 0, total_error_1 = 0, total_error_2 = 0, total_error_3 = 0, total_error_4 = 0;
290         unsigned i, order;
291
292         for(i = 0; i < data_len; i++) {
293                 error  = data[i]     ; total_error_0 += local_abs(error);                      save = error;
294                 error -= last_error_0; total_error_1 += local_abs(error); last_error_0 = save; save = error;
295                 error -= last_error_1; total_error_2 += local_abs(error); last_error_1 = save; save = error;
296                 error -= last_error_2; total_error_3 += local_abs(error); last_error_2 = save; save = error;
297                 error -= last_error_3; total_error_4 += local_abs(error); last_error_3 = save;
298         }
299
300         if(total_error_0 < flac_min(flac_min(flac_min(total_error_1, total_error_2), total_error_3), total_error_4))
301                 order = 0;
302         else if(total_error_1 < flac_min(flac_min(total_error_2, total_error_3), total_error_4))
303                 order = 1;
304         else if(total_error_2 < flac_min(total_error_3, total_error_4))
305                 order = 2;
306         else if(total_error_3 < total_error_4)
307                 order = 3;
308         else
309                 order = 4;
310
311         /* Estimate the expected number of bits per residual signal sample. */
312         /* 'total_error*' is linearly related to the variance of the residual */
313         /* signal, so we use it directly to compute E(|x|) */
314         FLAC__ASSERT(data_len > 0 || total_error_0 == 0);
315         FLAC__ASSERT(data_len > 0 || total_error_1 == 0);
316         FLAC__ASSERT(data_len > 0 || total_error_2 == 0);
317         FLAC__ASSERT(data_len > 0 || total_error_3 == 0);
318         FLAC__ASSERT(data_len > 0 || total_error_4 == 0);
319 #ifndef FLAC__INTEGER_ONLY_LIBRARY
320         residual_bits_per_sample[0] = (float)((total_error_0 > 0) ? log(M_LN2 * (double)total_error_0 / (double)data_len) / M_LN2 : 0.0);
321         residual_bits_per_sample[1] = (float)((total_error_1 > 0) ? log(M_LN2 * (double)total_error_1 / (double)data_len) / M_LN2 : 0.0);
322         residual_bits_per_sample[2] = (float)((total_error_2 > 0) ? log(M_LN2 * (double)total_error_2 / (double)data_len) / M_LN2 : 0.0);
323         residual_bits_per_sample[3] = (float)((total_error_3 > 0) ? log(M_LN2 * (double)total_error_3 / (double)data_len) / M_LN2 : 0.0);
324         residual_bits_per_sample[4] = (float)((total_error_4 > 0) ? log(M_LN2 * (double)total_error_4 / (double)data_len) / M_LN2 : 0.0);
325 #else
326         residual_bits_per_sample[0] = (total_error_0 > 0) ? local__compute_rbps_wide_integerized(total_error_0, data_len) : 0;
327         residual_bits_per_sample[1] = (total_error_1 > 0) ? local__compute_rbps_wide_integerized(total_error_1, data_len) : 0;
328         residual_bits_per_sample[2] = (total_error_2 > 0) ? local__compute_rbps_wide_integerized(total_error_2, data_len) : 0;
329         residual_bits_per_sample[3] = (total_error_3 > 0) ? local__compute_rbps_wide_integerized(total_error_3, data_len) : 0;
330         residual_bits_per_sample[4] = (total_error_4 > 0) ? local__compute_rbps_wide_integerized(total_error_4, data_len) : 0;
331 #endif
332
333         return order;
334 }
335
336 void FLAC__fixed_compute_residual(const FLAC__int32 data[], unsigned data_len, unsigned order, FLAC__int32 residual[])
337 {
338         const int idata_len = (int)data_len;
339         int i;
340
341         switch(order) {
342                 case 0:
343                         FLAC__ASSERT(sizeof(residual[0]) == sizeof(data[0]));
344                         memcpy(residual, data, sizeof(residual[0])*data_len);
345                         break;
346                 case 1:
347                         for(i = 0; i < idata_len; i++)
348                                 residual[i] = data[i] - data[i-1];
349                         break;
350                 case 2:
351                         for(i = 0; i < idata_len; i++)
352                                 residual[i] = data[i] - 2*data[i-1] + data[i-2];
353                         break;
354                 case 3:
355                         for(i = 0; i < idata_len; i++)
356                                 residual[i] = data[i] - 3*data[i-1] + 3*data[i-2] - data[i-3];
357                         break;
358                 case 4:
359                         for(i = 0; i < idata_len; i++)
360                                 residual[i] = data[i] - 4*data[i-1] + 6*data[i-2] - 4*data[i-3] + data[i-4];
361                         break;
362                 default:
363                         FLAC__ASSERT(0);
364         }
365 }
366
367 void FLAC__fixed_restore_signal(const FLAC__int32 residual[], unsigned data_len, unsigned order, FLAC__int32 data[])
368 {
369         int i, idata_len = (int)data_len;
370
371         switch(order) {
372                 case 0:
373                         FLAC__ASSERT(sizeof(residual[0]) == sizeof(data[0]));
374                         memcpy(data, residual, sizeof(residual[0])*data_len);
375                         break;
376                 case 1:
377                         for(i = 0; i < idata_len; i++)
378                                 data[i] = residual[i] + data[i-1];
379                         break;
380                 case 2:
381                         for(i = 0; i < idata_len; i++)
382                                 data[i] = residual[i] + 2*data[i-1] - data[i-2];
383                         break;
384                 case 3:
385                         for(i = 0; i < idata_len; i++)
386                                 data[i] = residual[i] + 3*data[i-1] - 3*data[i-2] + data[i-3];
387                         break;
388                 case 4:
389                         for(i = 0; i < idata_len; i++)
390                                 data[i] = residual[i] + 4*data[i-1] - 6*data[i-2] + 4*data[i-3] - data[i-4];
391                         break;
392                 default:
393                         FLAC__ASSERT(0);
394         }
395 }