annotate ffmpeg/libavcodec/jrevdct.c @ 13:844d341cf643 tip

Back up before ISMIR
author Yading Song <yading.song@eecs.qmul.ac.uk>
date Thu, 31 Oct 2013 13:17:06 +0000
parents 6840f77b83aa
children
rev   line source
yading@10 1 /*
yading@10 2 * This file is part of the Independent JPEG Group's software.
yading@10 3 *
yading@10 4 * The authors make NO WARRANTY or representation, either express or implied,
yading@10 5 * with respect to this software, its quality, accuracy, merchantability, or
yading@10 6 * fitness for a particular purpose. This software is provided "AS IS", and
yading@10 7 * you, its user, assume the entire risk as to its quality and accuracy.
yading@10 8 *
yading@10 9 * This software is copyright (C) 1991, 1992, Thomas G. Lane.
yading@10 10 * All Rights Reserved except as specified below.
yading@10 11 *
yading@10 12 * Permission is hereby granted to use, copy, modify, and distribute this
yading@10 13 * software (or portions thereof) for any purpose, without fee, subject to
yading@10 14 * these conditions:
yading@10 15 * (1) If any part of the source code for this software is distributed, then
yading@10 16 * this README file must be included, with this copyright and no-warranty
yading@10 17 * notice unaltered; and any additions, deletions, or changes to the original
yading@10 18 * files must be clearly indicated in accompanying documentation.
yading@10 19 * (2) If only executable code is distributed, then the accompanying
yading@10 20 * documentation must state that "this software is based in part on the work
yading@10 21 * of the Independent JPEG Group".
yading@10 22 * (3) Permission for use of this software is granted only if the user accepts
yading@10 23 * full responsibility for any undesirable consequences; the authors accept
yading@10 24 * NO LIABILITY for damages of any kind.
yading@10 25 *
yading@10 26 * These conditions apply to any software derived from or based on the IJG
yading@10 27 * code, not just to the unmodified library. If you use our work, you ought
yading@10 28 * to acknowledge us.
yading@10 29 *
yading@10 30 * Permission is NOT granted for the use of any IJG author's name or company
yading@10 31 * name in advertising or publicity relating to this software or products
yading@10 32 * derived from it. This software may be referred to only as "the Independent
yading@10 33 * JPEG Group's software".
yading@10 34 *
yading@10 35 * We specifically permit and encourage the use of this software as the basis
yading@10 36 * of commercial products, provided that all warranty or liability claims are
yading@10 37 * assumed by the product vendor.
yading@10 38 *
yading@10 39 * This file contains the basic inverse-DCT transformation subroutine.
yading@10 40 *
yading@10 41 * This implementation is based on an algorithm described in
yading@10 42 * C. Loeffler, A. Ligtenberg and G. Moschytz, "Practical Fast 1-D DCT
yading@10 43 * Algorithms with 11 Multiplications", Proc. Int'l. Conf. on Acoustics,
yading@10 44 * Speech, and Signal Processing 1989 (ICASSP '89), pp. 988-991.
yading@10 45 * The primary algorithm described there uses 11 multiplies and 29 adds.
yading@10 46 * We use their alternate method with 12 multiplies and 32 adds.
yading@10 47 * The advantage of this method is that no data path contains more than one
yading@10 48 * multiplication; this allows a very simple and accurate implementation in
yading@10 49 * scaled fixed-point arithmetic, with a minimal number of shifts.
yading@10 50 *
yading@10 51 * I've made lots of modifications to attempt to take advantage of the
yading@10 52 * sparse nature of the DCT matrices we're getting. Although the logic
yading@10 53 * is cumbersome, it's straightforward and the resulting code is much
yading@10 54 * faster.
yading@10 55 *
yading@10 56 * A better way to do this would be to pass in the DCT block as a sparse
yading@10 57 * matrix, perhaps with the difference cases encoded.
yading@10 58 */
yading@10 59
yading@10 60 /**
yading@10 61 * @file
yading@10 62 * Independent JPEG Group's LLM idct.
yading@10 63 */
yading@10 64
yading@10 65 #include "libavutil/common.h"
yading@10 66 #include "dct.h"
yading@10 67
yading@10 68 #define EIGHT_BIT_SAMPLES
yading@10 69
yading@10 70 #define DCTSIZE 8
yading@10 71 #define DCTSIZE2 64
yading@10 72
yading@10 73 #define GLOBAL
yading@10 74
yading@10 75 #define RIGHT_SHIFT(x, n) ((x) >> (n))
yading@10 76
yading@10 77 typedef int16_t DCTBLOCK[DCTSIZE2];
yading@10 78
yading@10 79 #define CONST_BITS 13
yading@10 80
yading@10 81 /*
yading@10 82 * This routine is specialized to the case DCTSIZE = 8.
yading@10 83 */
yading@10 84
yading@10 85 #if DCTSIZE != 8
yading@10 86 Sorry, this code only copes with 8x8 DCTs. /* deliberate syntax err */
yading@10 87 #endif
yading@10 88
yading@10 89
yading@10 90 /*
yading@10 91 * A 2-D IDCT can be done by 1-D IDCT on each row followed by 1-D IDCT
yading@10 92 * on each column. Direct algorithms are also available, but they are
yading@10 93 * much more complex and seem not to be any faster when reduced to code.
yading@10 94 *
yading@10 95 * The poop on this scaling stuff is as follows:
yading@10 96 *
yading@10 97 * Each 1-D IDCT step produces outputs which are a factor of sqrt(N)
yading@10 98 * larger than the true IDCT outputs. The final outputs are therefore
yading@10 99 * a factor of N larger than desired; since N=8 this can be cured by
yading@10 100 * a simple right shift at the end of the algorithm. The advantage of
yading@10 101 * this arrangement is that we save two multiplications per 1-D IDCT,
yading@10 102 * because the y0 and y4 inputs need not be divided by sqrt(N).
yading@10 103 *
yading@10 104 * We have to do addition and subtraction of the integer inputs, which
yading@10 105 * is no problem, and multiplication by fractional constants, which is
yading@10 106 * a problem to do in integer arithmetic. We multiply all the constants
yading@10 107 * by CONST_SCALE and convert them to integer constants (thus retaining
yading@10 108 * CONST_BITS bits of precision in the constants). After doing a
yading@10 109 * multiplication we have to divide the product by CONST_SCALE, with proper
yading@10 110 * rounding, to produce the correct output. This division can be done
yading@10 111 * cheaply as a right shift of CONST_BITS bits. We postpone shifting
yading@10 112 * as long as possible so that partial sums can be added together with
yading@10 113 * full fractional precision.
yading@10 114 *
yading@10 115 * The outputs of the first pass are scaled up by PASS1_BITS bits so that
yading@10 116 * they are represented to better-than-integral precision. These outputs
yading@10 117 * require BITS_IN_JSAMPLE + PASS1_BITS + 3 bits; this fits in a 16-bit word
yading@10 118 * with the recommended scaling. (To scale up 12-bit sample data further, an
yading@10 119 * intermediate int32 array would be needed.)
yading@10 120 *
yading@10 121 * To avoid overflow of the 32-bit intermediate results in pass 2, we must
yading@10 122 * have BITS_IN_JSAMPLE + CONST_BITS + PASS1_BITS <= 26. Error analysis
yading@10 123 * shows that the values given below are the most effective.
yading@10 124 */
yading@10 125
yading@10 126 #ifdef EIGHT_BIT_SAMPLES
yading@10 127 #define PASS1_BITS 2
yading@10 128 #else
yading@10 129 #define PASS1_BITS 1 /* lose a little precision to avoid overflow */
yading@10 130 #endif
yading@10 131
yading@10 132 #define ONE ((int32_t) 1)
yading@10 133
yading@10 134 #define CONST_SCALE (ONE << CONST_BITS)
yading@10 135
yading@10 136 /* Convert a positive real constant to an integer scaled by CONST_SCALE.
yading@10 137 * IMPORTANT: if your compiler doesn't do this arithmetic at compile time,
yading@10 138 * you will pay a significant penalty in run time. In that case, figure
yading@10 139 * the correct integer constant values and insert them by hand.
yading@10 140 */
yading@10 141
yading@10 142 /* Actually FIX is no longer used, we precomputed them all */
yading@10 143 #define FIX(x) ((int32_t) ((x) * CONST_SCALE + 0.5))
yading@10 144
yading@10 145 /* Descale and correctly round an int32_t value that's scaled by N bits.
yading@10 146 * We assume RIGHT_SHIFT rounds towards minus infinity, so adding
yading@10 147 * the fudge factor is correct for either sign of X.
yading@10 148 */
yading@10 149
yading@10 150 #define DESCALE(x,n) RIGHT_SHIFT((x) + (ONE << ((n)-1)), n)
yading@10 151
yading@10 152 /* Multiply an int32_t variable by an int32_t constant to yield an int32_t result.
yading@10 153 * For 8-bit samples with the recommended scaling, all the variable
yading@10 154 * and constant values involved are no more than 16 bits wide, so a
yading@10 155 * 16x16->32 bit multiply can be used instead of a full 32x32 multiply;
yading@10 156 * this provides a useful speedup on many machines.
yading@10 157 * There is no way to specify a 16x16->32 multiply in portable C, but
yading@10 158 * some C compilers will do the right thing if you provide the correct
yading@10 159 * combination of casts.
yading@10 160 * NB: for 12-bit samples, a full 32-bit multiplication will be needed.
yading@10 161 */
yading@10 162
yading@10 163 #ifdef EIGHT_BIT_SAMPLES
yading@10 164 #ifdef SHORTxSHORT_32 /* may work if 'int' is 32 bits */
yading@10 165 #define MULTIPLY(var,const) (((int16_t) (var)) * ((int16_t) (const)))
yading@10 166 #endif
yading@10 167 #ifdef SHORTxLCONST_32 /* known to work with Microsoft C 6.0 */
yading@10 168 #define MULTIPLY(var,const) (((int16_t) (var)) * ((int32_t) (const)))
yading@10 169 #endif
yading@10 170 #endif
yading@10 171
yading@10 172 #ifndef MULTIPLY /* default definition */
yading@10 173 #define MULTIPLY(var,const) ((var) * (const))
yading@10 174 #endif
yading@10 175
yading@10 176
yading@10 177 /*
yading@10 178 Unlike our decoder where we approximate the FIXes, we need to use exact
yading@10 179 ones here or successive P-frames will drift too much with Reference frame coding
yading@10 180 */
yading@10 181 #define FIX_0_211164243 1730
yading@10 182 #define FIX_0_275899380 2260
yading@10 183 #define FIX_0_298631336 2446
yading@10 184 #define FIX_0_390180644 3196
yading@10 185 #define FIX_0_509795579 4176
yading@10 186 #define FIX_0_541196100 4433
yading@10 187 #define FIX_0_601344887 4926
yading@10 188 #define FIX_0_765366865 6270
yading@10 189 #define FIX_0_785694958 6436
yading@10 190 #define FIX_0_899976223 7373
yading@10 191 #define FIX_1_061594337 8697
yading@10 192 #define FIX_1_111140466 9102
yading@10 193 #define FIX_1_175875602 9633
yading@10 194 #define FIX_1_306562965 10703
yading@10 195 #define FIX_1_387039845 11363
yading@10 196 #define FIX_1_451774981 11893
yading@10 197 #define FIX_1_501321110 12299
yading@10 198 #define FIX_1_662939225 13623
yading@10 199 #define FIX_1_847759065 15137
yading@10 200 #define FIX_1_961570560 16069
yading@10 201 #define FIX_2_053119869 16819
yading@10 202 #define FIX_2_172734803 17799
yading@10 203 #define FIX_2_562915447 20995
yading@10 204 #define FIX_3_072711026 25172
yading@10 205
yading@10 206 /*
yading@10 207 * Perform the inverse DCT on one block of coefficients.
yading@10 208 */
yading@10 209
yading@10 210 void ff_j_rev_dct(DCTBLOCK data)
yading@10 211 {
yading@10 212 int32_t tmp0, tmp1, tmp2, tmp3;
yading@10 213 int32_t tmp10, tmp11, tmp12, tmp13;
yading@10 214 int32_t z1, z2, z3, z4, z5;
yading@10 215 int32_t d0, d1, d2, d3, d4, d5, d6, d7;
yading@10 216 register int16_t *dataptr;
yading@10 217 int rowctr;
yading@10 218
yading@10 219 /* Pass 1: process rows. */
yading@10 220 /* Note results are scaled up by sqrt(8) compared to a true IDCT; */
yading@10 221 /* furthermore, we scale the results by 2**PASS1_BITS. */
yading@10 222
yading@10 223 dataptr = data;
yading@10 224
yading@10 225 for (rowctr = DCTSIZE-1; rowctr >= 0; rowctr--) {
yading@10 226 /* Due to quantization, we will usually find that many of the input
yading@10 227 * coefficients are zero, especially the AC terms. We can exploit this
yading@10 228 * by short-circuiting the IDCT calculation for any row in which all
yading@10 229 * the AC terms are zero. In that case each output is equal to the
yading@10 230 * DC coefficient (with scale factor as needed).
yading@10 231 * With typical images and quantization tables, half or more of the
yading@10 232 * row DCT calculations can be simplified this way.
yading@10 233 */
yading@10 234
yading@10 235 register int *idataptr = (int*)dataptr;
yading@10 236
yading@10 237 /* WARNING: we do the same permutation as MMX idct to simplify the
yading@10 238 video core */
yading@10 239 d0 = dataptr[0];
yading@10 240 d2 = dataptr[1];
yading@10 241 d4 = dataptr[2];
yading@10 242 d6 = dataptr[3];
yading@10 243 d1 = dataptr[4];
yading@10 244 d3 = dataptr[5];
yading@10 245 d5 = dataptr[6];
yading@10 246 d7 = dataptr[7];
yading@10 247
yading@10 248 if ((d1 | d2 | d3 | d4 | d5 | d6 | d7) == 0) {
yading@10 249 /* AC terms all zero */
yading@10 250 if (d0) {
yading@10 251 /* Compute a 32 bit value to assign. */
yading@10 252 int16_t dcval = (int16_t) (d0 << PASS1_BITS);
yading@10 253 register int v = (dcval & 0xffff) | ((dcval << 16) & 0xffff0000);
yading@10 254
yading@10 255 idataptr[0] = v;
yading@10 256 idataptr[1] = v;
yading@10 257 idataptr[2] = v;
yading@10 258 idataptr[3] = v;
yading@10 259 }
yading@10 260
yading@10 261 dataptr += DCTSIZE; /* advance pointer to next row */
yading@10 262 continue;
yading@10 263 }
yading@10 264
yading@10 265 /* Even part: reverse the even part of the forward DCT. */
yading@10 266 /* The rotator is sqrt(2)*c(-6). */
yading@10 267 {
yading@10 268 if (d6) {
yading@10 269 if (d2) {
yading@10 270 /* d0 != 0, d2 != 0, d4 != 0, d6 != 0 */
yading@10 271 z1 = MULTIPLY(d2 + d6, FIX_0_541196100);
yading@10 272 tmp2 = z1 + MULTIPLY(-d6, FIX_1_847759065);
yading@10 273 tmp3 = z1 + MULTIPLY(d2, FIX_0_765366865);
yading@10 274
yading@10 275 tmp0 = (d0 + d4) << CONST_BITS;
yading@10 276 tmp1 = (d0 - d4) << CONST_BITS;
yading@10 277
yading@10 278 tmp10 = tmp0 + tmp3;
yading@10 279 tmp13 = tmp0 - tmp3;
yading@10 280 tmp11 = tmp1 + tmp2;
yading@10 281 tmp12 = tmp1 - tmp2;
yading@10 282 } else {
yading@10 283 /* d0 != 0, d2 == 0, d4 != 0, d6 != 0 */
yading@10 284 tmp2 = MULTIPLY(-d6, FIX_1_306562965);
yading@10 285 tmp3 = MULTIPLY(d6, FIX_0_541196100);
yading@10 286
yading@10 287 tmp0 = (d0 + d4) << CONST_BITS;
yading@10 288 tmp1 = (d0 - d4) << CONST_BITS;
yading@10 289
yading@10 290 tmp10 = tmp0 + tmp3;
yading@10 291 tmp13 = tmp0 - tmp3;
yading@10 292 tmp11 = tmp1 + tmp2;
yading@10 293 tmp12 = tmp1 - tmp2;
yading@10 294 }
yading@10 295 } else {
yading@10 296 if (d2) {
yading@10 297 /* d0 != 0, d2 != 0, d4 != 0, d6 == 0 */
yading@10 298 tmp2 = MULTIPLY(d2, FIX_0_541196100);
yading@10 299 tmp3 = MULTIPLY(d2, FIX_1_306562965);
yading@10 300
yading@10 301 tmp0 = (d0 + d4) << CONST_BITS;
yading@10 302 tmp1 = (d0 - d4) << CONST_BITS;
yading@10 303
yading@10 304 tmp10 = tmp0 + tmp3;
yading@10 305 tmp13 = tmp0 - tmp3;
yading@10 306 tmp11 = tmp1 + tmp2;
yading@10 307 tmp12 = tmp1 - tmp2;
yading@10 308 } else {
yading@10 309 /* d0 != 0, d2 == 0, d4 != 0, d6 == 0 */
yading@10 310 tmp10 = tmp13 = (d0 + d4) << CONST_BITS;
yading@10 311 tmp11 = tmp12 = (d0 - d4) << CONST_BITS;
yading@10 312 }
yading@10 313 }
yading@10 314
yading@10 315 /* Odd part per figure 8; the matrix is unitary and hence its
yading@10 316 * transpose is its inverse. i0..i3 are y7,y5,y3,y1 respectively.
yading@10 317 */
yading@10 318
yading@10 319 if (d7) {
yading@10 320 if (d5) {
yading@10 321 if (d3) {
yading@10 322 if (d1) {
yading@10 323 /* d1 != 0, d3 != 0, d5 != 0, d7 != 0 */
yading@10 324 z1 = d7 + d1;
yading@10 325 z2 = d5 + d3;
yading@10 326 z3 = d7 + d3;
yading@10 327 z4 = d5 + d1;
yading@10 328 z5 = MULTIPLY(z3 + z4, FIX_1_175875602);
yading@10 329
yading@10 330 tmp0 = MULTIPLY(d7, FIX_0_298631336);
yading@10 331 tmp1 = MULTIPLY(d5, FIX_2_053119869);
yading@10 332 tmp2 = MULTIPLY(d3, FIX_3_072711026);
yading@10 333 tmp3 = MULTIPLY(d1, FIX_1_501321110);
yading@10 334 z1 = MULTIPLY(-z1, FIX_0_899976223);
yading@10 335 z2 = MULTIPLY(-z2, FIX_2_562915447);
yading@10 336 z3 = MULTIPLY(-z3, FIX_1_961570560);
yading@10 337 z4 = MULTIPLY(-z4, FIX_0_390180644);
yading@10 338
yading@10 339 z3 += z5;
yading@10 340 z4 += z5;
yading@10 341
yading@10 342 tmp0 += z1 + z3;
yading@10 343 tmp1 += z2 + z4;
yading@10 344 tmp2 += z2 + z3;
yading@10 345 tmp3 += z1 + z4;
yading@10 346 } else {
yading@10 347 /* d1 == 0, d3 != 0, d5 != 0, d7 != 0 */
yading@10 348 z2 = d5 + d3;
yading@10 349 z3 = d7 + d3;
yading@10 350 z5 = MULTIPLY(z3 + d5, FIX_1_175875602);
yading@10 351
yading@10 352 tmp0 = MULTIPLY(d7, FIX_0_298631336);
yading@10 353 tmp1 = MULTIPLY(d5, FIX_2_053119869);
yading@10 354 tmp2 = MULTIPLY(d3, FIX_3_072711026);
yading@10 355 z1 = MULTIPLY(-d7, FIX_0_899976223);
yading@10 356 z2 = MULTIPLY(-z2, FIX_2_562915447);
yading@10 357 z3 = MULTIPLY(-z3, FIX_1_961570560);
yading@10 358 z4 = MULTIPLY(-d5, FIX_0_390180644);
yading@10 359
yading@10 360 z3 += z5;
yading@10 361 z4 += z5;
yading@10 362
yading@10 363 tmp0 += z1 + z3;
yading@10 364 tmp1 += z2 + z4;
yading@10 365 tmp2 += z2 + z3;
yading@10 366 tmp3 = z1 + z4;
yading@10 367 }
yading@10 368 } else {
yading@10 369 if (d1) {
yading@10 370 /* d1 != 0, d3 == 0, d5 != 0, d7 != 0 */
yading@10 371 z1 = d7 + d1;
yading@10 372 z4 = d5 + d1;
yading@10 373 z5 = MULTIPLY(d7 + z4, FIX_1_175875602);
yading@10 374
yading@10 375 tmp0 = MULTIPLY(d7, FIX_0_298631336);
yading@10 376 tmp1 = MULTIPLY(d5, FIX_2_053119869);
yading@10 377 tmp3 = MULTIPLY(d1, FIX_1_501321110);
yading@10 378 z1 = MULTIPLY(-z1, FIX_0_899976223);
yading@10 379 z2 = MULTIPLY(-d5, FIX_2_562915447);
yading@10 380 z3 = MULTIPLY(-d7, FIX_1_961570560);
yading@10 381 z4 = MULTIPLY(-z4, FIX_0_390180644);
yading@10 382
yading@10 383 z3 += z5;
yading@10 384 z4 += z5;
yading@10 385
yading@10 386 tmp0 += z1 + z3;
yading@10 387 tmp1 += z2 + z4;
yading@10 388 tmp2 = z2 + z3;
yading@10 389 tmp3 += z1 + z4;
yading@10 390 } else {
yading@10 391 /* d1 == 0, d3 == 0, d5 != 0, d7 != 0 */
yading@10 392 tmp0 = MULTIPLY(-d7, FIX_0_601344887);
yading@10 393 z1 = MULTIPLY(-d7, FIX_0_899976223);
yading@10 394 z3 = MULTIPLY(-d7, FIX_1_961570560);
yading@10 395 tmp1 = MULTIPLY(-d5, FIX_0_509795579);
yading@10 396 z2 = MULTIPLY(-d5, FIX_2_562915447);
yading@10 397 z4 = MULTIPLY(-d5, FIX_0_390180644);
yading@10 398 z5 = MULTIPLY(d5 + d7, FIX_1_175875602);
yading@10 399
yading@10 400 z3 += z5;
yading@10 401 z4 += z5;
yading@10 402
yading@10 403 tmp0 += z3;
yading@10 404 tmp1 += z4;
yading@10 405 tmp2 = z2 + z3;
yading@10 406 tmp3 = z1 + z4;
yading@10 407 }
yading@10 408 }
yading@10 409 } else {
yading@10 410 if (d3) {
yading@10 411 if (d1) {
yading@10 412 /* d1 != 0, d3 != 0, d5 == 0, d7 != 0 */
yading@10 413 z1 = d7 + d1;
yading@10 414 z3 = d7 + d3;
yading@10 415 z5 = MULTIPLY(z3 + d1, FIX_1_175875602);
yading@10 416
yading@10 417 tmp0 = MULTIPLY(d7, FIX_0_298631336);
yading@10 418 tmp2 = MULTIPLY(d3, FIX_3_072711026);
yading@10 419 tmp3 = MULTIPLY(d1, FIX_1_501321110);
yading@10 420 z1 = MULTIPLY(-z1, FIX_0_899976223);
yading@10 421 z2 = MULTIPLY(-d3, FIX_2_562915447);
yading@10 422 z3 = MULTIPLY(-z3, FIX_1_961570560);
yading@10 423 z4 = MULTIPLY(-d1, FIX_0_390180644);
yading@10 424
yading@10 425 z3 += z5;
yading@10 426 z4 += z5;
yading@10 427
yading@10 428 tmp0 += z1 + z3;
yading@10 429 tmp1 = z2 + z4;
yading@10 430 tmp2 += z2 + z3;
yading@10 431 tmp3 += z1 + z4;
yading@10 432 } else {
yading@10 433 /* d1 == 0, d3 != 0, d5 == 0, d7 != 0 */
yading@10 434 z3 = d7 + d3;
yading@10 435
yading@10 436 tmp0 = MULTIPLY(-d7, FIX_0_601344887);
yading@10 437 z1 = MULTIPLY(-d7, FIX_0_899976223);
yading@10 438 tmp2 = MULTIPLY(d3, FIX_0_509795579);
yading@10 439 z2 = MULTIPLY(-d3, FIX_2_562915447);
yading@10 440 z5 = MULTIPLY(z3, FIX_1_175875602);
yading@10 441 z3 = MULTIPLY(-z3, FIX_0_785694958);
yading@10 442
yading@10 443 tmp0 += z3;
yading@10 444 tmp1 = z2 + z5;
yading@10 445 tmp2 += z3;
yading@10 446 tmp3 = z1 + z5;
yading@10 447 }
yading@10 448 } else {
yading@10 449 if (d1) {
yading@10 450 /* d1 != 0, d3 == 0, d5 == 0, d7 != 0 */
yading@10 451 z1 = d7 + d1;
yading@10 452 z5 = MULTIPLY(z1, FIX_1_175875602);
yading@10 453
yading@10 454 z1 = MULTIPLY(z1, FIX_0_275899380);
yading@10 455 z3 = MULTIPLY(-d7, FIX_1_961570560);
yading@10 456 tmp0 = MULTIPLY(-d7, FIX_1_662939225);
yading@10 457 z4 = MULTIPLY(-d1, FIX_0_390180644);
yading@10 458 tmp3 = MULTIPLY(d1, FIX_1_111140466);
yading@10 459
yading@10 460 tmp0 += z1;
yading@10 461 tmp1 = z4 + z5;
yading@10 462 tmp2 = z3 + z5;
yading@10 463 tmp3 += z1;
yading@10 464 } else {
yading@10 465 /* d1 == 0, d3 == 0, d5 == 0, d7 != 0 */
yading@10 466 tmp0 = MULTIPLY(-d7, FIX_1_387039845);
yading@10 467 tmp1 = MULTIPLY(d7, FIX_1_175875602);
yading@10 468 tmp2 = MULTIPLY(-d7, FIX_0_785694958);
yading@10 469 tmp3 = MULTIPLY(d7, FIX_0_275899380);
yading@10 470 }
yading@10 471 }
yading@10 472 }
yading@10 473 } else {
yading@10 474 if (d5) {
yading@10 475 if (d3) {
yading@10 476 if (d1) {
yading@10 477 /* d1 != 0, d3 != 0, d5 != 0, d7 == 0 */
yading@10 478 z2 = d5 + d3;
yading@10 479 z4 = d5 + d1;
yading@10 480 z5 = MULTIPLY(d3 + z4, FIX_1_175875602);
yading@10 481
yading@10 482 tmp1 = MULTIPLY(d5, FIX_2_053119869);
yading@10 483 tmp2 = MULTIPLY(d3, FIX_3_072711026);
yading@10 484 tmp3 = MULTIPLY(d1, FIX_1_501321110);
yading@10 485 z1 = MULTIPLY(-d1, FIX_0_899976223);
yading@10 486 z2 = MULTIPLY(-z2, FIX_2_562915447);
yading@10 487 z3 = MULTIPLY(-d3, FIX_1_961570560);
yading@10 488 z4 = MULTIPLY(-z4, FIX_0_390180644);
yading@10 489
yading@10 490 z3 += z5;
yading@10 491 z4 += z5;
yading@10 492
yading@10 493 tmp0 = z1 + z3;
yading@10 494 tmp1 += z2 + z4;
yading@10 495 tmp2 += z2 + z3;
yading@10 496 tmp3 += z1 + z4;
yading@10 497 } else {
yading@10 498 /* d1 == 0, d3 != 0, d5 != 0, d7 == 0 */
yading@10 499 z2 = d5 + d3;
yading@10 500
yading@10 501 z5 = MULTIPLY(z2, FIX_1_175875602);
yading@10 502 tmp1 = MULTIPLY(d5, FIX_1_662939225);
yading@10 503 z4 = MULTIPLY(-d5, FIX_0_390180644);
yading@10 504 z2 = MULTIPLY(-z2, FIX_1_387039845);
yading@10 505 tmp2 = MULTIPLY(d3, FIX_1_111140466);
yading@10 506 z3 = MULTIPLY(-d3, FIX_1_961570560);
yading@10 507
yading@10 508 tmp0 = z3 + z5;
yading@10 509 tmp1 += z2;
yading@10 510 tmp2 += z2;
yading@10 511 tmp3 = z4 + z5;
yading@10 512 }
yading@10 513 } else {
yading@10 514 if (d1) {
yading@10 515 /* d1 != 0, d3 == 0, d5 != 0, d7 == 0 */
yading@10 516 z4 = d5 + d1;
yading@10 517
yading@10 518 z5 = MULTIPLY(z4, FIX_1_175875602);
yading@10 519 z1 = MULTIPLY(-d1, FIX_0_899976223);
yading@10 520 tmp3 = MULTIPLY(d1, FIX_0_601344887);
yading@10 521 tmp1 = MULTIPLY(-d5, FIX_0_509795579);
yading@10 522 z2 = MULTIPLY(-d5, FIX_2_562915447);
yading@10 523 z4 = MULTIPLY(z4, FIX_0_785694958);
yading@10 524
yading@10 525 tmp0 = z1 + z5;
yading@10 526 tmp1 += z4;
yading@10 527 tmp2 = z2 + z5;
yading@10 528 tmp3 += z4;
yading@10 529 } else {
yading@10 530 /* d1 == 0, d3 == 0, d5 != 0, d7 == 0 */
yading@10 531 tmp0 = MULTIPLY(d5, FIX_1_175875602);
yading@10 532 tmp1 = MULTIPLY(d5, FIX_0_275899380);
yading@10 533 tmp2 = MULTIPLY(-d5, FIX_1_387039845);
yading@10 534 tmp3 = MULTIPLY(d5, FIX_0_785694958);
yading@10 535 }
yading@10 536 }
yading@10 537 } else {
yading@10 538 if (d3) {
yading@10 539 if (d1) {
yading@10 540 /* d1 != 0, d3 != 0, d5 == 0, d7 == 0 */
yading@10 541 z5 = d1 + d3;
yading@10 542 tmp3 = MULTIPLY(d1, FIX_0_211164243);
yading@10 543 tmp2 = MULTIPLY(-d3, FIX_1_451774981);
yading@10 544 z1 = MULTIPLY(d1, FIX_1_061594337);
yading@10 545 z2 = MULTIPLY(-d3, FIX_2_172734803);
yading@10 546 z4 = MULTIPLY(z5, FIX_0_785694958);
yading@10 547 z5 = MULTIPLY(z5, FIX_1_175875602);
yading@10 548
yading@10 549 tmp0 = z1 - z4;
yading@10 550 tmp1 = z2 + z4;
yading@10 551 tmp2 += z5;
yading@10 552 tmp3 += z5;
yading@10 553 } else {
yading@10 554 /* d1 == 0, d3 != 0, d5 == 0, d7 == 0 */
yading@10 555 tmp0 = MULTIPLY(-d3, FIX_0_785694958);
yading@10 556 tmp1 = MULTIPLY(-d3, FIX_1_387039845);
yading@10 557 tmp2 = MULTIPLY(-d3, FIX_0_275899380);
yading@10 558 tmp3 = MULTIPLY(d3, FIX_1_175875602);
yading@10 559 }
yading@10 560 } else {
yading@10 561 if (d1) {
yading@10 562 /* d1 != 0, d3 == 0, d5 == 0, d7 == 0 */
yading@10 563 tmp0 = MULTIPLY(d1, FIX_0_275899380);
yading@10 564 tmp1 = MULTIPLY(d1, FIX_0_785694958);
yading@10 565 tmp2 = MULTIPLY(d1, FIX_1_175875602);
yading@10 566 tmp3 = MULTIPLY(d1, FIX_1_387039845);
yading@10 567 } else {
yading@10 568 /* d1 == 0, d3 == 0, d5 == 0, d7 == 0 */
yading@10 569 tmp0 = tmp1 = tmp2 = tmp3 = 0;
yading@10 570 }
yading@10 571 }
yading@10 572 }
yading@10 573 }
yading@10 574 }
yading@10 575 /* Final output stage: inputs are tmp10..tmp13, tmp0..tmp3 */
yading@10 576
yading@10 577 dataptr[0] = (int16_t) DESCALE(tmp10 + tmp3, CONST_BITS-PASS1_BITS);
yading@10 578 dataptr[7] = (int16_t) DESCALE(tmp10 - tmp3, CONST_BITS-PASS1_BITS);
yading@10 579 dataptr[1] = (int16_t) DESCALE(tmp11 + tmp2, CONST_BITS-PASS1_BITS);
yading@10 580 dataptr[6] = (int16_t) DESCALE(tmp11 - tmp2, CONST_BITS-PASS1_BITS);
yading@10 581 dataptr[2] = (int16_t) DESCALE(tmp12 + tmp1, CONST_BITS-PASS1_BITS);
yading@10 582 dataptr[5] = (int16_t) DESCALE(tmp12 - tmp1, CONST_BITS-PASS1_BITS);
yading@10 583 dataptr[3] = (int16_t) DESCALE(tmp13 + tmp0, CONST_BITS-PASS1_BITS);
yading@10 584 dataptr[4] = (int16_t) DESCALE(tmp13 - tmp0, CONST_BITS-PASS1_BITS);
yading@10 585
yading@10 586 dataptr += DCTSIZE; /* advance pointer to next row */
yading@10 587 }
yading@10 588
yading@10 589 /* Pass 2: process columns. */
yading@10 590 /* Note that we must descale the results by a factor of 8 == 2**3, */
yading@10 591 /* and also undo the PASS1_BITS scaling. */
yading@10 592
yading@10 593 dataptr = data;
yading@10 594 for (rowctr = DCTSIZE-1; rowctr >= 0; rowctr--) {
yading@10 595 /* Columns of zeroes can be exploited in the same way as we did with rows.
yading@10 596 * However, the row calculation has created many nonzero AC terms, so the
yading@10 597 * simplification applies less often (typically 5% to 10% of the time).
yading@10 598 * On machines with very fast multiplication, it's possible that the
yading@10 599 * test takes more time than it's worth. In that case this section
yading@10 600 * may be commented out.
yading@10 601 */
yading@10 602
yading@10 603 d0 = dataptr[DCTSIZE*0];
yading@10 604 d1 = dataptr[DCTSIZE*1];
yading@10 605 d2 = dataptr[DCTSIZE*2];
yading@10 606 d3 = dataptr[DCTSIZE*3];
yading@10 607 d4 = dataptr[DCTSIZE*4];
yading@10 608 d5 = dataptr[DCTSIZE*5];
yading@10 609 d6 = dataptr[DCTSIZE*6];
yading@10 610 d7 = dataptr[DCTSIZE*7];
yading@10 611
yading@10 612 /* Even part: reverse the even part of the forward DCT. */
yading@10 613 /* The rotator is sqrt(2)*c(-6). */
yading@10 614 if (d6) {
yading@10 615 if (d2) {
yading@10 616 /* d0 != 0, d2 != 0, d4 != 0, d6 != 0 */
yading@10 617 z1 = MULTIPLY(d2 + d6, FIX_0_541196100);
yading@10 618 tmp2 = z1 + MULTIPLY(-d6, FIX_1_847759065);
yading@10 619 tmp3 = z1 + MULTIPLY(d2, FIX_0_765366865);
yading@10 620
yading@10 621 tmp0 = (d0 + d4) << CONST_BITS;
yading@10 622 tmp1 = (d0 - d4) << CONST_BITS;
yading@10 623
yading@10 624 tmp10 = tmp0 + tmp3;
yading@10 625 tmp13 = tmp0 - tmp3;
yading@10 626 tmp11 = tmp1 + tmp2;
yading@10 627 tmp12 = tmp1 - tmp2;
yading@10 628 } else {
yading@10 629 /* d0 != 0, d2 == 0, d4 != 0, d6 != 0 */
yading@10 630 tmp2 = MULTIPLY(-d6, FIX_1_306562965);
yading@10 631 tmp3 = MULTIPLY(d6, FIX_0_541196100);
yading@10 632
yading@10 633 tmp0 = (d0 + d4) << CONST_BITS;
yading@10 634 tmp1 = (d0 - d4) << CONST_BITS;
yading@10 635
yading@10 636 tmp10 = tmp0 + tmp3;
yading@10 637 tmp13 = tmp0 - tmp3;
yading@10 638 tmp11 = tmp1 + tmp2;
yading@10 639 tmp12 = tmp1 - tmp2;
yading@10 640 }
yading@10 641 } else {
yading@10 642 if (d2) {
yading@10 643 /* d0 != 0, d2 != 0, d4 != 0, d6 == 0 */
yading@10 644 tmp2 = MULTIPLY(d2, FIX_0_541196100);
yading@10 645 tmp3 = MULTIPLY(d2, FIX_1_306562965);
yading@10 646
yading@10 647 tmp0 = (d0 + d4) << CONST_BITS;
yading@10 648 tmp1 = (d0 - d4) << CONST_BITS;
yading@10 649
yading@10 650 tmp10 = tmp0 + tmp3;
yading@10 651 tmp13 = tmp0 - tmp3;
yading@10 652 tmp11 = tmp1 + tmp2;
yading@10 653 tmp12 = tmp1 - tmp2;
yading@10 654 } else {
yading@10 655 /* d0 != 0, d2 == 0, d4 != 0, d6 == 0 */
yading@10 656 tmp10 = tmp13 = (d0 + d4) << CONST_BITS;
yading@10 657 tmp11 = tmp12 = (d0 - d4) << CONST_BITS;
yading@10 658 }
yading@10 659 }
yading@10 660
yading@10 661 /* Odd part per figure 8; the matrix is unitary and hence its
yading@10 662 * transpose is its inverse. i0..i3 are y7,y5,y3,y1 respectively.
yading@10 663 */
yading@10 664 if (d7) {
yading@10 665 if (d5) {
yading@10 666 if (d3) {
yading@10 667 if (d1) {
yading@10 668 /* d1 != 0, d3 != 0, d5 != 0, d7 != 0 */
yading@10 669 z1 = d7 + d1;
yading@10 670 z2 = d5 + d3;
yading@10 671 z3 = d7 + d3;
yading@10 672 z4 = d5 + d1;
yading@10 673 z5 = MULTIPLY(z3 + z4, FIX_1_175875602);
yading@10 674
yading@10 675 tmp0 = MULTIPLY(d7, FIX_0_298631336);
yading@10 676 tmp1 = MULTIPLY(d5, FIX_2_053119869);
yading@10 677 tmp2 = MULTIPLY(d3, FIX_3_072711026);
yading@10 678 tmp3 = MULTIPLY(d1, FIX_1_501321110);
yading@10 679 z1 = MULTIPLY(-z1, FIX_0_899976223);
yading@10 680 z2 = MULTIPLY(-z2, FIX_2_562915447);
yading@10 681 z3 = MULTIPLY(-z3, FIX_1_961570560);
yading@10 682 z4 = MULTIPLY(-z4, FIX_0_390180644);
yading@10 683
yading@10 684 z3 += z5;
yading@10 685 z4 += z5;
yading@10 686
yading@10 687 tmp0 += z1 + z3;
yading@10 688 tmp1 += z2 + z4;
yading@10 689 tmp2 += z2 + z3;
yading@10 690 tmp3 += z1 + z4;
yading@10 691 } else {
yading@10 692 /* d1 == 0, d3 != 0, d5 != 0, d7 != 0 */
yading@10 693 z2 = d5 + d3;
yading@10 694 z3 = d7 + d3;
yading@10 695 z5 = MULTIPLY(z3 + d5, FIX_1_175875602);
yading@10 696
yading@10 697 tmp0 = MULTIPLY(d7, FIX_0_298631336);
yading@10 698 tmp1 = MULTIPLY(d5, FIX_2_053119869);
yading@10 699 tmp2 = MULTIPLY(d3, FIX_3_072711026);
yading@10 700 z1 = MULTIPLY(-d7, FIX_0_899976223);
yading@10 701 z2 = MULTIPLY(-z2, FIX_2_562915447);
yading@10 702 z3 = MULTIPLY(-z3, FIX_1_961570560);
yading@10 703 z4 = MULTIPLY(-d5, FIX_0_390180644);
yading@10 704
yading@10 705 z3 += z5;
yading@10 706 z4 += z5;
yading@10 707
yading@10 708 tmp0 += z1 + z3;
yading@10 709 tmp1 += z2 + z4;
yading@10 710 tmp2 += z2 + z3;
yading@10 711 tmp3 = z1 + z4;
yading@10 712 }
yading@10 713 } else {
yading@10 714 if (d1) {
yading@10 715 /* d1 != 0, d3 == 0, d5 != 0, d7 != 0 */
yading@10 716 z1 = d7 + d1;
yading@10 717 z3 = d7;
yading@10 718 z4 = d5 + d1;
yading@10 719 z5 = MULTIPLY(z3 + z4, FIX_1_175875602);
yading@10 720
yading@10 721 tmp0 = MULTIPLY(d7, FIX_0_298631336);
yading@10 722 tmp1 = MULTIPLY(d5, FIX_2_053119869);
yading@10 723 tmp3 = MULTIPLY(d1, FIX_1_501321110);
yading@10 724 z1 = MULTIPLY(-z1, FIX_0_899976223);
yading@10 725 z2 = MULTIPLY(-d5, FIX_2_562915447);
yading@10 726 z3 = MULTIPLY(-d7, FIX_1_961570560);
yading@10 727 z4 = MULTIPLY(-z4, FIX_0_390180644);
yading@10 728
yading@10 729 z3 += z5;
yading@10 730 z4 += z5;
yading@10 731
yading@10 732 tmp0 += z1 + z3;
yading@10 733 tmp1 += z2 + z4;
yading@10 734 tmp2 = z2 + z3;
yading@10 735 tmp3 += z1 + z4;
yading@10 736 } else {
yading@10 737 /* d1 == 0, d3 == 0, d5 != 0, d7 != 0 */
yading@10 738 tmp0 = MULTIPLY(-d7, FIX_0_601344887);
yading@10 739 z1 = MULTIPLY(-d7, FIX_0_899976223);
yading@10 740 z3 = MULTIPLY(-d7, FIX_1_961570560);
yading@10 741 tmp1 = MULTIPLY(-d5, FIX_0_509795579);
yading@10 742 z2 = MULTIPLY(-d5, FIX_2_562915447);
yading@10 743 z4 = MULTIPLY(-d5, FIX_0_390180644);
yading@10 744 z5 = MULTIPLY(d5 + d7, FIX_1_175875602);
yading@10 745
yading@10 746 z3 += z5;
yading@10 747 z4 += z5;
yading@10 748
yading@10 749 tmp0 += z3;
yading@10 750 tmp1 += z4;
yading@10 751 tmp2 = z2 + z3;
yading@10 752 tmp3 = z1 + z4;
yading@10 753 }
yading@10 754 }
yading@10 755 } else {
yading@10 756 if (d3) {
yading@10 757 if (d1) {
yading@10 758 /* d1 != 0, d3 != 0, d5 == 0, d7 != 0 */
yading@10 759 z1 = d7 + d1;
yading@10 760 z3 = d7 + d3;
yading@10 761 z5 = MULTIPLY(z3 + d1, FIX_1_175875602);
yading@10 762
yading@10 763 tmp0 = MULTIPLY(d7, FIX_0_298631336);
yading@10 764 tmp2 = MULTIPLY(d3, FIX_3_072711026);
yading@10 765 tmp3 = MULTIPLY(d1, FIX_1_501321110);
yading@10 766 z1 = MULTIPLY(-z1, FIX_0_899976223);
yading@10 767 z2 = MULTIPLY(-d3, FIX_2_562915447);
yading@10 768 z3 = MULTIPLY(-z3, FIX_1_961570560);
yading@10 769 z4 = MULTIPLY(-d1, FIX_0_390180644);
yading@10 770
yading@10 771 z3 += z5;
yading@10 772 z4 += z5;
yading@10 773
yading@10 774 tmp0 += z1 + z3;
yading@10 775 tmp1 = z2 + z4;
yading@10 776 tmp2 += z2 + z3;
yading@10 777 tmp3 += z1 + z4;
yading@10 778 } else {
yading@10 779 /* d1 == 0, d3 != 0, d5 == 0, d7 != 0 */
yading@10 780 z3 = d7 + d3;
yading@10 781
yading@10 782 tmp0 = MULTIPLY(-d7, FIX_0_601344887);
yading@10 783 z1 = MULTIPLY(-d7, FIX_0_899976223);
yading@10 784 tmp2 = MULTIPLY(d3, FIX_0_509795579);
yading@10 785 z2 = MULTIPLY(-d3, FIX_2_562915447);
yading@10 786 z5 = MULTIPLY(z3, FIX_1_175875602);
yading@10 787 z3 = MULTIPLY(-z3, FIX_0_785694958);
yading@10 788
yading@10 789 tmp0 += z3;
yading@10 790 tmp1 = z2 + z5;
yading@10 791 tmp2 += z3;
yading@10 792 tmp3 = z1 + z5;
yading@10 793 }
yading@10 794 } else {
yading@10 795 if (d1) {
yading@10 796 /* d1 != 0, d3 == 0, d5 == 0, d7 != 0 */
yading@10 797 z1 = d7 + d1;
yading@10 798 z5 = MULTIPLY(z1, FIX_1_175875602);
yading@10 799
yading@10 800 z1 = MULTIPLY(z1, FIX_0_275899380);
yading@10 801 z3 = MULTIPLY(-d7, FIX_1_961570560);
yading@10 802 tmp0 = MULTIPLY(-d7, FIX_1_662939225);
yading@10 803 z4 = MULTIPLY(-d1, FIX_0_390180644);
yading@10 804 tmp3 = MULTIPLY(d1, FIX_1_111140466);
yading@10 805
yading@10 806 tmp0 += z1;
yading@10 807 tmp1 = z4 + z5;
yading@10 808 tmp2 = z3 + z5;
yading@10 809 tmp3 += z1;
yading@10 810 } else {
yading@10 811 /* d1 == 0, d3 == 0, d5 == 0, d7 != 0 */
yading@10 812 tmp0 = MULTIPLY(-d7, FIX_1_387039845);
yading@10 813 tmp1 = MULTIPLY(d7, FIX_1_175875602);
yading@10 814 tmp2 = MULTIPLY(-d7, FIX_0_785694958);
yading@10 815 tmp3 = MULTIPLY(d7, FIX_0_275899380);
yading@10 816 }
yading@10 817 }
yading@10 818 }
yading@10 819 } else {
yading@10 820 if (d5) {
yading@10 821 if (d3) {
yading@10 822 if (d1) {
yading@10 823 /* d1 != 0, d3 != 0, d5 != 0, d7 == 0 */
yading@10 824 z2 = d5 + d3;
yading@10 825 z4 = d5 + d1;
yading@10 826 z5 = MULTIPLY(d3 + z4, FIX_1_175875602);
yading@10 827
yading@10 828 tmp1 = MULTIPLY(d5, FIX_2_053119869);
yading@10 829 tmp2 = MULTIPLY(d3, FIX_3_072711026);
yading@10 830 tmp3 = MULTIPLY(d1, FIX_1_501321110);
yading@10 831 z1 = MULTIPLY(-d1, FIX_0_899976223);
yading@10 832 z2 = MULTIPLY(-z2, FIX_2_562915447);
yading@10 833 z3 = MULTIPLY(-d3, FIX_1_961570560);
yading@10 834 z4 = MULTIPLY(-z4, FIX_0_390180644);
yading@10 835
yading@10 836 z3 += z5;
yading@10 837 z4 += z5;
yading@10 838
yading@10 839 tmp0 = z1 + z3;
yading@10 840 tmp1 += z2 + z4;
yading@10 841 tmp2 += z2 + z3;
yading@10 842 tmp3 += z1 + z4;
yading@10 843 } else {
yading@10 844 /* d1 == 0, d3 != 0, d5 != 0, d7 == 0 */
yading@10 845 z2 = d5 + d3;
yading@10 846
yading@10 847 z5 = MULTIPLY(z2, FIX_1_175875602);
yading@10 848 tmp1 = MULTIPLY(d5, FIX_1_662939225);
yading@10 849 z4 = MULTIPLY(-d5, FIX_0_390180644);
yading@10 850 z2 = MULTIPLY(-z2, FIX_1_387039845);
yading@10 851 tmp2 = MULTIPLY(d3, FIX_1_111140466);
yading@10 852 z3 = MULTIPLY(-d3, FIX_1_961570560);
yading@10 853
yading@10 854 tmp0 = z3 + z5;
yading@10 855 tmp1 += z2;
yading@10 856 tmp2 += z2;
yading@10 857 tmp3 = z4 + z5;
yading@10 858 }
yading@10 859 } else {
yading@10 860 if (d1) {
yading@10 861 /* d1 != 0, d3 == 0, d5 != 0, d7 == 0 */
yading@10 862 z4 = d5 + d1;
yading@10 863
yading@10 864 z5 = MULTIPLY(z4, FIX_1_175875602);
yading@10 865 z1 = MULTIPLY(-d1, FIX_0_899976223);
yading@10 866 tmp3 = MULTIPLY(d1, FIX_0_601344887);
yading@10 867 tmp1 = MULTIPLY(-d5, FIX_0_509795579);
yading@10 868 z2 = MULTIPLY(-d5, FIX_2_562915447);
yading@10 869 z4 = MULTIPLY(z4, FIX_0_785694958);
yading@10 870
yading@10 871 tmp0 = z1 + z5;
yading@10 872 tmp1 += z4;
yading@10 873 tmp2 = z2 + z5;
yading@10 874 tmp3 += z4;
yading@10 875 } else {
yading@10 876 /* d1 == 0, d3 == 0, d5 != 0, d7 == 0 */
yading@10 877 tmp0 = MULTIPLY(d5, FIX_1_175875602);
yading@10 878 tmp1 = MULTIPLY(d5, FIX_0_275899380);
yading@10 879 tmp2 = MULTIPLY(-d5, FIX_1_387039845);
yading@10 880 tmp3 = MULTIPLY(d5, FIX_0_785694958);
yading@10 881 }
yading@10 882 }
yading@10 883 } else {
yading@10 884 if (d3) {
yading@10 885 if (d1) {
yading@10 886 /* d1 != 0, d3 != 0, d5 == 0, d7 == 0 */
yading@10 887 z5 = d1 + d3;
yading@10 888 tmp3 = MULTIPLY(d1, FIX_0_211164243);
yading@10 889 tmp2 = MULTIPLY(-d3, FIX_1_451774981);
yading@10 890 z1 = MULTIPLY(d1, FIX_1_061594337);
yading@10 891 z2 = MULTIPLY(-d3, FIX_2_172734803);
yading@10 892 z4 = MULTIPLY(z5, FIX_0_785694958);
yading@10 893 z5 = MULTIPLY(z5, FIX_1_175875602);
yading@10 894
yading@10 895 tmp0 = z1 - z4;
yading@10 896 tmp1 = z2 + z4;
yading@10 897 tmp2 += z5;
yading@10 898 tmp3 += z5;
yading@10 899 } else {
yading@10 900 /* d1 == 0, d3 != 0, d5 == 0, d7 == 0 */
yading@10 901 tmp0 = MULTIPLY(-d3, FIX_0_785694958);
yading@10 902 tmp1 = MULTIPLY(-d3, FIX_1_387039845);
yading@10 903 tmp2 = MULTIPLY(-d3, FIX_0_275899380);
yading@10 904 tmp3 = MULTIPLY(d3, FIX_1_175875602);
yading@10 905 }
yading@10 906 } else {
yading@10 907 if (d1) {
yading@10 908 /* d1 != 0, d3 == 0, d5 == 0, d7 == 0 */
yading@10 909 tmp0 = MULTIPLY(d1, FIX_0_275899380);
yading@10 910 tmp1 = MULTIPLY(d1, FIX_0_785694958);
yading@10 911 tmp2 = MULTIPLY(d1, FIX_1_175875602);
yading@10 912 tmp3 = MULTIPLY(d1, FIX_1_387039845);
yading@10 913 } else {
yading@10 914 /* d1 == 0, d3 == 0, d5 == 0, d7 == 0 */
yading@10 915 tmp0 = tmp1 = tmp2 = tmp3 = 0;
yading@10 916 }
yading@10 917 }
yading@10 918 }
yading@10 919 }
yading@10 920
yading@10 921 /* Final output stage: inputs are tmp10..tmp13, tmp0..tmp3 */
yading@10 922
yading@10 923 dataptr[DCTSIZE*0] = (int16_t) DESCALE(tmp10 + tmp3,
yading@10 924 CONST_BITS+PASS1_BITS+3);
yading@10 925 dataptr[DCTSIZE*7] = (int16_t) DESCALE(tmp10 - tmp3,
yading@10 926 CONST_BITS+PASS1_BITS+3);
yading@10 927 dataptr[DCTSIZE*1] = (int16_t) DESCALE(tmp11 + tmp2,
yading@10 928 CONST_BITS+PASS1_BITS+3);
yading@10 929 dataptr[DCTSIZE*6] = (int16_t) DESCALE(tmp11 - tmp2,
yading@10 930 CONST_BITS+PASS1_BITS+3);
yading@10 931 dataptr[DCTSIZE*2] = (int16_t) DESCALE(tmp12 + tmp1,
yading@10 932 CONST_BITS+PASS1_BITS+3);
yading@10 933 dataptr[DCTSIZE*5] = (int16_t) DESCALE(tmp12 - tmp1,
yading@10 934 CONST_BITS+PASS1_BITS+3);
yading@10 935 dataptr[DCTSIZE*3] = (int16_t) DESCALE(tmp13 + tmp0,
yading@10 936 CONST_BITS+PASS1_BITS+3);
yading@10 937 dataptr[DCTSIZE*4] = (int16_t) DESCALE(tmp13 - tmp0,
yading@10 938 CONST_BITS+PASS1_BITS+3);
yading@10 939
yading@10 940 dataptr++; /* advance pointer to next column */
yading@10 941 }
yading@10 942 }
yading@10 943
yading@10 944 #undef DCTSIZE
yading@10 945 #define DCTSIZE 4
yading@10 946 #define DCTSTRIDE 8
yading@10 947
yading@10 948 void ff_j_rev_dct4(DCTBLOCK data)
yading@10 949 {
yading@10 950 int32_t tmp0, tmp1, tmp2, tmp3;
yading@10 951 int32_t tmp10, tmp11, tmp12, tmp13;
yading@10 952 int32_t z1;
yading@10 953 int32_t d0, d2, d4, d6;
yading@10 954 register int16_t *dataptr;
yading@10 955 int rowctr;
yading@10 956
yading@10 957 /* Pass 1: process rows. */
yading@10 958 /* Note results are scaled up by sqrt(8) compared to a true IDCT; */
yading@10 959 /* furthermore, we scale the results by 2**PASS1_BITS. */
yading@10 960
yading@10 961 data[0] += 4;
yading@10 962
yading@10 963 dataptr = data;
yading@10 964
yading@10 965 for (rowctr = DCTSIZE-1; rowctr >= 0; rowctr--) {
yading@10 966 /* Due to quantization, we will usually find that many of the input
yading@10 967 * coefficients are zero, especially the AC terms. We can exploit this
yading@10 968 * by short-circuiting the IDCT calculation for any row in which all
yading@10 969 * the AC terms are zero. In that case each output is equal to the
yading@10 970 * DC coefficient (with scale factor as needed).
yading@10 971 * With typical images and quantization tables, half or more of the
yading@10 972 * row DCT calculations can be simplified this way.
yading@10 973 */
yading@10 974
yading@10 975 register int *idataptr = (int*)dataptr;
yading@10 976
yading@10 977 d0 = dataptr[0];
yading@10 978 d2 = dataptr[1];
yading@10 979 d4 = dataptr[2];
yading@10 980 d6 = dataptr[3];
yading@10 981
yading@10 982 if ((d2 | d4 | d6) == 0) {
yading@10 983 /* AC terms all zero */
yading@10 984 if (d0) {
yading@10 985 /* Compute a 32 bit value to assign. */
yading@10 986 int16_t dcval = (int16_t) (d0 << PASS1_BITS);
yading@10 987 register int v = (dcval & 0xffff) | ((dcval << 16) & 0xffff0000);
yading@10 988
yading@10 989 idataptr[0] = v;
yading@10 990 idataptr[1] = v;
yading@10 991 }
yading@10 992
yading@10 993 dataptr += DCTSTRIDE; /* advance pointer to next row */
yading@10 994 continue;
yading@10 995 }
yading@10 996
yading@10 997 /* Even part: reverse the even part of the forward DCT. */
yading@10 998 /* The rotator is sqrt(2)*c(-6). */
yading@10 999 if (d6) {
yading@10 1000 if (d2) {
yading@10 1001 /* d0 != 0, d2 != 0, d4 != 0, d6 != 0 */
yading@10 1002 z1 = MULTIPLY(d2 + d6, FIX_0_541196100);
yading@10 1003 tmp2 = z1 + MULTIPLY(-d6, FIX_1_847759065);
yading@10 1004 tmp3 = z1 + MULTIPLY(d2, FIX_0_765366865);
yading@10 1005
yading@10 1006 tmp0 = (d0 + d4) << CONST_BITS;
yading@10 1007 tmp1 = (d0 - d4) << CONST_BITS;
yading@10 1008
yading@10 1009 tmp10 = tmp0 + tmp3;
yading@10 1010 tmp13 = tmp0 - tmp3;
yading@10 1011 tmp11 = tmp1 + tmp2;
yading@10 1012 tmp12 = tmp1 - tmp2;
yading@10 1013 } else {
yading@10 1014 /* d0 != 0, d2 == 0, d4 != 0, d6 != 0 */
yading@10 1015 tmp2 = MULTIPLY(-d6, FIX_1_306562965);
yading@10 1016 tmp3 = MULTIPLY(d6, FIX_0_541196100);
yading@10 1017
yading@10 1018 tmp0 = (d0 + d4) << CONST_BITS;
yading@10 1019 tmp1 = (d0 - d4) << CONST_BITS;
yading@10 1020
yading@10 1021 tmp10 = tmp0 + tmp3;
yading@10 1022 tmp13 = tmp0 - tmp3;
yading@10 1023 tmp11 = tmp1 + tmp2;
yading@10 1024 tmp12 = tmp1 - tmp2;
yading@10 1025 }
yading@10 1026 } else {
yading@10 1027 if (d2) {
yading@10 1028 /* d0 != 0, d2 != 0, d4 != 0, d6 == 0 */
yading@10 1029 tmp2 = MULTIPLY(d2, FIX_0_541196100);
yading@10 1030 tmp3 = MULTIPLY(d2, FIX_1_306562965);
yading@10 1031
yading@10 1032 tmp0 = (d0 + d4) << CONST_BITS;
yading@10 1033 tmp1 = (d0 - d4) << CONST_BITS;
yading@10 1034
yading@10 1035 tmp10 = tmp0 + tmp3;
yading@10 1036 tmp13 = tmp0 - tmp3;
yading@10 1037 tmp11 = tmp1 + tmp2;
yading@10 1038 tmp12 = tmp1 - tmp2;
yading@10 1039 } else {
yading@10 1040 /* d0 != 0, d2 == 0, d4 != 0, d6 == 0 */
yading@10 1041 tmp10 = tmp13 = (d0 + d4) << CONST_BITS;
yading@10 1042 tmp11 = tmp12 = (d0 - d4) << CONST_BITS;
yading@10 1043 }
yading@10 1044 }
yading@10 1045
yading@10 1046 /* Final output stage: inputs are tmp10..tmp13, tmp0..tmp3 */
yading@10 1047
yading@10 1048 dataptr[0] = (int16_t) DESCALE(tmp10, CONST_BITS-PASS1_BITS);
yading@10 1049 dataptr[1] = (int16_t) DESCALE(tmp11, CONST_BITS-PASS1_BITS);
yading@10 1050 dataptr[2] = (int16_t) DESCALE(tmp12, CONST_BITS-PASS1_BITS);
yading@10 1051 dataptr[3] = (int16_t) DESCALE(tmp13, CONST_BITS-PASS1_BITS);
yading@10 1052
yading@10 1053 dataptr += DCTSTRIDE; /* advance pointer to next row */
yading@10 1054 }
yading@10 1055
yading@10 1056 /* Pass 2: process columns. */
yading@10 1057 /* Note that we must descale the results by a factor of 8 == 2**3, */
yading@10 1058 /* and also undo the PASS1_BITS scaling. */
yading@10 1059
yading@10 1060 dataptr = data;
yading@10 1061 for (rowctr = DCTSIZE-1; rowctr >= 0; rowctr--) {
yading@10 1062 /* Columns of zeroes can be exploited in the same way as we did with rows.
yading@10 1063 * However, the row calculation has created many nonzero AC terms, so the
yading@10 1064 * simplification applies less often (typically 5% to 10% of the time).
yading@10 1065 * On machines with very fast multiplication, it's possible that the
yading@10 1066 * test takes more time than it's worth. In that case this section
yading@10 1067 * may be commented out.
yading@10 1068 */
yading@10 1069
yading@10 1070 d0 = dataptr[DCTSTRIDE*0];
yading@10 1071 d2 = dataptr[DCTSTRIDE*1];
yading@10 1072 d4 = dataptr[DCTSTRIDE*2];
yading@10 1073 d6 = dataptr[DCTSTRIDE*3];
yading@10 1074
yading@10 1075 /* Even part: reverse the even part of the forward DCT. */
yading@10 1076 /* The rotator is sqrt(2)*c(-6). */
yading@10 1077 if (d6) {
yading@10 1078 if (d2) {
yading@10 1079 /* d0 != 0, d2 != 0, d4 != 0, d6 != 0 */
yading@10 1080 z1 = MULTIPLY(d2 + d6, FIX_0_541196100);
yading@10 1081 tmp2 = z1 + MULTIPLY(-d6, FIX_1_847759065);
yading@10 1082 tmp3 = z1 + MULTIPLY(d2, FIX_0_765366865);
yading@10 1083
yading@10 1084 tmp0 = (d0 + d4) << CONST_BITS;
yading@10 1085 tmp1 = (d0 - d4) << CONST_BITS;
yading@10 1086
yading@10 1087 tmp10 = tmp0 + tmp3;
yading@10 1088 tmp13 = tmp0 - tmp3;
yading@10 1089 tmp11 = tmp1 + tmp2;
yading@10 1090 tmp12 = tmp1 - tmp2;
yading@10 1091 } else {
yading@10 1092 /* d0 != 0, d2 == 0, d4 != 0, d6 != 0 */
yading@10 1093 tmp2 = MULTIPLY(-d6, FIX_1_306562965);
yading@10 1094 tmp3 = MULTIPLY(d6, FIX_0_541196100);
yading@10 1095
yading@10 1096 tmp0 = (d0 + d4) << CONST_BITS;
yading@10 1097 tmp1 = (d0 - d4) << CONST_BITS;
yading@10 1098
yading@10 1099 tmp10 = tmp0 + tmp3;
yading@10 1100 tmp13 = tmp0 - tmp3;
yading@10 1101 tmp11 = tmp1 + tmp2;
yading@10 1102 tmp12 = tmp1 - tmp2;
yading@10 1103 }
yading@10 1104 } else {
yading@10 1105 if (d2) {
yading@10 1106 /* d0 != 0, d2 != 0, d4 != 0, d6 == 0 */
yading@10 1107 tmp2 = MULTIPLY(d2, FIX_0_541196100);
yading@10 1108 tmp3 = MULTIPLY(d2, FIX_1_306562965);
yading@10 1109
yading@10 1110 tmp0 = (d0 + d4) << CONST_BITS;
yading@10 1111 tmp1 = (d0 - d4) << CONST_BITS;
yading@10 1112
yading@10 1113 tmp10 = tmp0 + tmp3;
yading@10 1114 tmp13 = tmp0 - tmp3;
yading@10 1115 tmp11 = tmp1 + tmp2;
yading@10 1116 tmp12 = tmp1 - tmp2;
yading@10 1117 } else {
yading@10 1118 /* d0 != 0, d2 == 0, d4 != 0, d6 == 0 */
yading@10 1119 tmp10 = tmp13 = (d0 + d4) << CONST_BITS;
yading@10 1120 tmp11 = tmp12 = (d0 - d4) << CONST_BITS;
yading@10 1121 }
yading@10 1122 }
yading@10 1123
yading@10 1124 /* Final output stage: inputs are tmp10..tmp13, tmp0..tmp3 */
yading@10 1125
yading@10 1126 dataptr[DCTSTRIDE*0] = tmp10 >> (CONST_BITS+PASS1_BITS+3);
yading@10 1127 dataptr[DCTSTRIDE*1] = tmp11 >> (CONST_BITS+PASS1_BITS+3);
yading@10 1128 dataptr[DCTSTRIDE*2] = tmp12 >> (CONST_BITS+PASS1_BITS+3);
yading@10 1129 dataptr[DCTSTRIDE*3] = tmp13 >> (CONST_BITS+PASS1_BITS+3);
yading@10 1130
yading@10 1131 dataptr++; /* advance pointer to next column */
yading@10 1132 }
yading@10 1133 }
yading@10 1134
yading@10 1135 void ff_j_rev_dct2(DCTBLOCK data){
yading@10 1136 int d00, d01, d10, d11;
yading@10 1137
yading@10 1138 data[0] += 4;
yading@10 1139 d00 = data[0+0*DCTSTRIDE] + data[1+0*DCTSTRIDE];
yading@10 1140 d01 = data[0+0*DCTSTRIDE] - data[1+0*DCTSTRIDE];
yading@10 1141 d10 = data[0+1*DCTSTRIDE] + data[1+1*DCTSTRIDE];
yading@10 1142 d11 = data[0+1*DCTSTRIDE] - data[1+1*DCTSTRIDE];
yading@10 1143
yading@10 1144 data[0+0*DCTSTRIDE]= (d00 + d10)>>3;
yading@10 1145 data[1+0*DCTSTRIDE]= (d01 + d11)>>3;
yading@10 1146 data[0+1*DCTSTRIDE]= (d00 - d10)>>3;
yading@10 1147 data[1+1*DCTSTRIDE]= (d01 - d11)>>3;
yading@10 1148 }
yading@10 1149
yading@10 1150 void ff_j_rev_dct1(DCTBLOCK data){
yading@10 1151 data[0] = (data[0] + 4)>>3;
yading@10 1152 }
yading@10 1153
yading@10 1154 #undef FIX
yading@10 1155 #undef CONST_BITS