C++ DESCALE函数代码示例

int main(int argc, char *argv[]) {
  printf("Linear Regression\n");
  printf("=================\n\n");

  // Check for a filename
  if (argc == 2) {
    FILE *input = fopen(argv[1], "r");
    
    if (input == NULL) {
      perror(argv[1]);
      exit(1);

    } else {
      DataSet theData = load_data(input);
      fclose(input);
     
      LinRegResult linReg = linear_regression(theData);
      clean(theData);
      printf("\nSlope   \tm = %15.6e\n", DESCALE(linReg.m)); // print slope
      printf("y-intercept\tb = %15.6e\n", DESCALE(linReg.b)); // print y-intercept
      printf("Correlation\tr = %15.6e\n", DESCALE(linReg.r)); // print correlation
    }

  } else {
    printf("ERROR: Must enter filename after ./linReg\n");

  }
  return 0;
}

int dotProd(int *a, int *b, int n) {
  double dotProd = 0;
  int result = 0;
  int i;
  for (i = 0; i < n; i++) {
    dotProd += DESCALE(a[i]) * DESCALE(b[i]);
  }
  result = dotProd * SCALE;
  return result;
}

ff_fdct248_islow (DCTELEM * data)
{
	int_fast32_t tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7;
	int_fast32_t tmp10, tmp11, tmp12, tmp13;
	int_fast32_t z1;
	DCTELEM *dataptr;
	int ctr;

	row_fdct(data);

	/* Pass 2: process columns.
	 * We remove the PASS1_BITS scaling, but leave the results scaled up
	 * by an overall factor of 8.
	 */

	dataptr = data;
	for (ctr = DCTSIZE-1; ctr >= 0; ctr--)
	{
		tmp0 = dataptr[DCTSIZE*0] + dataptr[DCTSIZE*1];
		tmp1 = dataptr[DCTSIZE*2] + dataptr[DCTSIZE*3];
		tmp2 = dataptr[DCTSIZE*4] + dataptr[DCTSIZE*5];
		tmp3 = dataptr[DCTSIZE*6] + dataptr[DCTSIZE*7];
		tmp4 = dataptr[DCTSIZE*0] - dataptr[DCTSIZE*1];
		tmp5 = dataptr[DCTSIZE*2] - dataptr[DCTSIZE*3];
		tmp6 = dataptr[DCTSIZE*4] - dataptr[DCTSIZE*5];
		tmp7 = dataptr[DCTSIZE*6] - dataptr[DCTSIZE*7];

		tmp10 = tmp0 + tmp3;
		tmp11 = tmp1 + tmp2;
		tmp12 = tmp1 - tmp2;
		tmp13 = tmp0 - tmp3;

		dataptr[DCTSIZE*0] = (DCTELEM) DESCALE(tmp10 + tmp11, PASS1_BITS);
		dataptr[DCTSIZE*4] = (DCTELEM) DESCALE(tmp10 - tmp11, PASS1_BITS);

		z1 = MULTIPLY(tmp12 + tmp13, FIX_0_541196100);
		dataptr[DCTSIZE*2] = (DCTELEM) DESCALE(z1 + MULTIPLY(tmp13, FIX_0_765366865),
		                                       CONST_BITS+PASS1_BITS);
		dataptr[DCTSIZE*6] = (DCTELEM) DESCALE(z1 + MULTIPLY(tmp12, - FIX_1_847759065),
		                                       CONST_BITS+PASS1_BITS);

		tmp10 = tmp4 + tmp7;
		tmp11 = tmp5 + tmp6;
		tmp12 = tmp5 - tmp6;
		tmp13 = tmp4 - tmp7;

		dataptr[DCTSIZE*1] = (DCTELEM) DESCALE(tmp10 + tmp11, PASS1_BITS);
		dataptr[DCTSIZE*5] = (DCTELEM) DESCALE(tmp10 - tmp11, PASS1_BITS);

		z1 = MULTIPLY(tmp12 + tmp13, FIX_0_541196100);
		dataptr[DCTSIZE*3] = (DCTELEM) DESCALE(z1 + MULTIPLY(tmp13, FIX_0_765366865),
		                                       CONST_BITS+PASS1_BITS);
		dataptr[DCTSIZE*7] = (DCTELEM) DESCALE(z1 + MULTIPLY(tmp12, - FIX_1_847759065),
		                                       CONST_BITS+PASS1_BITS);

		dataptr++;                 /* advance pointer to next column */
	}
}

