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