text: 
  /*
   * mfwddct.c (derived from jfwddct.c, which carries the following info)
   *
   * Copyright (C) 1991, 1992, Thomas G. Lane. This file is part of the
   * Independent JPEG Group's software. For conditions of distribution and use,
   * see the accompanying README file.
   *
   * This file contains the basic DCT (Discrete Cosine Transform) transformation
   * subroutine.
   *
   * This implementation is based on Appendix A.2 of the book "Discrete Cosine
   * Transform---Algorithms, Advantages, Applications" by K.R. Rao and P. Yip
   * (Academic Press, Inc, London, 1990). It uses scaled fixed-point arithmetic
   * instead of floating point.
   */
  
  /*************************************************************************
                               R C S Information
  **************************************************************************/
  
  /* $Log: mfwddct.c,v $
   * Revision 1.2  1994/01/12  19:07:46  daf
   * Added RCS info.
   *
  
  * revision 1.1    locked by: daf;
  * date: 1994/01/07 17:05:21;  author: daf;  state: Exp;
  * Initial revision
  */
  
  #include "all.h"
  
  #include "dct.h"
  #include "mtypes.h"
  
  /*
   * The poop on this scaling stuff is as follows:
   *
   * We have to do addition and subtraction of the integer inputs, which is no
   * problem, and multiplication by fractional constants, which is a problem to
   * do in integer arithmetic.  We multiply all the constants by DCT_SCALE and
   * convert them to integer constants (thus retaining LG2_DCT_SCALE bits of
   * precision in the constants).  After doing a multiplication we have to
   * divide the product by DCT_SCALE, with proper rounding, to produce the
   * correct output.  The division can be implemented cheaply as a right shift
   * of LG2_DCT_SCALE bits.  The DCT equations also specify an additional
   * division by 2 on the final outputs; this can be folded into the
   * right-shift by shifting one more bit (see UNFIXH).
   *
   * If you are planning to recode this in assembler, you might want to set
   * LG2_DCT_SCALE to 15.  This loses a bit of precision, but then all the
   * multiplications are between 16-bit quantities (given 8-bit JSAMPLEs!) so
   * you could use a signed 16x16=>32 bit multiply instruction instead of full
   * 32x32 multiply.  Unfortunately there's no way to describe such a multiply
   * portably in C, so we've gone for the extra bit of accuracy here.
   */
  
  #define EIGHT_BIT_SAMPLES
  #ifdef EIGHT_BIT_SAMPLES
  #define LG2_DCT_SCALE 16
  #else
  #define LG2_DCT_SCALE 15	/* lose a little precision to avoid overflow */
  #endif
  
  #define ONE	((int32) 1)
  
  #define DCT_SCALE (ONE << LG2_DCT_SCALE)
  
  /* In some places we shift the inputs left by a couple more bits, */
  /* so that they can be added to fractional results without too much */
  /* loss of precision. */
  #define LG2_OVERSCALE 2
  #define OVERSCALE  (ONE << LG2_OVERSCALE)
  #define OVERSHIFT(x)  ((x) <<= LG2_OVERSCALE)
  
  /* Scale a fractional constant by DCT_SCALE */
  #define FIX(x)	((int32) ((x) * DCT_SCALE + 0.5))
  
  /* Scale a fractional constant by DCT_SCALE/OVERSCALE */
  /* Such a constant can be multiplied with an overscaled input */
  /* to produce something that's scaled by DCT_SCALE */
  #define FIXO(x)  ((int32) ((x) * DCT_SCALE / OVERSCALE + 0.5))
  
  /* Descale and correctly round a value that's scaled by DCT_SCALE */
  #define UNFIX(x)   RIGHT_SHIFT((x) + (ONE << (LG2_DCT_SCALE-1)), LG2_DCT_SCALE)
  
  /* Same with an additional division by 2, ie, correctly rounded UNFIX(x/2) */
  #define UNFIXH(x)  RIGHT_SHIFT((x) + (ONE << LG2_DCT_SCALE), LG2_DCT_SCALE+1)
  
  /* Take a value scaled by DCT_SCALE and round to integer scaled by OVERSCALE */
  #define UNFIXO(x)  RIGHT_SHIFT((x) + (ONE << (LG2_DCT_SCALE-1-LG2_OVERSCALE)),\
  			       LG2_DCT_SCALE-LG2_OVERSCALE)
  
  /* Here are the constants we need */
  /* SIN_i_j is sine of i*pi/j, scaled by DCT_SCALE */
  /* COS_i_j is cosine of i*pi/j, scaled by DCT_SCALE */
  
  #define SIN_1_4 FIX(0.707106781)
  #define COS_1_4 SIN_1_4
  
  #define SIN_1_8 FIX(0.382683432)
  #define COS_1_8 FIX(0.923879533)
  #define SIN_3_8 COS_1_8
  #define COS_3_8 SIN_1_8
  
  #define SIN_1_16 FIX(0.195090322)
  #define COS_1_16 FIX(0.980785280)
  #define SIN_7_16 COS_1_16
  #define COS_7_16 SIN_1_16
  
  #define SIN_3_16 FIX(0.555570233)
  #define COS_3_16 FIX(0.831469612)
  #define SIN_5_16 COS_3_16
  #define COS_5_16 SIN_3_16
  
  /* OSIN_i_j is sine of i*pi/j, scaled by DCT_SCALE/OVERSCALE */
  /* OCOS_i_j is cosine of i*pi/j, scaled by DCT_SCALE/OVERSCALE */
  
  #define OSIN_1_4 FIXO(0.707106781)
  #define OCOS_1_4 OSIN_1_4
  
  #define OSIN_1_8 FIXO(0.382683432)
  #define OCOS_1_8 FIXO(0.923879533)
  #define OSIN_3_8 OCOS_1_8
  #define OCOS_3_8 OSIN_1_8
  
  #define OSIN_1_16 FIXO(0.195090322)
  #define OCOS_1_16 FIXO(0.980785280)
  #define OSIN_7_16 OCOS_1_16
  #define OCOS_7_16 OSIN_1_16
  
  #define OSIN_3_16 FIXO(0.555570233)
  #define OCOS_3_16 FIXO(0.831469612)
  #define OSIN_5_16 OCOS_3_16
  #define OCOS_5_16 OSIN_3_16
  
  
  
  
  /*
   * --------------------------------------------------------------
   *
   * mp_fwd_dct_block --
   *
   * Perform the forward DCT on one block of samples.
   *
   * A 2-D DCT can be done by 1-D DCT on each row followed by 1-D DCT on each
   * column.
   *
   * Results: None
   *
   * Side effects: Overwrites the input data
   *
   * --------------------------------------------------------------
   */
  
  void
  mp_fwd_dct_block(data2d)
      Block data2d;
  {
      int16 *data = (int16 *) data2d;	/* this algorithm wants
  					 * a 1-d array */
      int pass, rowctr;
      register int16 *inptr, *outptr;
      int16 workspace[DCTSIZE_SQ];
      SHIFT_TEMPS
  
  #ifdef ndef
      {
  	int y;
  
  	printf("fwd_dct (beforehand):\n");
  	for (y = 0; y < 8; y++)
  	    printf("%4d %4d %4d %4d %4d %4d %4d %4d\n",
  		   data2d[y][0], data2d[y][1],
  		   data2d[y][2], data2d[y][3],
  		   data2d[y][4], data2d[y][5],
  		   data2d[y][6], data2d[y][7]);
      }
  #endif
  
      /*
       * Each iteration of the inner loop performs one 8-point 1-D DCT. It
       * reads from a *row* of the input matrix and stores into a *column*
       * of the output matrix.  In the first pass, we read from the data[]
       * array and store into the local workspace[].  In the second pass,
       * we read from the workspace[] array and store into data[], thus
       * performing the equivalent of a columnar DCT pass with no variable
       * array indexing.
       */
  
      inptr = data;		/* initialize pointers for first pass */
      outptr = workspace;
      for (pass = 1; pass >= 0; pass--) {
  	for (rowctr = DCTSIZE - 1; rowctr >= 0; rowctr--) {
  	    /*
  	     * many tmps have nonoverlapping lifetime -- flashy
  	     * register colourers should be able to do this lot
  	     * very well
  	     */
  	    int32 tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7;
  	    int32 tmp10, tmp11, tmp12, tmp13;
  	    int32 tmp14, tmp15, tmp16, tmp17;
  	    int32 tmp25, tmp26;
  	    /* SHIFT_TEMPS */
  
  	    /* temp0 through tmp7:  -512 to +512 */
  	    /* if I-block, then -256 to +256 */
  	    tmp0 = inptr[7] + inptr[0];
  	    tmp1 = inptr[6] + inptr[1];
  	    tmp2 = inptr[5] + inptr[2];
  	    tmp3 = inptr[4] + inptr[3];
  	    tmp4 = inptr[3] - inptr[4];
  	    tmp5 = inptr[2] - inptr[5];
  	    tmp6 = inptr[1] - inptr[6];
  	    tmp7 = inptr[0] - inptr[7];
  
  	    /* tmp10 through tmp13:  -1024 to +1024 */
  	    /* if I-block, then -512 to +512 */
  	    tmp10 = tmp3 + tmp0;
  	    tmp11 = tmp2 + tmp1;
  	    tmp12 = tmp1 - tmp2;
  	    tmp13 = tmp0 - tmp3;
  
  	    outptr[0] = (int16) UNFIXH((tmp10 + tmp11) * SIN_1_4);
  	    outptr[DCTSIZE * 4] = (int16) UNFIXH((tmp10 - tmp11) * COS_1_4);
  
  	    outptr[DCTSIZE * 2] = (int16) UNFIXH(tmp13 * COS_1_8 + tmp12 * SIN_1_8);
  	    outptr[DCTSIZE * 6] = (int16) UNFIXH(tmp13 * SIN_1_8 - tmp12 * COS_1_8);
  
  	    tmp16 = UNFIXO((tmp6 + tmp5) * SIN_1_4);
  	    tmp15 = UNFIXO((tmp6 - tmp5) * COS_1_4);
  
  	    OVERSHIFT(tmp4);
  	    OVERSHIFT(tmp7);
  
  	    /*
  	     * tmp4, tmp7, tmp15, tmp16 are overscaled by
  	     * OVERSCALE
  	     */
  
  	    tmp14 = tmp4 + tmp15;
  	    tmp25 = tmp4 - tmp15;
  	    tmp26 = tmp7 - tmp16;
  	    tmp17 = tmp7 + tmp16;
  
  	    outptr[DCTSIZE] = (int16) UNFIXH(tmp17 * OCOS_1_16 + tmp14 * OSIN_1_16);
  	    outptr[DCTSIZE * 7] = (int16) UNFIXH(tmp17 * OCOS_7_16 - tmp14 * OSIN_7_16);
  	    outptr[DCTSIZE * 5] = (int16) UNFIXH(tmp26 * OCOS_5_16 + tmp25 * OSIN_5_16);
  	    outptr[DCTSIZE * 3] = (int16) UNFIXH(tmp26 * OCOS_3_16 - tmp25 * OSIN_3_16);
  
  	    inptr += DCTSIZE;	/* advance inptr to next row */
  	    outptr++;		/* advance outptr to next column */
  	}
  	/* end of pass; in case it was pass 1, set up for pass 2 */
  	inptr = workspace;
  	outptr = data;
      }
  #ifdef ndef
      {
  	int y;
  
  	printf("fwd_dct (afterward):\n");
  	for (y = 0; y < 8; y++)
  	    printf("%4d %4d %4d %4d %4d %4d %4d %4d\n",
  		   data2d[y][0], data2d[y][1],
  		   data2d[y][2], data2d[y][3],
  		   data2d[y][4], data2d[y][5],
  		   data2d[y][6], data2d[y][7]);
      }
  #endif
  }
  