/* Inverse 2-D Discrete Cosine Transform. */
void IDCT(const FBlock * input, PBlock * output)
{
	int Y[64];
	int k, l;

	/* Pass 1: process rows. */
	for (k = 0; k < 8; k++) {

		/* Prescale k-th row: */
		for (l = 0; l < 8; l++)
			Y(k, l) = SCALE(input->block[k][l], S_BITS);

		/* 1-D IDCT on k-th row: */
		idct_1d(&Y(k, 0));
		/* Result Y is scaled up by factor sqrt(8)*2^S_BITS. */
	}

	/* Pass 2: process columns. */
	for (l = 0; l < 8; l++) {
		int Yc[8];

		for (k = 0; k < 8; k++)
			Yc[k] = Y(k, l);
		/* 1-D IDCT on l-th column: */
		idct_1d(Yc);
		/* Result is once more scaled up by a factor sqrt(8). */
		for (k = 0; k < 8; k++) {
			int r = 128 + DESCALE(Yc[k], S_BITS + 3);	/* includes level shift */

			/* Clip to 8 bits unsigned: */
			r = r > 0 ? (r < 255 ? r : 255) : 0;
			X(k, l) = r;
		}
	}
}

static void idct1(int *block)
{
	int i, val;

	val = RANGE(DESCALE(block[0], PASS1_BITS+3));

	for(i=0;i<DCTSIZE2;i++)
		block[i]=val;
}

LinRegResult linear_regression(DataSet theData) {
  
  LinRegResult result;
  int n = theData.n; // number of data points
  double sumx = DESCALE(sum(theData.x, n)); // sum of x
  double sumxx = DESCALE(dotProd(theData.x, theData.x, n)); // sum of each x squared
  double sumy = DESCALE(sum(theData.y, n)); // sum of y
  double sumyy = DESCALE(dotProd(theData.y, theData.y, n)); // sum of each y squared
  double sumxy = DESCALE(dotProd(theData.x, theData.y, n)); // sum of each x * y
  
  double m, b, r;
  // Compute least-squares best fit straight line
  m = (n * sumxy - sumx * sumy) / (n * sumxx - sqr(sumx)); // slope
  b = (sumy * sumxx - (sumx * sumxy)) / (n * sumxx - sqr(sumx)); // y-intercept
  r = (sumxy - sumx * sumy / n) / sqrt((sumxx - sqr(sumx) / n) * (sumyy - sqr(sumy)/ n)); // correlation

  result.m = m * SCALE;
  result.b = b * SCALE;
  result.r = r * SCALE;
  return result;
}

/*
 * Perform dequantization and inverse DCT on one block of coefficients,
 * producing a reduced-size 1x1 output block.
 */
