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