be more precise in long long int literal specification for ANSI compilers (gcc3 requi...
[flac.git] / src / libFLAC / fixed.c
1 /* libFLAC - Free Lossless Audio Codec library
2  * Copyright (C) 2000,2001,2002,2003,2004  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 #include <math.h>
33 #include "private/fixed.h"
34 #include "FLAC/assert.h"
35
36 #ifndef M_LN2
37 /* math.h in VC++ doesn't seem to have this (how Microsoft is that?) */
38 #define M_LN2 0.69314718055994530942
39 #endif
40
41 #ifdef min
42 #undef min
43 #endif
44 #define min(x,y) ((x) < (y)? (x) : (y))
45
46 #ifdef local_abs
47 #undef local_abs
48 #endif
49 #define local_abs(x) ((unsigned)((x)<0? -(x) : (x)))
50
51 unsigned FLAC__fixed_compute_best_predictor(const FLAC__int32 data[], unsigned data_len, FLAC__real residual_bits_per_sample[FLAC__MAX_FIXED_ORDER+1])
52 {
53         FLAC__int32 last_error_0 = data[-1];
54         FLAC__int32 last_error_1 = data[-1] - data[-2];
55         FLAC__int32 last_error_2 = last_error_1 - (data[-2] - data[-3]);
56         FLAC__int32 last_error_3 = last_error_2 - (data[-2] - 2*data[-3] + data[-4]);
57         FLAC__int32 error, save;
58         FLAC__uint32 total_error_0 = 0, total_error_1 = 0, total_error_2 = 0, total_error_3 = 0, total_error_4 = 0;
59         unsigned i, order;
60
61         for(i = 0; i < data_len; i++) {
62                 error  = data[i]     ; total_error_0 += local_abs(error);                      save = error;
63                 error -= last_error_0; total_error_1 += local_abs(error); last_error_0 = save; save = error;
64                 error -= last_error_1; total_error_2 += local_abs(error); last_error_1 = save; save = error;
65                 error -= last_error_2; total_error_3 += local_abs(error); last_error_2 = save; save = error;
66                 error -= last_error_3; total_error_4 += local_abs(error); last_error_3 = save;
67         }
68
69         if(total_error_0 < min(min(min(total_error_1, total_error_2), total_error_3), total_error_4))
70                 order = 0;
71         else if(total_error_1 < min(min(total_error_2, total_error_3), total_error_4))
72                 order = 1;
73         else if(total_error_2 < min(total_error_3, total_error_4))
74                 order = 2;
75         else if(total_error_3 < total_error_4)
76                 order = 3;
77         else
78                 order = 4;
79
80         /* Estimate the expected number of bits per residual signal sample. */
81         /* 'total_error*' is linearly related to the variance of the residual */
82         /* signal, so we use it directly to compute E(|x|) */
83         FLAC__ASSERT(data_len > 0 || total_error_0 == 0);
84         FLAC__ASSERT(data_len > 0 || total_error_1 == 0);
85         FLAC__ASSERT(data_len > 0 || total_error_2 == 0);
86         FLAC__ASSERT(data_len > 0 || total_error_3 == 0);
87         FLAC__ASSERT(data_len > 0 || total_error_4 == 0);
88         residual_bits_per_sample[0] = (FLAC__real)((total_error_0 > 0) ? log(M_LN2 * (double)total_error_0 / (double)data_len) / M_LN2 : 0.0);
89         residual_bits_per_sample[1] = (FLAC__real)((total_error_1 > 0) ? log(M_LN2 * (double)total_error_1 / (double)data_len) / M_LN2 : 0.0);
90         residual_bits_per_sample[2] = (FLAC__real)((total_error_2 > 0) ? log(M_LN2 * (double)total_error_2 / (double)data_len) / M_LN2 : 0.0);
91         residual_bits_per_sample[3] = (FLAC__real)((total_error_3 > 0) ? log(M_LN2 * (double)total_error_3 / (double)data_len) / M_LN2 : 0.0);
92         residual_bits_per_sample[4] = (FLAC__real)((total_error_4 > 0) ? log(M_LN2 * (double)total_error_4 / (double)data_len) / M_LN2 : 0.0);
93
94         return order;
95 }
96
97 unsigned FLAC__fixed_compute_best_predictor_wide(const FLAC__int32 data[], unsigned data_len, FLAC__real residual_bits_per_sample[FLAC__MAX_FIXED_ORDER+1])
98 {
99         FLAC__int32 last_error_0 = data[-1];
100         FLAC__int32 last_error_1 = data[-1] - data[-2];
101         FLAC__int32 last_error_2 = last_error_1 - (data[-2] - data[-3]);
102         FLAC__int32 last_error_3 = last_error_2 - (data[-2] - 2*data[-3] + data[-4]);
103         FLAC__int32 error, save;
104         /* total_error_* are 64-bits to avoid overflow when encoding
105          * erratic signals when the bits-per-sample and blocksize are
106          * large.
107          */
108         FLAC__uint64 total_error_0 = 0, total_error_1 = 0, total_error_2 = 0, total_error_3 = 0, total_error_4 = 0;
109         unsigned i, order;
110
111         for(i = 0; i < data_len; i++) {
112                 error  = data[i]     ; total_error_0 += local_abs(error);                      save = error;
113                 error -= last_error_0; total_error_1 += local_abs(error); last_error_0 = save; save = error;
114                 error -= last_error_1; total_error_2 += local_abs(error); last_error_1 = save; save = error;
115                 error -= last_error_2; total_error_3 += local_abs(error); last_error_2 = save; save = error;
116                 error -= last_error_3; total_error_4 += local_abs(error); last_error_3 = save;
117         }
118
119         if(total_error_0 < min(min(min(total_error_1, total_error_2), total_error_3), total_error_4))
120                 order = 0;
121         else if(total_error_1 < min(min(total_error_2, total_error_3), total_error_4))
122                 order = 1;
123         else if(total_error_2 < min(total_error_3, total_error_4))
124                 order = 2;
125         else if(total_error_3 < total_error_4)
126                 order = 3;
127         else
128                 order = 4;
129
130         /* Estimate the expected number of bits per residual signal sample. */
131         /* 'total_error*' is linearly related to the variance of the residual */
132         /* signal, so we use it directly to compute E(|x|) */
133         FLAC__ASSERT(data_len > 0 || total_error_0 == 0);
134         FLAC__ASSERT(data_len > 0 || total_error_1 == 0);
135         FLAC__ASSERT(data_len > 0 || total_error_2 == 0);
136         FLAC__ASSERT(data_len > 0 || total_error_3 == 0);
137         FLAC__ASSERT(data_len > 0 || total_error_4 == 0);
138 #if defined _MSC_VER || defined __MINGW32__
139         /* with VC++ you have to spoon feed it the casting */
140         residual_bits_per_sample[0] = (FLAC__real)((total_error_0 > 0) ? log(M_LN2 * (double)(FLAC__int64)total_error_0 / (double)data_len) / M_LN2 : 0.0);
141         residual_bits_per_sample[1] = (FLAC__real)((total_error_1 > 0) ? log(M_LN2 * (double)(FLAC__int64)total_error_1 / (double)data_len) / M_LN2 : 0.0);
142         residual_bits_per_sample[2] = (FLAC__real)((total_error_2 > 0) ? log(M_LN2 * (double)(FLAC__int64)total_error_2 / (double)data_len) / M_LN2 : 0.0);
143         residual_bits_per_sample[3] = (FLAC__real)((total_error_3 > 0) ? log(M_LN2 * (double)(FLAC__int64)total_error_3 / (double)data_len) / M_LN2 : 0.0);
144         residual_bits_per_sample[4] = (FLAC__real)((total_error_4 > 0) ? log(M_LN2 * (double)(FLAC__int64)total_error_4 / (double)data_len) / M_LN2 : 0.0);
145 #else
146         residual_bits_per_sample[0] = (FLAC__real)((total_error_0 > 0) ? log(M_LN2 * (double)total_error_0 / (double)data_len) / M_LN2 : 0.0);
147         residual_bits_per_sample[1] = (FLAC__real)((total_error_1 > 0) ? log(M_LN2 * (double)total_error_1 / (double)data_len) / M_LN2 : 0.0);
148         residual_bits_per_sample[2] = (FLAC__real)((total_error_2 > 0) ? log(M_LN2 * (double)total_error_2 / (double)data_len) / M_LN2 : 0.0);
149         residual_bits_per_sample[3] = (FLAC__real)((total_error_3 > 0) ? log(M_LN2 * (double)total_error_3 / (double)data_len) / M_LN2 : 0.0);
150         residual_bits_per_sample[4] = (FLAC__real)((total_error_4 > 0) ? log(M_LN2 * (double)total_error_4 / (double)data_len) / M_LN2 : 0.0);
151 #endif
152
153         return order;
154 }
155
156 void FLAC__fixed_compute_residual(const FLAC__int32 data[], unsigned data_len, unsigned order, FLAC__int32 residual[])
157 {
158         const int idata_len = (int)data_len;
159         int i;
160
161         switch(order) {
162                 case 0:
163                         for(i = 0; i < idata_len; i++) {
164                                 residual[i] = data[i];
165                         }
166                         break;
167                 case 1:
168                         for(i = 0; i < idata_len; i++) {
169                                 residual[i] = data[i] - data[i-1];
170                         }
171                         break;
172                 case 2:
173                         for(i = 0; i < idata_len; i++) {
174                                 /* == data[i] - 2*data[i-1] + data[i-2] */
175                                 residual[i] = data[i] - (data[i-1] << 1) + data[i-2];
176                         }
177                         break;
178                 case 3:
179                         for(i = 0; i < idata_len; i++) {
180                                 /* == data[i] - 3*data[i-1] + 3*data[i-2] - data[i-3] */
181                                 residual[i] = data[i] - (((data[i-1]-data[i-2])<<1) + (data[i-1]-data[i-2])) - data[i-3];
182                         }
183                         break;
184                 case 4:
185                         for(i = 0; i < idata_len; i++) {
186                                 /* == data[i] - 4*data[i-1] + 6*data[i-2] - 4*data[i-3] + data[i-4] */
187                                 residual[i] = data[i] - ((data[i-1]+data[i-3])<<2) + ((data[i-2]<<2) + (data[i-2]<<1)) + data[i-4];
188                         }
189                         break;
190                 default:
191                         FLAC__ASSERT(0);
192         }
193 }
194
195 void FLAC__fixed_restore_signal(const FLAC__int32 residual[], unsigned data_len, unsigned order, FLAC__int32 data[])
196 {
197         int i, idata_len = (int)data_len;
198
199         switch(order) {
200                 case 0:
201                         for(i = 0; i < idata_len; i++) {
202                                 data[i] = residual[i];
203                         }
204                         break;
205                 case 1:
206                         for(i = 0; i < idata_len; i++) {
207                                 data[i] = residual[i] + data[i-1];
208                         }
209                         break;
210                 case 2:
211                         for(i = 0; i < idata_len; i++) {
212                                 /* == residual[i] + 2*data[i-1] - data[i-2] */
213                                 data[i] = residual[i] + (data[i-1]<<1) - data[i-2];
214                         }
215                         break;
216                 case 3:
217                         for(i = 0; i < idata_len; i++) {
218                                 /* residual[i] + 3*data[i-1] - 3*data[i-2]) + data[i-3] */
219                                 data[i] = residual[i] + (((data[i-1]-data[i-2])<<1) + (data[i-1]-data[i-2])) + data[i-3];
220                         }
221                         break;
222                 case 4:
223                         for(i = 0; i < idata_len; i++) {
224                                 /* == residual[i] + 4*data[i-1] - 6*data[i-2] + 4*data[i-3] - data[i-4] */
225                                 data[i] = residual[i] + ((data[i-1]+data[i-3])<<2) - ((data[i-2]<<2) + (data[i-2]<<1)) - data[i-4];
226                         }
227                         break;
228                 default:
229                         FLAC__ASSERT(0);
230         }
231 }