void JPEG_SWIDCT_1X1(int16 *coef_block, uint8 *output_buf, const int32 *quantptr)
{
	int32 dcval;

	/* We hardly need an inverse DCT routine for this: just take the
   * average pixel value, which is one-eighth of the DC coefficient.
   */
	dcval = DEQUANTIZE(coef_block[0], quantptr[0]);
	dcval = (int32) DESCALE((int32) dcval, 3);

	output_buf[0]/*[0]*/ = s_pClip_table[dcval+128];
}

  static void idct(int* tempptr, const jpgd_block_t* dataptr)
  {
    // will be optimized at compile time to either an array access, or 0
    #define ACCESS_COL(x) (((x) < NONZERO_COLS) ? (int)dataptr[x] : 0)
      
    const int z2 = ACCESS_COL(2);
    const int z3 = ACCESS_COL(6);

    const int z1 = MULTIPLY(z2 + z3, FIX_0_541196100);
    const int tmp2 = z1 + MULTIPLY(z3, - FIX_1_847759065);
    const int tmp3 = z1 + MULTIPLY(z2, FIX_0_765366865);

    const int tmp0 = (ACCESS_COL(0) + ACCESS_COL(4)) << CONST_BITS;
    const int tmp1 = (ACCESS_COL(0) - ACCESS_COL(4)) << CONST_BITS;

    const int tmp10 = tmp0 + tmp3;
    const int tmp13 = tmp0 - tmp3;
    const int tmp11 = tmp1 + tmp2;
    const int tmp12 = tmp1 - tmp2;

    const int atmp0 = ACCESS_COL(7);
    const int atmp1 = ACCESS_COL(5);
    const int atmp2 = ACCESS_COL(3);
    const int atmp3 = ACCESS_COL(1);

    const int bz1 = atmp0 + atmp3;
    const int bz2 = atmp1 + atmp2;
    const int bz3 = atmp0 + atmp2;
    const int bz4 = atmp1 + atmp3;
    const int bz5 = MULTIPLY(bz3 + bz4, FIX_1_175875602);
       
    const int az1 = MULTIPLY(bz1, - FIX_0_899976223);
    const int az2 = MULTIPLY(bz2, - FIX_2_562915447);
    const int az3 = MULTIPLY(bz3, - FIX_1_961570560) + bz5;
    const int az4 = MULTIPLY(bz4, - FIX_0_390180644) + bz5;
    
    const int btmp0 = MULTIPLY(atmp0, FIX_0_298631336) + az1 + az3;
    const int btmp1 = MULTIPLY(atmp1, FIX_2_053119869) + az2 + az4;
    const int btmp2 = MULTIPLY(atmp2, FIX_3_072711026) + az2 + az3;
    const int btmp3 = MULTIPLY(atmp3, FIX_1_501321110) + az1 + az4;

    tempptr[0] = DESCALE(tmp10 + btmp3, CONST_BITS-PASS1_BITS);
    tempptr[7] = DESCALE(tmp10 - btmp3, CONST_BITS-PASS1_BITS);
    tempptr[1] = DESCALE(tmp11 + btmp2, CONST_BITS-PASS1_BITS);
    tempptr[6] = DESCALE(tmp11 - btmp2, CONST_BITS-PASS1_BITS);
    tempptr[2] = DESCALE(tmp12 + btmp1, CONST_BITS-PASS1_BITS);
    tempptr[5] = DESCALE(tmp12 - btmp1, CONST_BITS-PASS1_BITS);
    tempptr[3] = DESCALE(tmp13 + btmp0, CONST_BITS-PASS1_BITS);
    tempptr[4] = DESCALE(tmp13 - btmp0, CONST_BITS-PASS1_BITS);
  }

static void dct_1d(dct_data_t *src, dct_data_t *dst)
{
	unsigned int k, n;
	int tmp;
	const dct_data_t dct_coeff_table[DCT_SIZE][DCT_SIZE] = {
#include "dct_coeff_table.txt"
	};

	for (k = 0; k < DCT_SIZE; k++) {
		for(n = 0, tmp = 0; n < DCT_SIZE; n++) {
			int coeff = (int)dct_coeff_table[k][n];
			tmp += src[n] * coeff;
		}
		dst[k] = DESCALE(tmp, CONST_BITS);
	}
}

openexif_jpeg_idct_1x1 (oe_j_decompress_ptr cinfo, openexif_jpeg_component_info * compptr,
	       OE_JCOEFPTR coef_block,
	       OE_JSAMPARRAY output_buf, OE_JDIMENSION output_col)
{
  int dcval;
  ISLOW_MULT_TYPE * quantptr;
  OE_JSAMPLE *range_limit = IDCT_range_limit(cinfo);
  SHIFT_TEMPS

  /* We hardly need an inverse DCT routine for this: just take the
   * average pixel value, which is one-eighth of the DC coefficient.
   */
  quantptr = (ISLOW_MULT_TYPE *) compptr->dct_table;
  dcval = DEQUANTIZE(coef_block[0], quantptr[0]);
  dcval = (int) DESCALE((INT32) dcval, 3);

  output_buf[0][output_col] = range_limit[dcval & RANGE_MASK];
}

