1 /* libFLAC - Free Lossless Audio Codec library
2 * Copyright (C) 2000-2009 Josh Coalson
3 * Copyright (C) 2011-2023 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 "private/cpu.h"
38
39 #ifndef FLAC__INTEGER_ONLY_LIBRARY
40 #ifndef FLAC__NO_ASM
41 #if (defined FLAC__CPU_IA32 || defined FLAC__CPU_X86_64) && FLAC__HAS_X86INTRIN
42 #include "private/fixed.h"
43 #ifdef FLAC__AVX2_SUPPORTED
44
45 #include <immintrin.h>
46 #include <math.h>
47 #include "private/macros.h"
48 #include "share/compat.h"
49 #include "FLAC/assert.h"
50
51 #ifdef local_abs
52 #undef local_abs
53 #endif
54 #define local_abs(x) ((uint32_t)((x)<0? -(x) : (x)))
55
56 FLAC__SSE_TARGET("avx2")
FLAC__fixed_compute_best_predictor_wide_intrin_avx2(const FLAC__int32 data[],uint32_t data_len,float residual_bits_per_sample[FLAC__MAX_FIXED_ORDER+1])57 uint32_t FLAC__fixed_compute_best_predictor_wide_intrin_avx2(const FLAC__int32 data[], uint32_t data_len, float residual_bits_per_sample[FLAC__MAX_FIXED_ORDER + 1])
58 {
59 FLAC__uint64 total_error_0, total_error_1, total_error_2, total_error_3, total_error_4;
60 FLAC__int32 i, data_len_int;
61 uint32_t order;
62 __m256i total_err0, total_err1, total_err2, total_err3, total_err4;
63 __m256i prev_err0, prev_err1, prev_err2, prev_err3;
64 __m256i tempA, tempB, bitmask;
65 FLAC__int64 data_scalar[4];
66 FLAC__int64 prev_err0_scalar[4];
67 FLAC__int64 prev_err1_scalar[4];
68 FLAC__int64 prev_err2_scalar[4];
69 FLAC__int64 prev_err3_scalar[4];
70 total_err0 = _mm256_setzero_si256();
71 total_err1 = _mm256_setzero_si256();
72 total_err2 = _mm256_setzero_si256();
73 total_err3 = _mm256_setzero_si256();
74 total_err4 = _mm256_setzero_si256();
75 data_len_int = data_len;
76
77 for(i = 0; i < 4; i++){
78 prev_err0_scalar[i] = data[-1+i*(data_len_int/4)];
79 prev_err1_scalar[i] = data[-1+i*(data_len_int/4)] - data[-2+i*(data_len_int/4)];
80 prev_err2_scalar[i] = prev_err1_scalar[i] - (data[-2+i*(data_len_int/4)] - data[-3+i*(data_len_int/4)]);
81 prev_err3_scalar[i] = prev_err2_scalar[i] - (data[-2+i*(data_len_int/4)] - 2*data[-3+i*(data_len_int/4)] + data[-4+i*(data_len_int/4)]);
82 }
83 prev_err0 = _mm256_loadu_si256((const __m256i*)(void*)prev_err0_scalar);
84 prev_err1 = _mm256_loadu_si256((const __m256i*)(void*)prev_err1_scalar);
85 prev_err2 = _mm256_loadu_si256((const __m256i*)(void*)prev_err2_scalar);
86 prev_err3 = _mm256_loadu_si256((const __m256i*)(void*)prev_err3_scalar);
87 for(i = 0; i < data_len_int / 4; i++){
88 data_scalar[0] = data[i];
89 data_scalar[1] = data[i+data_len/4];
90 data_scalar[2] = data[i+2*data_len/4];
91 data_scalar[3] = data[i+3*data_len/4];
92 tempA = _mm256_loadu_si256((const __m256i*)(void*)data_scalar);
93 /* Next three intrinsics calculate tempB as abs of tempA */
94 bitmask = _mm256_cmpgt_epi64(_mm256_set1_epi64x(0), tempA);
95 tempB = _mm256_xor_si256(tempA, bitmask);
96 tempB = _mm256_sub_epi64(tempB, bitmask);
97 total_err0 = _mm256_add_epi64(total_err0,tempB);
98 tempB = _mm256_sub_epi64(tempA,prev_err0);
99 prev_err0 = tempA;
100 /* Next three intrinsics calculate tempA as abs of tempB */
101 bitmask = _mm256_cmpgt_epi64(_mm256_set1_epi64x(0), tempB);
102 tempA = _mm256_xor_si256(tempB, bitmask);
103 tempA = _mm256_sub_epi64(tempA, bitmask);
104 total_err1 = _mm256_add_epi64(total_err1,tempA);
105 tempA = _mm256_sub_epi64(tempB,prev_err1);
106 prev_err1 = tempB;
107 /* Next three intrinsics calculate tempB as abs of tempA */
108 bitmask = _mm256_cmpgt_epi64(_mm256_set1_epi64x(0), tempA);
109 tempB = _mm256_xor_si256(tempA, bitmask);
110 tempB = _mm256_sub_epi64(tempB, bitmask);
111 total_err2 = _mm256_add_epi64(total_err2,tempB);
112 tempB = _mm256_sub_epi64(tempA,prev_err2);
113 prev_err2 = tempA;
114 /* Next three intrinsics calculate tempA as abs of tempB */
115 bitmask = _mm256_cmpgt_epi64(_mm256_set1_epi64x(0), tempB);
116 tempA = _mm256_xor_si256(tempB, bitmask);
117 tempA = _mm256_sub_epi64(tempA, bitmask);
118 total_err3 = _mm256_add_epi64(total_err3,tempA);
119 tempA = _mm256_sub_epi64(tempB,prev_err3);
120 prev_err3 = tempB;
121 /* Next three intrinsics calculate tempB as abs of tempA */
122 bitmask = _mm256_cmpgt_epi64(_mm256_set1_epi64x(0), tempA);
123 tempB = _mm256_xor_si256(tempA, bitmask);
124 tempB = _mm256_sub_epi64(tempB, bitmask);
125 total_err4 = _mm256_add_epi64(total_err4,tempB);
126 }
127 _mm256_storeu_si256((__m256i*)(void*)data_scalar,total_err0);
128 total_error_0 = data_scalar[0] + data_scalar[1] + data_scalar[2] + data_scalar[3];
129 _mm256_storeu_si256((__m256i*)(void*)data_scalar,total_err1);
130 total_error_1 = data_scalar[0] + data_scalar[1] + data_scalar[2] + data_scalar[3];
131 _mm256_storeu_si256((__m256i*)(void*)data_scalar,total_err2);
132 total_error_2 = data_scalar[0] + data_scalar[1] + data_scalar[2] + data_scalar[3];
133 _mm256_storeu_si256((__m256i*)(void*)data_scalar,total_err3);
134 total_error_3 = data_scalar[0] + data_scalar[1] + data_scalar[2] + data_scalar[3];
135 _mm256_storeu_si256((__m256i*)(void*)data_scalar,total_err4);
136 total_error_4 = data_scalar[0] + data_scalar[1] + data_scalar[2] + data_scalar[3];
137
138 /* Ignore the remainder, we're ignore the first few samples too */
139
140 /* prefer lower order */
141 if(total_error_0 <= flac_min(flac_min(flac_min(total_error_1, total_error_2), total_error_3), total_error_4))
142 order = 0;
143 else if(total_error_1 <= flac_min(flac_min(total_error_2, total_error_3), total_error_4))
144 order = 1;
145 else if(total_error_2 <= flac_min(total_error_3, total_error_4))
146 order = 2;
147 else if(total_error_3 <= total_error_4)
148 order = 3;
149 else
150 order = 4;
151
152 /* Estimate the expected number of bits per residual signal sample. */
153 /* 'total_error*' is linearly related to the variance of the residual */
154 /* signal, so we use it directly to compute E(|x|) */
155 FLAC__ASSERT(data_len > 0 || total_error_0 == 0);
156 FLAC__ASSERT(data_len > 0 || total_error_1 == 0);
157 FLAC__ASSERT(data_len > 0 || total_error_2 == 0);
158 FLAC__ASSERT(data_len > 0 || total_error_3 == 0);
159 FLAC__ASSERT(data_len > 0 || total_error_4 == 0);
160
161 residual_bits_per_sample[0] = (float)((total_error_0 > 0) ? log(M_LN2 * (double)total_error_0 / (double)data_len) / M_LN2 : 0.0);
162 residual_bits_per_sample[1] = (float)((total_error_1 > 0) ? log(M_LN2 * (double)total_error_1 / (double)data_len) / M_LN2 : 0.0);
163 residual_bits_per_sample[2] = (float)((total_error_2 > 0) ? log(M_LN2 * (double)total_error_2 / (double)data_len) / M_LN2 : 0.0);
164 residual_bits_per_sample[3] = (float)((total_error_3 > 0) ? log(M_LN2 * (double)total_error_3 / (double)data_len) / M_LN2 : 0.0);
165 residual_bits_per_sample[4] = (float)((total_error_4 > 0) ? log(M_LN2 * (double)total_error_4 / (double)data_len) / M_LN2 : 0.0);
166
167 return order;
168 }
169
170 #ifdef local_abs64
171 #undef local_abs64
172 #endif
173 #define local_abs64(x) ((uint64_t)((x)<0? -(x) : (x)))
174
175 #define CHECK_ORDER_IS_VALID(macro_order) \
176 if(shadow_error_##macro_order <= INT32_MAX) { \
177 if(total_error_##macro_order < smallest_error) { \
178 order = macro_order; \
179 smallest_error = total_error_##macro_order ; \
180 } \
181 residual_bits_per_sample[ macro_order ] = (float)((total_error_0 > 0) ? log(M_LN2 * (double)total_error_0 / (double)data_len) / M_LN2 : 0.0); \
182 } \
183 else \
184 residual_bits_per_sample[ macro_order ] = 34.0f;
185
186 FLAC__SSE_TARGET("avx2")
FLAC__fixed_compute_best_predictor_limit_residual_intrin_avx2(const FLAC__int32 data[],uint32_t data_len,float residual_bits_per_sample[FLAC__MAX_FIXED_ORDER+1])187 uint32_t FLAC__fixed_compute_best_predictor_limit_residual_intrin_avx2(const FLAC__int32 data[], uint32_t data_len, float residual_bits_per_sample[FLAC__MAX_FIXED_ORDER + 1])
188 {
189 FLAC__uint64 total_error_0 = 0, total_error_1 = 0, total_error_2 = 0, total_error_3 = 0, total_error_4 = 0, smallest_error = UINT64_MAX;
190 FLAC__uint64 shadow_error_0 = 0, shadow_error_1 = 0, shadow_error_2 = 0, shadow_error_3 = 0, shadow_error_4 = 0;
191 FLAC__uint64 error_0, error_1, error_2, error_3, error_4;
192 FLAC__int32 i, data_len_int;
193 uint32_t order = 0;
194 __m256i total_err0, total_err1, total_err2, total_err3, total_err4;
195 __m256i shadow_err0, shadow_err1, shadow_err2, shadow_err3, shadow_err4;
196 __m256i prev_err0, prev_err1, prev_err2, prev_err3;
197 __m256i tempA, tempB, bitmask;
198 FLAC__int64 data_scalar[4];
199 FLAC__int64 prev_err0_scalar[4];
200 FLAC__int64 prev_err1_scalar[4];
201 FLAC__int64 prev_err2_scalar[4];
202 FLAC__int64 prev_err3_scalar[4];
203 total_err0 = _mm256_setzero_si256();
204 total_err1 = _mm256_setzero_si256();
205 total_err2 = _mm256_setzero_si256();
206 total_err3 = _mm256_setzero_si256();
207 total_err4 = _mm256_setzero_si256();
208 shadow_err0 = _mm256_setzero_si256();
209 shadow_err1 = _mm256_setzero_si256();
210 shadow_err2 = _mm256_setzero_si256();
211 shadow_err3 = _mm256_setzero_si256();
212 shadow_err4 = _mm256_setzero_si256();
213 data_len_int = data_len;
214
215 /* First take care of preceding samples */
216 for(i = -4; i < 0; i++) {
217 error_0 = local_abs64((FLAC__int64)data[i]);
218 error_1 = (i > -4) ? local_abs64((FLAC__int64)data[i] - data[i-1]) : 0 ;
219 error_2 = (i > -3) ? local_abs64((FLAC__int64)data[i] - 2 * (FLAC__int64)data[i-1] + data[i-2]) : 0;
220 error_3 = (i > -2) ? local_abs64((FLAC__int64)data[i] - 3 * (FLAC__int64)data[i-1] + 3 * (FLAC__int64)data[i-2] - data[i-3]) : 0;
221
222 total_error_0 += error_0;
223 total_error_1 += error_1;
224 total_error_2 += error_2;
225 total_error_3 += error_3;
226
227 shadow_error_0 |= error_0;
228 shadow_error_1 |= error_1;
229 shadow_error_2 |= error_2;
230 shadow_error_3 |= error_3;
231 }
232
233 for(i = 0; i < 4; i++){
234 prev_err0_scalar[i] = data[-1+i*(data_len_int/4)];
235 prev_err1_scalar[i] = (FLAC__int64)(data[-1+i*(data_len_int/4)]) - data[-2+i*(data_len_int/4)];
236 prev_err2_scalar[i] = prev_err1_scalar[i] - ((FLAC__int64)(data[-2+i*(data_len_int/4)]) - data[-3+i*(data_len_int/4)]);
237 prev_err3_scalar[i] = prev_err2_scalar[i] - ((FLAC__int64)(data[-2+i*(data_len_int/4)]) - 2*(FLAC__int64)(data[-3+i*(data_len_int/4)]) + data[-4+i*(data_len_int/4)]);
238 }
239 prev_err0 = _mm256_loadu_si256((const __m256i*)(void*)prev_err0_scalar);
240 prev_err1 = _mm256_loadu_si256((const __m256i*)(void*)prev_err1_scalar);
241 prev_err2 = _mm256_loadu_si256((const __m256i*)(void*)prev_err2_scalar);
242 prev_err3 = _mm256_loadu_si256((const __m256i*)(void*)prev_err3_scalar);
243 for(i = 0; i < data_len_int / 4; i++){
244 data_scalar[0] = data[i];
245 data_scalar[1] = data[i+data_len/4];
246 data_scalar[2] = data[i+2*data_len/4];
247 data_scalar[3] = data[i+3*data_len/4];
248 tempA = _mm256_loadu_si256((const __m256i*)(void*)data_scalar);
249 /* Next three intrinsics calculate tempB as abs of tempA */
250 bitmask = _mm256_cmpgt_epi64(_mm256_set1_epi64x(0), tempA);
251 tempB = _mm256_xor_si256(tempA, bitmask);
252 tempB = _mm256_sub_epi64(tempB, bitmask);
253 total_err0 = _mm256_add_epi64(total_err0,tempB);
254 shadow_err0 = _mm256_or_si256(shadow_err0,tempB);
255 tempB = _mm256_sub_epi64(tempA,prev_err0);
256 prev_err0 = tempA;
257 /* Next three intrinsics calculate tempA as abs of tempB */
258 bitmask = _mm256_cmpgt_epi64(_mm256_set1_epi64x(0), tempB);
259 tempA = _mm256_xor_si256(tempB, bitmask);
260 tempA = _mm256_sub_epi64(tempA, bitmask);
261 total_err1 = _mm256_add_epi64(total_err1,tempA);
262 shadow_err1 = _mm256_or_si256(shadow_err1,tempA);
263 tempA = _mm256_sub_epi64(tempB,prev_err1);
264 prev_err1 = tempB;
265 /* Next three intrinsics calculate tempB as abs of tempA */
266 bitmask = _mm256_cmpgt_epi64(_mm256_set1_epi64x(0), tempA);
267 tempB = _mm256_xor_si256(tempA, bitmask);
268 tempB = _mm256_sub_epi64(tempB, bitmask);
269 total_err2 = _mm256_add_epi64(total_err2,tempB);
270 shadow_err2 = _mm256_or_si256(shadow_err2,tempB);
271 tempB = _mm256_sub_epi64(tempA,prev_err2);
272 prev_err2 = tempA;
273 /* Next three intrinsics calculate tempA as abs of tempB */
274 bitmask = _mm256_cmpgt_epi64(_mm256_set1_epi64x(0), tempB);
275 tempA = _mm256_xor_si256(tempB, bitmask);
276 tempA = _mm256_sub_epi64(tempA, bitmask);
277 total_err3 = _mm256_add_epi64(total_err3,tempA);
278 shadow_err3 = _mm256_or_si256(shadow_err3,tempA);
279 tempA = _mm256_sub_epi64(tempB,prev_err3);
280 prev_err3 = tempB;
281 /* Next three intrinsics calculate tempB as abs of tempA */
282 bitmask = _mm256_cmpgt_epi64(_mm256_set1_epi64x(0), tempA);
283 tempB = _mm256_xor_si256(tempA, bitmask);
284 tempB = _mm256_sub_epi64(tempB, bitmask);
285 total_err4 = _mm256_add_epi64(total_err4,tempB);
286 shadow_err4 = _mm256_or_si256(shadow_err4,tempB);
287 }
288 _mm256_storeu_si256((__m256i*)(void*)data_scalar,total_err0);
289 total_error_0 += data_scalar[0] + data_scalar[1] + data_scalar[2] + data_scalar[3];
290 _mm256_storeu_si256((__m256i*)(void*)data_scalar,total_err1);
291 total_error_1 += data_scalar[0] + data_scalar[1] + data_scalar[2] + data_scalar[3];
292 _mm256_storeu_si256((__m256i*)(void*)data_scalar,total_err2);
293 total_error_2 += data_scalar[0] + data_scalar[1] + data_scalar[2] + data_scalar[3];
294 _mm256_storeu_si256((__m256i*)(void*)data_scalar,total_err3);
295 total_error_3 += data_scalar[0] + data_scalar[1] + data_scalar[2] + data_scalar[3];
296 _mm256_storeu_si256((__m256i*)(void*)data_scalar,total_err4);
297 total_error_4 += data_scalar[0] + data_scalar[1] + data_scalar[2] + data_scalar[3];
298 _mm256_storeu_si256((__m256i*)(void*)data_scalar,shadow_err0);
299 shadow_error_0 |= data_scalar[0] | data_scalar[1] | data_scalar[2] | data_scalar[3];
300 _mm256_storeu_si256((__m256i*)(void*)data_scalar,shadow_err1);
301 shadow_error_1 |= data_scalar[0] | data_scalar[1] | data_scalar[2] | data_scalar[3];
302 _mm256_storeu_si256((__m256i*)(void*)data_scalar,shadow_err2);
303 shadow_error_2 |= data_scalar[0] | data_scalar[1] | data_scalar[2] | data_scalar[3];
304 _mm256_storeu_si256((__m256i*)(void*)data_scalar,shadow_err3);
305 shadow_error_3 |= data_scalar[0] | data_scalar[1] | data_scalar[2] | data_scalar[3];
306 _mm256_storeu_si256((__m256i*)(void*)data_scalar,shadow_err4);
307 shadow_error_4 |= data_scalar[0] | data_scalar[1] | data_scalar[2] | data_scalar[3];
308
309 /* Take care of remaining sample */
310 for(i = (data_len/4)*4; i < data_len_int; i++) {
311 error_0 = local_abs64((FLAC__int64)data[i]);
312 error_1 = local_abs64((FLAC__int64)data[i] - data[i-1]);
313 error_2 = local_abs64((FLAC__int64)data[i] - 2 * (FLAC__int64)data[i-1] + data[i-2]);
314 error_3 = local_abs64((FLAC__int64)data[i] - 3 * (FLAC__int64)data[i-1] + 3 * (FLAC__int64)data[i-2] - data[i-3]);
315 error_4 = local_abs64((FLAC__int64)data[i] - 4 * (FLAC__int64)data[i-1] + 6 * (FLAC__int64)data[i-2] - 4 * (FLAC__int64)data[i-3] + data[i-4]);
316
317 total_error_0 += error_0;
318 total_error_1 += error_1;
319 total_error_2 += error_2;
320 total_error_3 += error_3;
321 total_error_4 += error_4;
322
323 shadow_error_0 |= error_0;
324 shadow_error_1 |= error_1;
325 shadow_error_2 |= error_2;
326 shadow_error_3 |= error_3;
327 shadow_error_4 |= error_4;
328 }
329
330
331 CHECK_ORDER_IS_VALID(0);
332 CHECK_ORDER_IS_VALID(1);
333 CHECK_ORDER_IS_VALID(2);
334 CHECK_ORDER_IS_VALID(3);
335 CHECK_ORDER_IS_VALID(4);
336
337 return order;
338 }
339
340 #endif /* FLAC__AVX2_SUPPORTED */
341 #endif /* (FLAC__CPU_IA32 || FLAC__CPU_X86_64) && FLAC__HAS_X86INTRIN */
342 #endif /* FLAC__NO_ASM */
343 #endif /* FLAC__INTEGER_ONLY_LIBRARY */
344