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