METHODDEF void
start_pass_fdctmgr (j_compress_ptr cinfo)
{
  my_fdct_ptr fdct = (my_fdct_ptr) cinfo->fdct;
  int ci, qtblno, i;
  jpeg_component_info *compptr;
  JQUANT_TBL * qtbl;
#ifdef DCT_ISLOW_SUPPORTED
  DCTELEM * dtbl;
#endif

  for (ci = 0, compptr = cinfo->comp_info; ci < cinfo->num_components;
       ci++, compptr++) {
    qtblno = compptr->quant_tbl_no;
    /* Make sure specified quantization table is present */
    if (qtblno < 0 || qtblno >= NUM_QUANT_TBLS ||
	cinfo->quant_tbl_ptrs[qtblno] == NULL)
      ERREXIT1(cinfo, JERR_NO_QUANT_TABLE, qtblno);
    qtbl = cinfo->quant_tbl_ptrs[qtblno];
    /* Compute divisors for this quant table */
    /* We may do this more than once for same table, but it's not a big deal */
    switch (cinfo->dct_method) {
#ifdef DCT_ISLOW_SUPPORTED
    case JDCT_ISLOW:
      /* For LL&M IDCT method, divisors are equal to raw quantization
       * coefficients multiplied by 8 (to counteract scaling).
       */
      if (fdct->divisors[qtblno] == NULL) {
	fdct->divisors[qtblno] = (DCTELEM *)
	  (*cinfo->mem->alloc_small) ((j_common_ptr) cinfo, JPOOL_IMAGE,
				      DCTSIZE2 * SIZEOF(DCTELEM));
      }
      dtbl = fdct->divisors[qtblno];
      for (i = 0; i < DCTSIZE2; i++) {
	dtbl[i] = ((DCTELEM) qtbl->quantval[jpeg_zigzag_order[i]]) << 3;
      }
      break;
#endif
#ifdef DCT_IFAST_SUPPORTED
    case JDCT_IFAST:
      {
	/* For AA&N IDCT method, divisors are equal to quantization
	 * coefficients scaled by scalefactor[row]*scalefactor[col], where
	 *   scalefactor[0] = 1
	 *   scalefactor[k] = cos(k*PI/16) * sqrt(2)    for k=1..7
	 * We apply a further scale factor of 8.
	 */
#define CONST_BITS 14
	static const INT16 aanscales[DCTSIZE2] = {
	  /* precomputed values scaled up by 14 bits: in natural order */
	  16384, 22725, 21407, 19266, 16384, 12873,  8867,  4520,
	  22725, 31521, 29692, 26722, 22725, 17855, 12299,  6270,
	  21407, 29692, 27969, 25172, 21407, 16819, 11585,  5906,
	  19266, 26722, 25172, 22654, 19266, 15137, 10426,  5315,
	  16384, 22725, 21407, 19266, 16384, 12873,  8867,  4520,
	  12873, 17855, 16819, 15137, 12873, 10114,  6967,  3552,
	   8867, 12299, 11585, 10426,  8867,  6967,  4799,  2446,
	   4520,  6270,  5906,  5315,  4520,  3552,  2446,  1247
	};
	SHIFT_TEMPS

	if (fdct->divisors[qtblno] == NULL) {
	  fdct->divisors[qtblno] = (DCTELEM *)
	    (*cinfo->mem->alloc_small) ((j_common_ptr) cinfo, JPOOL_IMAGE,
					DCTSIZE2 * SIZEOF(DCTELEM));
	}
	dtbl = fdct->divisors[qtblno];
	for (i = 0; i < DCTSIZE2; i++) {
	  dtbl[i] = (DCTELEM)
	    DESCALE(MULTIPLY16V16((INT32) qtbl->quantval[jpeg_zigzag_order[i]],
				  (INT32) aanscales[i]),
		    CONST_BITS-3);
	}
      }
      break;
#endif
#ifdef DCT_FLOAT_SUPPORTED
    case JDCT_FLOAT:
      {
	/* For float AA&N IDCT method, divisors are equal to quantization
	 * coefficients scaled by scalefactor[row]*scalefactor[col], where
	 *   scalefactor[0] = 1
	 *   scalefactor[k] = cos(k*PI/16) * sqrt(2)    for k=1..7
	 * We apply a further scale factor of 8.
	 * What's actually stored is 1/divisor so that the inner loop can
	 * use a multiplication rather than a division.
	 */
	FAST_FLOAT * fdtbl;
	int row, col;
	static const double aanscalefactor[DCTSIZE] = {
	  1.0, 1.387039845, 1.306562965, 1.175875602,
	  1.0, 0.785694958, 0.541196100, 0.275899379
	};

	if (fdct->float_divisors[qtblno] == NULL) {
	  fdct->float_divisors[qtblno] = (FAST_FLOAT *)
	    (*cinfo->mem->alloc_small) ((j_common_ptr) cinfo, JPOOL_IMAGE,
					DCTSIZE2 * SIZEOF(FAST_FLOAT));
	}
	fdtbl = fdct->float_divisors[qtblno];
//.........这里部分代码省略.........

static av_always_inline void row_fdct(DCTELEM * data)
{
	int_fast32_t tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7;
	int_fast32_t tmp10, tmp11, tmp12, tmp13;
	int_fast32_t z1, z2, z3, z4, z5;
	DCTELEM *dataptr;
	int ctr;

	/* Pass 1: process rows. */
	/* Note results are scaled up by sqrt(8) compared to a true DCT; */
	/* furthermore, we scale the results by 2**PASS1_BITS. */

	dataptr = data;
	for (ctr = DCTSIZE-1; ctr >= 0; ctr--)
	{
		tmp0 = dataptr[0] + dataptr[7];
		tmp7 = dataptr[0] - dataptr[7];
		tmp1 = dataptr[1] + dataptr[6];
		tmp6 = dataptr[1] - dataptr[6];
		tmp2 = dataptr[2] + dataptr[5];
		tmp5 = dataptr[2] - dataptr[5];
		tmp3 = dataptr[3] + dataptr[4];
		tmp4 = dataptr[3] - dataptr[4];

		/* Even part per LL&M figure 1 --- note that published figure is faulty;
		 * rotator "sqrt(2)*c1" should be "sqrt(2)*c6".
		 */

		tmp10 = tmp0 + tmp3;
		tmp13 = tmp0 - tmp3;
		tmp11 = tmp1 + tmp2;
		tmp12 = tmp1 - tmp2;

		dataptr[0] = (DCTELEM) ((tmp10 + tmp11) << PASS1_BITS);
		dataptr[4] = (DCTELEM) ((tmp10 - tmp11) << PASS1_BITS);

		z1 = MULTIPLY(tmp12 + tmp13, FIX_0_541196100);
		dataptr[2] = (DCTELEM) DESCALE(z1 + MULTIPLY(tmp13, FIX_0_765366865),
		                               CONST_BITS-PASS1_BITS);
		dataptr[6] = (DCTELEM) DESCALE(z1 + MULTIPLY(tmp12, - FIX_1_847759065),
		                               CONST_BITS-PASS1_BITS);

		/* Odd part per figure 8 --- note paper omits factor of sqrt(2).
		 * cK represents cos(K*pi/16).
		 * i0..i3 in the paper are tmp4..tmp7 here.
		 */

		z1 = tmp4 + tmp7;
		z2 = tmp5 + tmp6;
		z3 = tmp4 + tmp6;
		z4 = tmp5 + tmp7;
		z5 = MULTIPLY(z3 + z4, FIX_1_175875602); /* sqrt(2) * c3 */

		tmp4 = MULTIPLY(tmp4, FIX_0_298631336); /* sqrt(2) * (-c1+c3+c5-c7) */
		tmp5 = MULTIPLY(tmp5, FIX_2_053119869); /* sqrt(2) * ( c1+c3-c5+c7) */
		tmp6 = MULTIPLY(tmp6, FIX_3_072711026); /* sqrt(2) * ( c1+c3+c5-c7) */
		tmp7 = MULTIPLY(tmp7, FIX_1_501321110); /* sqrt(2) * ( c1+c3-c5-c7) */
		z1 = MULTIPLY(z1, - FIX_0_899976223); /* sqrt(2) * (c7-c3) */
		z2 = MULTIPLY(z2, - FIX_2_562915447); /* sqrt(2) * (-c1-c3) */
		z3 = MULTIPLY(z3, - FIX_1_961570560); /* sqrt(2) * (-c3-c5) */
		z4 = MULTIPLY(z4, - FIX_0_390180644); /* sqrt(2) * (c5-c3) */

		z3 += z5;
		z4 += z5;

		dataptr[7] = (DCTELEM) DESCALE(tmp4 + z1 + z3, CONST_BITS-PASS1_BITS);
		dataptr[5] = (DCTELEM) DESCALE(tmp5 + z2 + z4, CONST_BITS-PASS1_BITS);
		dataptr[3] = (DCTELEM) DESCALE(tmp6 + z2 + z3, CONST_BITS-PASS1_BITS);
		dataptr[1] = (DCTELEM) DESCALE(tmp7 + z1 + z4, CONST_BITS-PASS1_BITS);

		dataptr += DCTSIZE;         /* advance pointer to next row */
	}
}

//.........这里部分代码省略.........
                *--fblock = fr2;
                *--fblock = fr1;
                *--fblock = fr0;
                fblock+=8+4;
        } while(--i);
        block-=8*8;
        fblock-=8*8+4;

        load_matrix(odd_table);

        i = 8;

//        ofs1 = sizeof(float)*1;
//        ofs2 = sizeof(float)*2;
//        ofs3 = sizeof(float)*3;

        do {
                float t0,t1,t2,t3;
                fr0 = block[1];
                fr1 = block[3];
                fr2 = block[5];
                fr3 = block[7];
                block+=8;
                ftrv();
                t0 = *fblock++;
                t1 = *fblock++;
                t2 = *fblock++;
                t3 = *fblock++;
                fblock+=4;
                *--fblock = t0 - fr0;
                *--fblock = t1 - fr1;
                *--fblock = t2 - fr2;
                *--fblock = t3 - fr3;
                *--fblock = t3 + fr3;
                *--fblock = t2 + fr2;
                *--fblock = t1 + fr1;
                *--fblock = t0 + fr0;
                fblock+=8;
        } while(--i);
        block-=8*8;
        fblock-=8*8;

        /* col */

        /* even part */
        load_matrix(even_table);

        ofs1 = sizeof(float)*2*8;
        ofs2 = sizeof(float)*4*8;
        ofs3 = sizeof(float)*6*8;

        i = 8;

#define        OA(fblock,ofs)   *(float*)((char*)fblock + ofs)

        do {
                fr0 = OA(fblock,   0);
                fr1 = OA(fblock,ofs1);
                fr2 = OA(fblock,ofs2);
                fr3 = OA(fblock,ofs3);
                ftrv();
                OA(fblock,0   ) = fr0;
                OA(fblock,ofs1) = fr1;
                OA(fblock,ofs2) = fr2;
                OA(fblock,ofs3) = fr3;
                fblock++;
        } while(--i);
        fblock-=8;

        load_matrix(odd_table);

        i=8;
        do {
                float t0,t1,t2,t3;
                t0 = OA(fblock,   0); /* [8*0] */
                t1 = OA(fblock,ofs1); /* [8*2] */
                t2 = OA(fblock,ofs2); /* [8*4] */
                t3 = OA(fblock,ofs3); /* [8*6] */
                fblock+=8;
                fr0 = OA(fblock,   0); /* [8*1] */
                fr1 = OA(fblock,ofs1); /* [8*3] */
                fr2 = OA(fblock,ofs2); /* [8*5] */
                fr3 = OA(fblock,ofs3); /* [8*7] */
                fblock+=-8+1;
                ftrv();
                block[8*0] = DESCALE(t0 + fr0,3);
                block[8*7] = DESCALE(t0 - fr0,3);
                block[8*1] = DESCALE(t1 + fr1,3);
                block[8*6] = DESCALE(t1 - fr1,3);
                block[8*2] = DESCALE(t2 + fr2,3);
                block[8*5] = DESCALE(t2 - fr2,3);
                block[8*3] = DESCALE(t3 + fr3,3);
                block[8*4] = DESCALE(t3 - fr3,3);
                block++;
        } while(--i);

#if defined(__SH4__)
#error  "FIXME!! change to double"
#endif
}

void j_rev_dct(DCTBLOCK data)
{
  INT32 tmp0, tmp1, tmp2, tmp3;
  INT32 tmp10, tmp11, tmp12, tmp13;
  INT32 z1, z2, z3, z4, z5;
  register DCTELEM *dataptr;
  int rowctr;

  /* Pass 1: process rows. */
  /* Note results are scaled up by sqrt(8) compared to a true IDCT; */
  /* furthermore, we scale the results by 2**PASS1_BITS. */

  dataptr = data;
  for (rowctr = DCTSIZE-1; rowctr >= 0; rowctr--) 
  {
    /* Even part: reverse the even part of the forward DCT. */
    /* The rotator is sqrt(2)*c(-6). */

    z2 = (INT32) dataptr[2];
    z3 = (INT32) dataptr[6];

    z1 = MULTIPLY(z2 + z3, FIX_0_541196100);
    tmp2 = z1 + MULTIPLY(z3, - FIX_1_847759065);
    tmp3 = z1 + MULTIPLY(z2, FIX_0_765366865);

    tmp0 = ((INT32) dataptr[0] + (INT32) dataptr[4]) << CONST_BITS;
    tmp1 = ((INT32) dataptr[0] - (INT32) dataptr[4]) << CONST_BITS;

    tmp10 = tmp0 + tmp3;
    tmp13 = tmp0 - tmp3;
    tmp11 = tmp1 + tmp2;
    tmp12 = tmp1 - tmp2;
    
    /* Odd part per figure 8; the matrix is unitary and hence its
     * transpose is its inverse.  i0..i3 are y7,y5,y3,y1 respectively.
     */

    tmp0 = (INT32) dataptr[7];
    tmp1 = (INT32) dataptr[5];
    tmp2 = (INT32) dataptr[3];
    tmp3 = (INT32) dataptr[1];

    z1 = tmp0 + tmp3;
    z2 = tmp1 + tmp2;
    z3 = tmp0 + tmp2;
    z4 = tmp1 + tmp3;
    z5 = MULTIPLY(z3 + z4, FIX_1_175875602); /* sqrt(2) * c3 */
    
    tmp0 = MULTIPLY(tmp0, FIX_0_298631336); /* sqrt(2) * (-c1+c3+c5-c7) */
    tmp1 = MULTIPLY(tmp1, FIX_2_053119869); /* sqrt(2) * ( c1+c3-c5+c7) */
    tmp2 = MULTIPLY(tmp2, FIX_3_072711026); /* sqrt(2) * ( c1+c3+c5-c7) */
    tmp3 = MULTIPLY(tmp3, FIX_1_501321110); /* sqrt(2) * ( c1+c3-c5-c7) */
    z1 = MULTIPLY(z1, - FIX_0_899976223); /* sqrt(2) * (c7-c3) */
    z2 = MULTIPLY(z2, - FIX_2_562915447); /* sqrt(2) * (-c1-c3) */
    z3 = MULTIPLY(z3, - FIX_1_961570560); /* sqrt(2) * (-c3-c5) */
    z4 = MULTIPLY(z4, - FIX_0_390180644); /* sqrt(2) * (c5-c3) */
    
    z3 += z5;
    z4 += z5;
    
    tmp0 += z1 + z3;
    tmp1 += z2 + z4;
    tmp2 += z2 + z3;
    tmp3 += z1 + z4;

    /* Final output stage: inputs are tmp10..tmp13, tmp0..tmp3 */

    dataptr[0] = (DCTELEM) DESCALE(tmp10 + tmp3, CONST_BITS-PASS1_BITS);
    dataptr[7] = (DCTELEM) DESCALE(tmp10 - tmp3, CONST_BITS-PASS1_BITS);
    dataptr[1] = (DCTELEM) DESCALE(tmp11 + tmp2, CONST_BITS-PASS1_BITS);
    dataptr[6] = (DCTELEM) DESCALE(tmp11 - tmp2, CONST_BITS-PASS1_BITS);
    dataptr[2] = (DCTELEM) DESCALE(tmp12 + tmp1, CONST_BITS-PASS1_BITS);
    dataptr[5] = (DCTELEM) DESCALE(tmp12 - tmp1, CONST_BITS-PASS1_BITS);
    dataptr[3] = (DCTELEM) DESCALE(tmp13 + tmp0, CONST_BITS-PASS1_BITS);
    dataptr[4] = (DCTELEM) DESCALE(tmp13 - tmp0, CONST_BITS-PASS1_BITS);

    dataptr += DCTSIZE;		/* advance pointer to next row */
  }

  /* Pass 2: process columns. */
  /* Note that we must descale the results by a factor of 8 == 2**3, */
  /* and also undo the PASS1_BITS scaling. */

  dataptr = data;
  for (rowctr = DCTSIZE-1; rowctr >= 0; rowctr--) {
    /* Even part: reverse the even part of the forward DCT. */
    /* The rotator is sqrt(2)*c(-6). */

    z2 = (INT32) dataptr[DCTSIZE*2];
    z3 = (INT32) dataptr[DCTSIZE*6];

    z1 = MULTIPLY(z2 + z3, FIX_0_541196100);
    tmp2 = z1 + MULTIPLY(z3, - FIX_1_847759065);
    tmp3 = z1 + MULTIPLY(z2, FIX_0_765366865);

    tmp0 = ((INT32) dataptr[DCTSIZE*0] + (INT32) dataptr[DCTSIZE*4]) << CONST_BITS;
    tmp1 = ((INT32) dataptr[DCTSIZE*0] - (INT32) dataptr[DCTSIZE*4]) << CONST_BITS;

    tmp10 = tmp0 + tmp3;
    tmp13 = tmp0 - tmp3;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
                    tmp2 = MULTIPLY(-d6, FIX_1_306562965);
                    tmp3 = MULTIPLY(d6, FIX_0_541196100);

                    tmp0 = (d0 + d4) << CONST_BITS;
                    tmp1 = (d0 - d4) << CONST_BITS;

                    tmp10 = tmp0 + tmp3;
                    tmp13 = tmp0 - tmp3;
                    tmp11 = tmp1 + tmp2;
                    tmp12 = tmp1 - tmp2;
            }
    } else {
            if (d2) {
                    /* d0 != 0, d2 != 0, d4 != 0, d6 == 0 */
                    tmp2 = MULTIPLY(d2, FIX_0_541196100);
                    tmp3 = MULTIPLY(d2, FIX_1_306562965);

                    tmp0 = (d0 + d4) << CONST_BITS;
                    tmp1 = (d0 - d4) << CONST_BITS;

                    tmp10 = tmp0 + tmp3;
                    tmp13 = tmp0 - tmp3;
                    tmp11 = tmp1 + tmp2;
                    tmp12 = tmp1 - tmp2;
            } else {
                    /* d0 != 0, d2 == 0, d4 != 0, d6 == 0 */
                    tmp10 = tmp13 = (d0 + d4) << CONST_BITS;
                    tmp11 = tmp12 = (d0 - d4) << CONST_BITS;
            }
      }

    /* Final output stage: inputs are tmp10..tmp13, tmp0..tmp3 */

    dataptr[0] = (DCTELEM) DESCALE(tmp10, CONST_BITS-PASS1_BITS);
    dataptr[1] = (DCTELEM) DESCALE(tmp11, CONST_BITS-PASS1_BITS);
    dataptr[2] = (DCTELEM) DESCALE(tmp12, CONST_BITS-PASS1_BITS);
    dataptr[3] = (DCTELEM) DESCALE(tmp13, CONST_BITS-PASS1_BITS);

    dataptr += DCTSTRIDE;       /* advance pointer to next row */
  }

  /* Pass 2: process columns. */
  /* Note that we must descale the results by a factor of 8 == 2**3, */
  /* and also undo the PASS1_BITS scaling. */

  dataptr = data;
  for (rowctr = DCTSIZE-1; rowctr >= 0; rowctr--) {
    /* Columns of zeroes can be exploited in the same way as we did with rows.
     * However, the row calculation has created many nonzero AC terms, so the
     * simplification applies less often (typically 5% to 10% of the time).
     * On machines with very fast multiplication, it's possible that the
     * test takes more time than it's worth.  In that case this section
     * may be commented out.
     */

    d0 = dataptr[DCTSTRIDE*0];
    d2 = dataptr[DCTSTRIDE*1];
    d4 = dataptr[DCTSTRIDE*2];
    d6 = dataptr[DCTSTRIDE*3];

    /* Even part: reverse the even part of the forward DCT. */
    /* The rotator is sqrt(2)*c(-6). */
    if (d6) {
            if (d2) {
                    /* d0 != 0, d2 != 0, d4 != 0, d6 != 0 */
                    z1 = MULTIPLY(d2 + d6, FIX_0_541196100);