jidctint-neon.c (35406B)
1 /* 2 * jidctint-neon.c - accurate integer IDCT (Arm Neon) 3 * 4 * Copyright (C) 2020, Arm Limited. All Rights Reserved. 5 * Copyright (C) 2020, D. R. Commander. All Rights Reserved. 6 * 7 * This software is provided 'as-is', without any express or implied 8 * warranty. In no event will the authors be held liable for any damages 9 * arising from the use of this software. 10 * 11 * Permission is granted to anyone to use this software for any purpose, 12 * including commercial applications, and to alter it and redistribute it 13 * freely, subject to the following restrictions: 14 * 15 * 1. The origin of this software must not be misrepresented; you must not 16 * claim that you wrote the original software. If you use this software 17 * in a product, an acknowledgment in the product documentation would be 18 * appreciated but is not required. 19 * 2. Altered source versions must be plainly marked as such, and must not be 20 * misrepresented as being the original software. 21 * 3. This notice may not be removed or altered from any source distribution. 22 */ 23 24 #define JPEG_INTERNALS 25 #include "../../jinclude.h" 26 #include "../../jpeglib.h" 27 #include "../../jsimd.h" 28 #include "../../jdct.h" 29 #include "../../jsimddct.h" 30 #include "../jsimd.h" 31 #include "align.h" 32 #include "neon-compat.h" 33 34 #include <arm_neon.h> 35 36 37 #define CONST_BITS 13 38 #define PASS1_BITS 2 39 40 #define DESCALE_P1 (CONST_BITS - PASS1_BITS) 41 #define DESCALE_P2 (CONST_BITS + PASS1_BITS + 3) 42 43 /* The computation of the inverse DCT requires the use of constants known at 44 * compile time. Scaled integer constants are used to avoid floating-point 45 * arithmetic: 46 * 0.298631336 = 2446 * 2^-13 47 * 0.390180644 = 3196 * 2^-13 48 * 0.541196100 = 4433 * 2^-13 49 * 0.765366865 = 6270 * 2^-13 50 * 0.899976223 = 7373 * 2^-13 51 * 1.175875602 = 9633 * 2^-13 52 * 1.501321110 = 12299 * 2^-13 53 * 1.847759065 = 15137 * 2^-13 54 * 1.961570560 = 16069 * 2^-13 55 * 2.053119869 = 16819 * 2^-13 56 * 2.562915447 = 20995 * 2^-13 57 * 3.072711026 = 25172 * 2^-13 58 */ 59 60 #define F_0_298 2446 61 #define F_0_390 3196 62 #define F_0_541 4433 63 #define F_0_765 6270 64 #define F_0_899 7373 65 #define F_1_175 9633 66 #define F_1_501 12299 67 #define F_1_847 15137 68 #define F_1_961 16069 69 #define F_2_053 16819 70 #define F_2_562 20995 71 #define F_3_072 25172 72 73 #define F_1_175_MINUS_1_961 (F_1_175 - F_1_961) 74 #define F_1_175_MINUS_0_390 (F_1_175 - F_0_390) 75 #define F_0_541_MINUS_1_847 (F_0_541 - F_1_847) 76 #define F_3_072_MINUS_2_562 (F_3_072 - F_2_562) 77 #define F_0_298_MINUS_0_899 (F_0_298 - F_0_899) 78 #define F_1_501_MINUS_0_899 (F_1_501 - F_0_899) 79 #define F_2_053_MINUS_2_562 (F_2_053 - F_2_562) 80 #define F_0_541_PLUS_0_765 (F_0_541 + F_0_765) 81 82 83 ALIGN(16) static const int16_t jsimd_idct_islow_neon_consts[] = { 84 F_0_899, F_0_541, 85 F_2_562, F_0_298_MINUS_0_899, 86 F_1_501_MINUS_0_899, F_2_053_MINUS_2_562, 87 F_0_541_PLUS_0_765, F_1_175, 88 F_1_175_MINUS_0_390, F_0_541_MINUS_1_847, 89 F_3_072_MINUS_2_562, F_1_175_MINUS_1_961, 90 0, 0, 0, 0 91 }; 92 93 94 /* Forward declaration of regular and sparse IDCT helper functions */ 95 96 static INLINE void jsimd_idct_islow_pass1_regular(int16x4_t row0, 97 int16x4_t row1, 98 int16x4_t row2, 99 int16x4_t row3, 100 int16x4_t row4, 101 int16x4_t row5, 102 int16x4_t row6, 103 int16x4_t row7, 104 int16x4_t quant_row0, 105 int16x4_t quant_row1, 106 int16x4_t quant_row2, 107 int16x4_t quant_row3, 108 int16x4_t quant_row4, 109 int16x4_t quant_row5, 110 int16x4_t quant_row6, 111 int16x4_t quant_row7, 112 int16_t *workspace_1, 113 int16_t *workspace_2); 114 115 static INLINE void jsimd_idct_islow_pass1_sparse(int16x4_t row0, 116 int16x4_t row1, 117 int16x4_t row2, 118 int16x4_t row3, 119 int16x4_t quant_row0, 120 int16x4_t quant_row1, 121 int16x4_t quant_row2, 122 int16x4_t quant_row3, 123 int16_t *workspace_1, 124 int16_t *workspace_2); 125 126 static INLINE void jsimd_idct_islow_pass2_regular(int16_t *workspace, 127 JSAMPARRAY output_buf, 128 JDIMENSION output_col, 129 unsigned buf_offset); 130 131 static INLINE void jsimd_idct_islow_pass2_sparse(int16_t *workspace, 132 JSAMPARRAY output_buf, 133 JDIMENSION output_col, 134 unsigned buf_offset); 135 136 137 /* Perform dequantization and inverse DCT on one block of coefficients. For 138 * reference, the C implementation (jpeg_idct_slow()) can be found in 139 * jidctint.c. 140 * 141 * Optimization techniques used for fast data access: 142 * 143 * In each pass, the inverse DCT is computed for the left and right 4x8 halves 144 * of the DCT block. This avoids spilling due to register pressure, and the 145 * increased granularity allows for an optimized calculation depending on the 146 * values of the DCT coefficients. Between passes, intermediate data is stored 147 * in 4x8 workspace buffers. 148 * 149 * Transposing the 8x8 DCT block after each pass can be achieved by transposing 150 * each of the four 4x4 quadrants and swapping quadrants 1 and 2 (refer to the 151 * diagram below.) Swapping quadrants is cheap, since the second pass can just 152 * swap the workspace buffer pointers. 153 * 154 * +-------+-------+ +-------+-------+ 155 * | | | | | | 156 * | 0 | 1 | | 0 | 2 | 157 * | | | transpose | | | 158 * +-------+-------+ ------> +-------+-------+ 159 * | | | | | | 160 * | 2 | 3 | | 1 | 3 | 161 * | | | | | | 162 * +-------+-------+ +-------+-------+ 163 * 164 * Optimization techniques used to accelerate the inverse DCT calculation: 165 * 166 * In a DCT coefficient block, the coefficients are increasingly likely to be 0 167 * as you move diagonally from top left to bottom right. If whole rows of 168 * coefficients are 0, then the inverse DCT calculation can be simplified. On 169 * the first pass of the inverse DCT, we test for three special cases before 170 * defaulting to a full "regular" inverse DCT: 171 * 172 * 1) Coefficients in rows 4-7 are all zero. In this case, we perform a 173 * "sparse" simplified inverse DCT on rows 0-3. 174 * 2) AC coefficients (rows 1-7) are all zero. In this case, the inverse DCT 175 * result is equal to the dequantized DC coefficients. 176 * 3) AC and DC coefficients are all zero. In this case, the inverse DCT 177 * result is all zero. For the left 4x8 half, this is handled identically 178 * to Case 2 above. For the right 4x8 half, we do no work and signal that 179 * the "sparse" algorithm is required for the second pass. 180 * 181 * In the second pass, only a single special case is tested: whether the AC and 182 * DC coefficients were all zero in the right 4x8 block during the first pass 183 * (refer to Case 3 above.) If this is the case, then a "sparse" variant of 184 * the second pass is performed for both the left and right halves of the DCT 185 * block. (The transposition after the first pass means that the right 4x8 186 * block during the first pass becomes rows 4-7 during the second pass.) 187 */ 188 189 void jsimd_idct_islow_neon(void *dct_table, JCOEFPTR coef_block, 190 JSAMPARRAY output_buf, JDIMENSION output_col) 191 { 192 ISLOW_MULT_TYPE *quantptr = dct_table; 193 194 int16_t workspace_l[8 * DCTSIZE / 2]; 195 int16_t workspace_r[8 * DCTSIZE / 2]; 196 197 /* Compute IDCT first pass on left 4x8 coefficient block. */ 198 199 /* Load DCT coefficients in left 4x8 block. */ 200 int16x4_t row0 = vld1_s16(coef_block + 0 * DCTSIZE); 201 int16x4_t row1 = vld1_s16(coef_block + 1 * DCTSIZE); 202 int16x4_t row2 = vld1_s16(coef_block + 2 * DCTSIZE); 203 int16x4_t row3 = vld1_s16(coef_block + 3 * DCTSIZE); 204 int16x4_t row4 = vld1_s16(coef_block + 4 * DCTSIZE); 205 int16x4_t row5 = vld1_s16(coef_block + 5 * DCTSIZE); 206 int16x4_t row6 = vld1_s16(coef_block + 6 * DCTSIZE); 207 int16x4_t row7 = vld1_s16(coef_block + 7 * DCTSIZE); 208 209 /* Load quantization table for left 4x8 block. */ 210 int16x4_t quant_row0 = vld1_s16(quantptr + 0 * DCTSIZE); 211 int16x4_t quant_row1 = vld1_s16(quantptr + 1 * DCTSIZE); 212 int16x4_t quant_row2 = vld1_s16(quantptr + 2 * DCTSIZE); 213 int16x4_t quant_row3 = vld1_s16(quantptr + 3 * DCTSIZE); 214 int16x4_t quant_row4 = vld1_s16(quantptr + 4 * DCTSIZE); 215 int16x4_t quant_row5 = vld1_s16(quantptr + 5 * DCTSIZE); 216 int16x4_t quant_row6 = vld1_s16(quantptr + 6 * DCTSIZE); 217 int16x4_t quant_row7 = vld1_s16(quantptr + 7 * DCTSIZE); 218 219 /* Construct bitmap to test if DCT coefficients in left 4x8 block are 0. */ 220 int16x4_t bitmap = vorr_s16(row7, row6); 221 bitmap = vorr_s16(bitmap, row5); 222 bitmap = vorr_s16(bitmap, row4); 223 int64_t bitmap_rows_4567 = vget_lane_s64(vreinterpret_s64_s16(bitmap), 0); 224 225 if (bitmap_rows_4567 == 0) { 226 bitmap = vorr_s16(bitmap, row3); 227 bitmap = vorr_s16(bitmap, row2); 228 bitmap = vorr_s16(bitmap, row1); 229 int64_t left_ac_bitmap = vget_lane_s64(vreinterpret_s64_s16(bitmap), 0); 230 231 if (left_ac_bitmap == 0) { 232 int16x4_t dcval = vshl_n_s16(vmul_s16(row0, quant_row0), PASS1_BITS); 233 int16x4x4_t quadrant = { { dcval, dcval, dcval, dcval } }; 234 /* Store 4x4 blocks to workspace, transposing in the process. */ 235 vst4_s16(workspace_l, quadrant); 236 vst4_s16(workspace_r, quadrant); 237 } else { 238 jsimd_idct_islow_pass1_sparse(row0, row1, row2, row3, quant_row0, 239 quant_row1, quant_row2, quant_row3, 240 workspace_l, workspace_r); 241 } 242 } else { 243 jsimd_idct_islow_pass1_regular(row0, row1, row2, row3, row4, row5, 244 row6, row7, quant_row0, quant_row1, 245 quant_row2, quant_row3, quant_row4, 246 quant_row5, quant_row6, quant_row7, 247 workspace_l, workspace_r); 248 } 249 250 /* Compute IDCT first pass on right 4x8 coefficient block. */ 251 252 /* Load DCT coefficients in right 4x8 block. */ 253 row0 = vld1_s16(coef_block + 0 * DCTSIZE + 4); 254 row1 = vld1_s16(coef_block + 1 * DCTSIZE + 4); 255 row2 = vld1_s16(coef_block + 2 * DCTSIZE + 4); 256 row3 = vld1_s16(coef_block + 3 * DCTSIZE + 4); 257 row4 = vld1_s16(coef_block + 4 * DCTSIZE + 4); 258 row5 = vld1_s16(coef_block + 5 * DCTSIZE + 4); 259 row6 = vld1_s16(coef_block + 6 * DCTSIZE + 4); 260 row7 = vld1_s16(coef_block + 7 * DCTSIZE + 4); 261 262 /* Load quantization table for right 4x8 block. */ 263 quant_row0 = vld1_s16(quantptr + 0 * DCTSIZE + 4); 264 quant_row1 = vld1_s16(quantptr + 1 * DCTSIZE + 4); 265 quant_row2 = vld1_s16(quantptr + 2 * DCTSIZE + 4); 266 quant_row3 = vld1_s16(quantptr + 3 * DCTSIZE + 4); 267 quant_row4 = vld1_s16(quantptr + 4 * DCTSIZE + 4); 268 quant_row5 = vld1_s16(quantptr + 5 * DCTSIZE + 4); 269 quant_row6 = vld1_s16(quantptr + 6 * DCTSIZE + 4); 270 quant_row7 = vld1_s16(quantptr + 7 * DCTSIZE + 4); 271 272 /* Construct bitmap to test if DCT coefficients in right 4x8 block are 0. */ 273 bitmap = vorr_s16(row7, row6); 274 bitmap = vorr_s16(bitmap, row5); 275 bitmap = vorr_s16(bitmap, row4); 276 bitmap_rows_4567 = vget_lane_s64(vreinterpret_s64_s16(bitmap), 0); 277 bitmap = vorr_s16(bitmap, row3); 278 bitmap = vorr_s16(bitmap, row2); 279 bitmap = vorr_s16(bitmap, row1); 280 int64_t right_ac_bitmap = vget_lane_s64(vreinterpret_s64_s16(bitmap), 0); 281 282 /* If this remains non-zero, a "regular" second pass will be performed. */ 283 int64_t right_ac_dc_bitmap = 1; 284 285 if (right_ac_bitmap == 0) { 286 bitmap = vorr_s16(bitmap, row0); 287 right_ac_dc_bitmap = vget_lane_s64(vreinterpret_s64_s16(bitmap), 0); 288 289 if (right_ac_dc_bitmap != 0) { 290 int16x4_t dcval = vshl_n_s16(vmul_s16(row0, quant_row0), PASS1_BITS); 291 int16x4x4_t quadrant = { { dcval, dcval, dcval, dcval } }; 292 /* Store 4x4 blocks to workspace, transposing in the process. */ 293 vst4_s16(workspace_l + 4 * DCTSIZE / 2, quadrant); 294 vst4_s16(workspace_r + 4 * DCTSIZE / 2, quadrant); 295 } 296 } else { 297 if (bitmap_rows_4567 == 0) { 298 jsimd_idct_islow_pass1_sparse(row0, row1, row2, row3, quant_row0, 299 quant_row1, quant_row2, quant_row3, 300 workspace_l + 4 * DCTSIZE / 2, 301 workspace_r + 4 * DCTSIZE / 2); 302 } else { 303 jsimd_idct_islow_pass1_regular(row0, row1, row2, row3, row4, row5, 304 row6, row7, quant_row0, quant_row1, 305 quant_row2, quant_row3, quant_row4, 306 quant_row5, quant_row6, quant_row7, 307 workspace_l + 4 * DCTSIZE / 2, 308 workspace_r + 4 * DCTSIZE / 2); 309 } 310 } 311 312 /* Second pass: compute IDCT on rows in workspace. */ 313 314 /* If all coefficients in right 4x8 block are 0, use "sparse" second pass. */ 315 if (right_ac_dc_bitmap == 0) { 316 jsimd_idct_islow_pass2_sparse(workspace_l, output_buf, output_col, 0); 317 jsimd_idct_islow_pass2_sparse(workspace_r, output_buf, output_col, 4); 318 } else { 319 jsimd_idct_islow_pass2_regular(workspace_l, output_buf, output_col, 0); 320 jsimd_idct_islow_pass2_regular(workspace_r, output_buf, output_col, 4); 321 } 322 } 323 324 325 /* Perform dequantization and the first pass of the accurate inverse DCT on a 326 * 4x8 block of coefficients. (To process the full 8x8 DCT block, this 327 * function-- or some other optimized variant-- needs to be called for both the 328 * left and right 4x8 blocks.) 329 * 330 * This "regular" version assumes that no optimization can be made to the IDCT 331 * calculation, since no useful set of AC coefficients is all 0. 332 * 333 * The original C implementation of the accurate IDCT (jpeg_idct_slow()) can be 334 * found in jidctint.c. Algorithmic changes made here are documented inline. 335 */ 336 337 static INLINE void jsimd_idct_islow_pass1_regular(int16x4_t row0, 338 int16x4_t row1, 339 int16x4_t row2, 340 int16x4_t row3, 341 int16x4_t row4, 342 int16x4_t row5, 343 int16x4_t row6, 344 int16x4_t row7, 345 int16x4_t quant_row0, 346 int16x4_t quant_row1, 347 int16x4_t quant_row2, 348 int16x4_t quant_row3, 349 int16x4_t quant_row4, 350 int16x4_t quant_row5, 351 int16x4_t quant_row6, 352 int16x4_t quant_row7, 353 int16_t *workspace_1, 354 int16_t *workspace_2) 355 { 356 /* Load constants for IDCT computation. */ 357 #ifdef HAVE_VLD1_S16_X3 358 const int16x4x3_t consts = vld1_s16_x3(jsimd_idct_islow_neon_consts); 359 #else 360 const int16x4_t consts1 = vld1_s16(jsimd_idct_islow_neon_consts); 361 const int16x4_t consts2 = vld1_s16(jsimd_idct_islow_neon_consts + 4); 362 const int16x4_t consts3 = vld1_s16(jsimd_idct_islow_neon_consts + 8); 363 const int16x4x3_t consts = { { consts1, consts2, consts3 } }; 364 #endif 365 366 /* Even part */ 367 int16x4_t z2_s16 = vmul_s16(row2, quant_row2); 368 int16x4_t z3_s16 = vmul_s16(row6, quant_row6); 369 370 int32x4_t tmp2 = vmull_lane_s16(z2_s16, consts.val[0], 1); 371 int32x4_t tmp3 = vmull_lane_s16(z2_s16, consts.val[1], 2); 372 tmp2 = vmlal_lane_s16(tmp2, z3_s16, consts.val[2], 1); 373 tmp3 = vmlal_lane_s16(tmp3, z3_s16, consts.val[0], 1); 374 375 z2_s16 = vmul_s16(row0, quant_row0); 376 z3_s16 = vmul_s16(row4, quant_row4); 377 378 int32x4_t tmp0 = vshll_n_s16(vadd_s16(z2_s16, z3_s16), CONST_BITS); 379 int32x4_t tmp1 = vshll_n_s16(vsub_s16(z2_s16, z3_s16), CONST_BITS); 380 381 int32x4_t tmp10 = vaddq_s32(tmp0, tmp3); 382 int32x4_t tmp13 = vsubq_s32(tmp0, tmp3); 383 int32x4_t tmp11 = vaddq_s32(tmp1, tmp2); 384 int32x4_t tmp12 = vsubq_s32(tmp1, tmp2); 385 386 /* Odd part */ 387 int16x4_t tmp0_s16 = vmul_s16(row7, quant_row7); 388 int16x4_t tmp1_s16 = vmul_s16(row5, quant_row5); 389 int16x4_t tmp2_s16 = vmul_s16(row3, quant_row3); 390 int16x4_t tmp3_s16 = vmul_s16(row1, quant_row1); 391 392 z3_s16 = vadd_s16(tmp0_s16, tmp2_s16); 393 int16x4_t z4_s16 = vadd_s16(tmp1_s16, tmp3_s16); 394 395 /* Implementation as per jpeg_idct_islow() in jidctint.c: 396 * z5 = (z3 + z4) * 1.175875602; 397 * z3 = z3 * -1.961570560; z4 = z4 * -0.390180644; 398 * z3 += z5; z4 += z5; 399 * 400 * This implementation: 401 * z3 = z3 * (1.175875602 - 1.961570560) + z4 * 1.175875602; 402 * z4 = z3 * 1.175875602 + z4 * (1.175875602 - 0.390180644); 403 */ 404 405 int32x4_t z3 = vmull_lane_s16(z3_s16, consts.val[2], 3); 406 int32x4_t z4 = vmull_lane_s16(z3_s16, consts.val[1], 3); 407 z3 = vmlal_lane_s16(z3, z4_s16, consts.val[1], 3); 408 z4 = vmlal_lane_s16(z4, z4_s16, consts.val[2], 0); 409 410 /* Implementation as per jpeg_idct_islow() in jidctint.c: 411 * z1 = tmp0 + tmp3; z2 = tmp1 + tmp2; 412 * tmp0 = tmp0 * 0.298631336; tmp1 = tmp1 * 2.053119869; 413 * tmp2 = tmp2 * 3.072711026; tmp3 = tmp3 * 1.501321110; 414 * z1 = z1 * -0.899976223; z2 = z2 * -2.562915447; 415 * tmp0 += z1 + z3; tmp1 += z2 + z4; 416 * tmp2 += z2 + z3; tmp3 += z1 + z4; 417 * 418 * This implementation: 419 * tmp0 = tmp0 * (0.298631336 - 0.899976223) + tmp3 * -0.899976223; 420 * tmp1 = tmp1 * (2.053119869 - 2.562915447) + tmp2 * -2.562915447; 421 * tmp2 = tmp1 * -2.562915447 + tmp2 * (3.072711026 - 2.562915447); 422 * tmp3 = tmp0 * -0.899976223 + tmp3 * (1.501321110 - 0.899976223); 423 * tmp0 += z3; tmp1 += z4; 424 * tmp2 += z3; tmp3 += z4; 425 */ 426 427 tmp0 = vmull_lane_s16(tmp0_s16, consts.val[0], 3); 428 tmp1 = vmull_lane_s16(tmp1_s16, consts.val[1], 1); 429 tmp2 = vmull_lane_s16(tmp2_s16, consts.val[2], 2); 430 tmp3 = vmull_lane_s16(tmp3_s16, consts.val[1], 0); 431 432 tmp0 = vmlsl_lane_s16(tmp0, tmp3_s16, consts.val[0], 0); 433 tmp1 = vmlsl_lane_s16(tmp1, tmp2_s16, consts.val[0], 2); 434 tmp2 = vmlsl_lane_s16(tmp2, tmp1_s16, consts.val[0], 2); 435 tmp3 = vmlsl_lane_s16(tmp3, tmp0_s16, consts.val[0], 0); 436 437 tmp0 = vaddq_s32(tmp0, z3); 438 tmp1 = vaddq_s32(tmp1, z4); 439 tmp2 = vaddq_s32(tmp2, z3); 440 tmp3 = vaddq_s32(tmp3, z4); 441 442 /* Final output stage: descale and narrow to 16-bit. */ 443 int16x4x4_t rows_0123 = { { 444 vrshrn_n_s32(vaddq_s32(tmp10, tmp3), DESCALE_P1), 445 vrshrn_n_s32(vaddq_s32(tmp11, tmp2), DESCALE_P1), 446 vrshrn_n_s32(vaddq_s32(tmp12, tmp1), DESCALE_P1), 447 vrshrn_n_s32(vaddq_s32(tmp13, tmp0), DESCALE_P1) 448 } }; 449 int16x4x4_t rows_4567 = { { 450 vrshrn_n_s32(vsubq_s32(tmp13, tmp0), DESCALE_P1), 451 vrshrn_n_s32(vsubq_s32(tmp12, tmp1), DESCALE_P1), 452 vrshrn_n_s32(vsubq_s32(tmp11, tmp2), DESCALE_P1), 453 vrshrn_n_s32(vsubq_s32(tmp10, tmp3), DESCALE_P1) 454 } }; 455 456 /* Store 4x4 blocks to the intermediate workspace, ready for the second pass. 457 * (VST4 transposes the blocks. We need to operate on rows in the next 458 * pass.) 459 */ 460 vst4_s16(workspace_1, rows_0123); 461 vst4_s16(workspace_2, rows_4567); 462 } 463 464 465 /* Perform dequantization and the first pass of the accurate inverse DCT on a 466 * 4x8 block of coefficients. 467 * 468 * This "sparse" version assumes that the AC coefficients in rows 4-7 are all 469 * 0. This simplifies the IDCT calculation, accelerating overall performance. 470 */ 471 472 static INLINE void jsimd_idct_islow_pass1_sparse(int16x4_t row0, 473 int16x4_t row1, 474 int16x4_t row2, 475 int16x4_t row3, 476 int16x4_t quant_row0, 477 int16x4_t quant_row1, 478 int16x4_t quant_row2, 479 int16x4_t quant_row3, 480 int16_t *workspace_1, 481 int16_t *workspace_2) 482 { 483 /* Load constants for IDCT computation. */ 484 #ifdef HAVE_VLD1_S16_X3 485 const int16x4x3_t consts = vld1_s16_x3(jsimd_idct_islow_neon_consts); 486 #else 487 const int16x4_t consts1 = vld1_s16(jsimd_idct_islow_neon_consts); 488 const int16x4_t consts2 = vld1_s16(jsimd_idct_islow_neon_consts + 4); 489 const int16x4_t consts3 = vld1_s16(jsimd_idct_islow_neon_consts + 8); 490 const int16x4x3_t consts = { { consts1, consts2, consts3 } }; 491 #endif 492 493 /* Even part (z3 is all 0) */ 494 int16x4_t z2_s16 = vmul_s16(row2, quant_row2); 495 496 int32x4_t tmp2 = vmull_lane_s16(z2_s16, consts.val[0], 1); 497 int32x4_t tmp3 = vmull_lane_s16(z2_s16, consts.val[1], 2); 498 499 z2_s16 = vmul_s16(row0, quant_row0); 500 int32x4_t tmp0 = vshll_n_s16(z2_s16, CONST_BITS); 501 int32x4_t tmp1 = vshll_n_s16(z2_s16, CONST_BITS); 502 503 int32x4_t tmp10 = vaddq_s32(tmp0, tmp3); 504 int32x4_t tmp13 = vsubq_s32(tmp0, tmp3); 505 int32x4_t tmp11 = vaddq_s32(tmp1, tmp2); 506 int32x4_t tmp12 = vsubq_s32(tmp1, tmp2); 507 508 /* Odd part (tmp0 and tmp1 are both all 0) */ 509 int16x4_t tmp2_s16 = vmul_s16(row3, quant_row3); 510 int16x4_t tmp3_s16 = vmul_s16(row1, quant_row1); 511 512 int16x4_t z3_s16 = tmp2_s16; 513 int16x4_t z4_s16 = tmp3_s16; 514 515 int32x4_t z3 = vmull_lane_s16(z3_s16, consts.val[2], 3); 516 int32x4_t z4 = vmull_lane_s16(z3_s16, consts.val[1], 3); 517 z3 = vmlal_lane_s16(z3, z4_s16, consts.val[1], 3); 518 z4 = vmlal_lane_s16(z4, z4_s16, consts.val[2], 0); 519 520 tmp0 = vmlsl_lane_s16(z3, tmp3_s16, consts.val[0], 0); 521 tmp1 = vmlsl_lane_s16(z4, tmp2_s16, consts.val[0], 2); 522 tmp2 = vmlal_lane_s16(z3, tmp2_s16, consts.val[2], 2); 523 tmp3 = vmlal_lane_s16(z4, tmp3_s16, consts.val[1], 0); 524 525 /* Final output stage: descale and narrow to 16-bit. */ 526 int16x4x4_t rows_0123 = { { 527 vrshrn_n_s32(vaddq_s32(tmp10, tmp3), DESCALE_P1), 528 vrshrn_n_s32(vaddq_s32(tmp11, tmp2), DESCALE_P1), 529 vrshrn_n_s32(vaddq_s32(tmp12, tmp1), DESCALE_P1), 530 vrshrn_n_s32(vaddq_s32(tmp13, tmp0), DESCALE_P1) 531 } }; 532 int16x4x4_t rows_4567 = { { 533 vrshrn_n_s32(vsubq_s32(tmp13, tmp0), DESCALE_P1), 534 vrshrn_n_s32(vsubq_s32(tmp12, tmp1), DESCALE_P1), 535 vrshrn_n_s32(vsubq_s32(tmp11, tmp2), DESCALE_P1), 536 vrshrn_n_s32(vsubq_s32(tmp10, tmp3), DESCALE_P1) 537 } }; 538 539 /* Store 4x4 blocks to the intermediate workspace, ready for the second pass. 540 * (VST4 transposes the blocks. We need to operate on rows in the next 541 * pass.) 542 */ 543 vst4_s16(workspace_1, rows_0123); 544 vst4_s16(workspace_2, rows_4567); 545 } 546 547 548 /* Perform the second pass of the accurate inverse DCT on a 4x8 block of 549 * coefficients. (To process the full 8x8 DCT block, this function-- or some 550 * other optimized variant-- needs to be called for both the right and left 4x8 551 * blocks.) 552 * 553 * This "regular" version assumes that no optimization can be made to the IDCT 554 * calculation, since no useful set of coefficient values are all 0 after the 555 * first pass. 556 * 557 * Again, the original C implementation of the accurate IDCT (jpeg_idct_slow()) 558 * can be found in jidctint.c. Algorithmic changes made here are documented 559 * inline. 560 */ 561 562 static INLINE void jsimd_idct_islow_pass2_regular(int16_t *workspace, 563 JSAMPARRAY output_buf, 564 JDIMENSION output_col, 565 unsigned buf_offset) 566 { 567 /* Load constants for IDCT computation. */ 568 #ifdef HAVE_VLD1_S16_X3 569 const int16x4x3_t consts = vld1_s16_x3(jsimd_idct_islow_neon_consts); 570 #else 571 const int16x4_t consts1 = vld1_s16(jsimd_idct_islow_neon_consts); 572 const int16x4_t consts2 = vld1_s16(jsimd_idct_islow_neon_consts + 4); 573 const int16x4_t consts3 = vld1_s16(jsimd_idct_islow_neon_consts + 8); 574 const int16x4x3_t consts = { { consts1, consts2, consts3 } }; 575 #endif 576 577 /* Even part */ 578 int16x4_t z2_s16 = vld1_s16(workspace + 2 * DCTSIZE / 2); 579 int16x4_t z3_s16 = vld1_s16(workspace + 6 * DCTSIZE / 2); 580 581 int32x4_t tmp2 = vmull_lane_s16(z2_s16, consts.val[0], 1); 582 int32x4_t tmp3 = vmull_lane_s16(z2_s16, consts.val[1], 2); 583 tmp2 = vmlal_lane_s16(tmp2, z3_s16, consts.val[2], 1); 584 tmp3 = vmlal_lane_s16(tmp3, z3_s16, consts.val[0], 1); 585 586 z2_s16 = vld1_s16(workspace + 0 * DCTSIZE / 2); 587 z3_s16 = vld1_s16(workspace + 4 * DCTSIZE / 2); 588 589 int32x4_t tmp0 = vshll_n_s16(vadd_s16(z2_s16, z3_s16), CONST_BITS); 590 int32x4_t tmp1 = vshll_n_s16(vsub_s16(z2_s16, z3_s16), CONST_BITS); 591 592 int32x4_t tmp10 = vaddq_s32(tmp0, tmp3); 593 int32x4_t tmp13 = vsubq_s32(tmp0, tmp3); 594 int32x4_t tmp11 = vaddq_s32(tmp1, tmp2); 595 int32x4_t tmp12 = vsubq_s32(tmp1, tmp2); 596 597 /* Odd part */ 598 int16x4_t tmp0_s16 = vld1_s16(workspace + 7 * DCTSIZE / 2); 599 int16x4_t tmp1_s16 = vld1_s16(workspace + 5 * DCTSIZE / 2); 600 int16x4_t tmp2_s16 = vld1_s16(workspace + 3 * DCTSIZE / 2); 601 int16x4_t tmp3_s16 = vld1_s16(workspace + 1 * DCTSIZE / 2); 602 603 z3_s16 = vadd_s16(tmp0_s16, tmp2_s16); 604 int16x4_t z4_s16 = vadd_s16(tmp1_s16, tmp3_s16); 605 606 /* Implementation as per jpeg_idct_islow() in jidctint.c: 607 * z5 = (z3 + z4) * 1.175875602; 608 * z3 = z3 * -1.961570560; z4 = z4 * -0.390180644; 609 * z3 += z5; z4 += z5; 610 * 611 * This implementation: 612 * z3 = z3 * (1.175875602 - 1.961570560) + z4 * 1.175875602; 613 * z4 = z3 * 1.175875602 + z4 * (1.175875602 - 0.390180644); 614 */ 615 616 int32x4_t z3 = vmull_lane_s16(z3_s16, consts.val[2], 3); 617 int32x4_t z4 = vmull_lane_s16(z3_s16, consts.val[1], 3); 618 z3 = vmlal_lane_s16(z3, z4_s16, consts.val[1], 3); 619 z4 = vmlal_lane_s16(z4, z4_s16, consts.val[2], 0); 620 621 /* Implementation as per jpeg_idct_islow() in jidctint.c: 622 * z1 = tmp0 + tmp3; z2 = tmp1 + tmp2; 623 * tmp0 = tmp0 * 0.298631336; tmp1 = tmp1 * 2.053119869; 624 * tmp2 = tmp2 * 3.072711026; tmp3 = tmp3 * 1.501321110; 625 * z1 = z1 * -0.899976223; z2 = z2 * -2.562915447; 626 * tmp0 += z1 + z3; tmp1 += z2 + z4; 627 * tmp2 += z2 + z3; tmp3 += z1 + z4; 628 * 629 * This implementation: 630 * tmp0 = tmp0 * (0.298631336 - 0.899976223) + tmp3 * -0.899976223; 631 * tmp1 = tmp1 * (2.053119869 - 2.562915447) + tmp2 * -2.562915447; 632 * tmp2 = tmp1 * -2.562915447 + tmp2 * (3.072711026 - 2.562915447); 633 * tmp3 = tmp0 * -0.899976223 + tmp3 * (1.501321110 - 0.899976223); 634 * tmp0 += z3; tmp1 += z4; 635 * tmp2 += z3; tmp3 += z4; 636 */ 637 638 tmp0 = vmull_lane_s16(tmp0_s16, consts.val[0], 3); 639 tmp1 = vmull_lane_s16(tmp1_s16, consts.val[1], 1); 640 tmp2 = vmull_lane_s16(tmp2_s16, consts.val[2], 2); 641 tmp3 = vmull_lane_s16(tmp3_s16, consts.val[1], 0); 642 643 tmp0 = vmlsl_lane_s16(tmp0, tmp3_s16, consts.val[0], 0); 644 tmp1 = vmlsl_lane_s16(tmp1, tmp2_s16, consts.val[0], 2); 645 tmp2 = vmlsl_lane_s16(tmp2, tmp1_s16, consts.val[0], 2); 646 tmp3 = vmlsl_lane_s16(tmp3, tmp0_s16, consts.val[0], 0); 647 648 tmp0 = vaddq_s32(tmp0, z3); 649 tmp1 = vaddq_s32(tmp1, z4); 650 tmp2 = vaddq_s32(tmp2, z3); 651 tmp3 = vaddq_s32(tmp3, z4); 652 653 /* Final output stage: descale and narrow to 16-bit. */ 654 int16x8_t cols_02_s16 = vcombine_s16(vaddhn_s32(tmp10, tmp3), 655 vaddhn_s32(tmp12, tmp1)); 656 int16x8_t cols_13_s16 = vcombine_s16(vaddhn_s32(tmp11, tmp2), 657 vaddhn_s32(tmp13, tmp0)); 658 int16x8_t cols_46_s16 = vcombine_s16(vsubhn_s32(tmp13, tmp0), 659 vsubhn_s32(tmp11, tmp2)); 660 int16x8_t cols_57_s16 = vcombine_s16(vsubhn_s32(tmp12, tmp1), 661 vsubhn_s32(tmp10, tmp3)); 662 /* Descale and narrow to 8-bit. */ 663 int8x8_t cols_02_s8 = vqrshrn_n_s16(cols_02_s16, DESCALE_P2 - 16); 664 int8x8_t cols_13_s8 = vqrshrn_n_s16(cols_13_s16, DESCALE_P2 - 16); 665 int8x8_t cols_46_s8 = vqrshrn_n_s16(cols_46_s16, DESCALE_P2 - 16); 666 int8x8_t cols_57_s8 = vqrshrn_n_s16(cols_57_s16, DESCALE_P2 - 16); 667 /* Clamp to range [0-255]. */ 668 uint8x8_t cols_02_u8 = vadd_u8(vreinterpret_u8_s8(cols_02_s8), 669 vdup_n_u8(CENTERJSAMPLE)); 670 uint8x8_t cols_13_u8 = vadd_u8(vreinterpret_u8_s8(cols_13_s8), 671 vdup_n_u8(CENTERJSAMPLE)); 672 uint8x8_t cols_46_u8 = vadd_u8(vreinterpret_u8_s8(cols_46_s8), 673 vdup_n_u8(CENTERJSAMPLE)); 674 uint8x8_t cols_57_u8 = vadd_u8(vreinterpret_u8_s8(cols_57_s8), 675 vdup_n_u8(CENTERJSAMPLE)); 676 677 /* Transpose 4x8 block and store to memory. (Zipping adjacent columns 678 * together allows us to store 16-bit elements.) 679 */ 680 uint8x8x2_t cols_01_23 = vzip_u8(cols_02_u8, cols_13_u8); 681 uint8x8x2_t cols_45_67 = vzip_u8(cols_46_u8, cols_57_u8); 682 uint16x4x4_t cols_01_23_45_67 = { { 683 vreinterpret_u16_u8(cols_01_23.val[0]), 684 vreinterpret_u16_u8(cols_01_23.val[1]), 685 vreinterpret_u16_u8(cols_45_67.val[0]), 686 vreinterpret_u16_u8(cols_45_67.val[1]) 687 } }; 688 689 JSAMPROW outptr0 = output_buf[buf_offset + 0] + output_col; 690 JSAMPROW outptr1 = output_buf[buf_offset + 1] + output_col; 691 JSAMPROW outptr2 = output_buf[buf_offset + 2] + output_col; 692 JSAMPROW outptr3 = output_buf[buf_offset + 3] + output_col; 693 /* VST4 of 16-bit elements completes the transpose. */ 694 vst4_lane_u16((uint16_t *)outptr0, cols_01_23_45_67, 0); 695 vst4_lane_u16((uint16_t *)outptr1, cols_01_23_45_67, 1); 696 vst4_lane_u16((uint16_t *)outptr2, cols_01_23_45_67, 2); 697 vst4_lane_u16((uint16_t *)outptr3, cols_01_23_45_67, 3); 698 } 699 700 701 /* Performs the second pass of the accurate inverse DCT on a 4x8 block 702 * of coefficients. 703 * 704 * This "sparse" version assumes that the coefficient values (after the first 705 * pass) in rows 4-7 are all 0. This simplifies the IDCT calculation, 706 * accelerating overall performance. 707 */ 708 709 static INLINE void jsimd_idct_islow_pass2_sparse(int16_t *workspace, 710 JSAMPARRAY output_buf, 711 JDIMENSION output_col, 712 unsigned buf_offset) 713 { 714 /* Load constants for IDCT computation. */ 715 #ifdef HAVE_VLD1_S16_X3 716 const int16x4x3_t consts = vld1_s16_x3(jsimd_idct_islow_neon_consts); 717 #else 718 const int16x4_t consts1 = vld1_s16(jsimd_idct_islow_neon_consts); 719 const int16x4_t consts2 = vld1_s16(jsimd_idct_islow_neon_consts + 4); 720 const int16x4_t consts3 = vld1_s16(jsimd_idct_islow_neon_consts + 8); 721 const int16x4x3_t consts = { { consts1, consts2, consts3 } }; 722 #endif 723 724 /* Even part (z3 is all 0) */ 725 int16x4_t z2_s16 = vld1_s16(workspace + 2 * DCTSIZE / 2); 726 727 int32x4_t tmp2 = vmull_lane_s16(z2_s16, consts.val[0], 1); 728 int32x4_t tmp3 = vmull_lane_s16(z2_s16, consts.val[1], 2); 729 730 z2_s16 = vld1_s16(workspace + 0 * DCTSIZE / 2); 731 int32x4_t tmp0 = vshll_n_s16(z2_s16, CONST_BITS); 732 int32x4_t tmp1 = vshll_n_s16(z2_s16, CONST_BITS); 733 734 int32x4_t tmp10 = vaddq_s32(tmp0, tmp3); 735 int32x4_t tmp13 = vsubq_s32(tmp0, tmp3); 736 int32x4_t tmp11 = vaddq_s32(tmp1, tmp2); 737 int32x4_t tmp12 = vsubq_s32(tmp1, tmp2); 738 739 /* Odd part (tmp0 and tmp1 are both all 0) */ 740 int16x4_t tmp2_s16 = vld1_s16(workspace + 3 * DCTSIZE / 2); 741 int16x4_t tmp3_s16 = vld1_s16(workspace + 1 * DCTSIZE / 2); 742 743 int16x4_t z3_s16 = tmp2_s16; 744 int16x4_t z4_s16 = tmp3_s16; 745 746 int32x4_t z3 = vmull_lane_s16(z3_s16, consts.val[2], 3); 747 z3 = vmlal_lane_s16(z3, z4_s16, consts.val[1], 3); 748 int32x4_t z4 = vmull_lane_s16(z3_s16, consts.val[1], 3); 749 z4 = vmlal_lane_s16(z4, z4_s16, consts.val[2], 0); 750 751 tmp0 = vmlsl_lane_s16(z3, tmp3_s16, consts.val[0], 0); 752 tmp1 = vmlsl_lane_s16(z4, tmp2_s16, consts.val[0], 2); 753 tmp2 = vmlal_lane_s16(z3, tmp2_s16, consts.val[2], 2); 754 tmp3 = vmlal_lane_s16(z4, tmp3_s16, consts.val[1], 0); 755 756 /* Final output stage: descale and narrow to 16-bit. */ 757 int16x8_t cols_02_s16 = vcombine_s16(vaddhn_s32(tmp10, tmp3), 758 vaddhn_s32(tmp12, tmp1)); 759 int16x8_t cols_13_s16 = vcombine_s16(vaddhn_s32(tmp11, tmp2), 760 vaddhn_s32(tmp13, tmp0)); 761 int16x8_t cols_46_s16 = vcombine_s16(vsubhn_s32(tmp13, tmp0), 762 vsubhn_s32(tmp11, tmp2)); 763 int16x8_t cols_57_s16 = vcombine_s16(vsubhn_s32(tmp12, tmp1), 764 vsubhn_s32(tmp10, tmp3)); 765 /* Descale and narrow to 8-bit. */ 766 int8x8_t cols_02_s8 = vqrshrn_n_s16(cols_02_s16, DESCALE_P2 - 16); 767 int8x8_t cols_13_s8 = vqrshrn_n_s16(cols_13_s16, DESCALE_P2 - 16); 768 int8x8_t cols_46_s8 = vqrshrn_n_s16(cols_46_s16, DESCALE_P2 - 16); 769 int8x8_t cols_57_s8 = vqrshrn_n_s16(cols_57_s16, DESCALE_P2 - 16); 770 /* Clamp to range [0-255]. */ 771 uint8x8_t cols_02_u8 = vadd_u8(vreinterpret_u8_s8(cols_02_s8), 772 vdup_n_u8(CENTERJSAMPLE)); 773 uint8x8_t cols_13_u8 = vadd_u8(vreinterpret_u8_s8(cols_13_s8), 774 vdup_n_u8(CENTERJSAMPLE)); 775 uint8x8_t cols_46_u8 = vadd_u8(vreinterpret_u8_s8(cols_46_s8), 776 vdup_n_u8(CENTERJSAMPLE)); 777 uint8x8_t cols_57_u8 = vadd_u8(vreinterpret_u8_s8(cols_57_s8), 778 vdup_n_u8(CENTERJSAMPLE)); 779 780 /* Transpose 4x8 block and store to memory. (Zipping adjacent columns 781 * together allows us to store 16-bit elements.) 782 */ 783 uint8x8x2_t cols_01_23 = vzip_u8(cols_02_u8, cols_13_u8); 784 uint8x8x2_t cols_45_67 = vzip_u8(cols_46_u8, cols_57_u8); 785 uint16x4x4_t cols_01_23_45_67 = { { 786 vreinterpret_u16_u8(cols_01_23.val[0]), 787 vreinterpret_u16_u8(cols_01_23.val[1]), 788 vreinterpret_u16_u8(cols_45_67.val[0]), 789 vreinterpret_u16_u8(cols_45_67.val[1]) 790 } }; 791 792 JSAMPROW outptr0 = output_buf[buf_offset + 0] + output_col; 793 JSAMPROW outptr1 = output_buf[buf_offset + 1] + output_col; 794 JSAMPROW outptr2 = output_buf[buf_offset + 2] + output_col; 795 JSAMPROW outptr3 = output_buf[buf_offset + 3] + output_col; 796 /* VST4 of 16-bit elements completes the transpose. */ 797 vst4_lane_u16((uint16_t *)outptr0, cols_01_23_45_67, 0); 798 vst4_lane_u16((uint16_t *)outptr1, cols_01_23_45_67, 1); 799 vst4_lane_u16((uint16_t *)outptr2, cols_01_23_45_67, 2); 800 vst4_lane_u16((uint16_t *)outptr3, cols_01_23_45_67, 3); 801 }