00001 /*
00002  * jrevdct.c
00003  *
00004  * Copyright (C) 1991, 1992, Thomas G. Lane.
00005  * This file is part of the Independent JPEG Group's software.
00006  * For conditions of distribution and use, see the accompanying README file.
00007  *
00008  * This file contains the basic inverse-DCT transformation subroutine.
00009  *
00010  * This implementation is based on an algorithm described in
00011  *   C. Loeffler, A. Ligtenberg and G. Moschytz, "Practical Fast 1-D DCT
00012  *   Algorithms with 11 Multiplications", Proc. Int'l. Conf. on Acoustics,
00013  *   Speech, and Signal Processing 1989 (ICASSP '89), pp. 988-991.
00014  * The primary algorithm described there uses 11 multiplies and 29 adds.
00015  * We use their alternate method with 12 multiplies and 32 adds.
00016  * The advantage of this method is that no data path contains more than one
00017  * multiplication; this allows a very simple and accurate implementation in
00018  * scaled fixed-point arithmetic, with a minimal number of shifts.
00019  * 
00020  * I've made lots of modifications to attempt to take advantage of the
00021  * sparse nature of the DCT matrices we're getting.  Although the logic
00022  * is cumbersome, it's straightforward and the resulting code is much
00023  * faster.
00024  *
00025  * A better way to do this would be to pass in the DCT block as a sparse
00026  * matrix, perhaps with the difference cases encoded.
00027  */
00028 
00029 #include <memory.h>
00030 #include "all.h"
00031 #include "ansi.h"
00032 #include "dct.h"
00033 
00034 
00035 #define CONST_BITS 13
00036 
00037 /*
00038  * This routine is specialized to the case DCTSIZE = 8.
00039  */
00040 
00041 #if DCTSIZE != 8
00042   Sorry, this code only copes with 8x8 DCTs. /* deliberate syntax err */
00043 #endif
00044 
00045 
00046 /*
00047  * A 2-D IDCT can be done by 1-D IDCT on each row followed by 1-D IDCT
00048  * on each column.  Direct algorithms are also available, but they are
00049  * much more complex and seem not to be any faster when reduced to code.
00050  *
00051  * The poop on this scaling stuff is as follows:
00052  *
00053  * Each 1-D IDCT step produces outputs which are a factor of sqrt(N)
00054  * larger than the true IDCT outputs.  The final outputs are therefore
00055  * a factor of N larger than desired; since N=8 this can be cured by
00056  * a simple right shift at the end of the algorithm.  The advantage of
00057  * this arrangement is that we save two multiplications per 1-D IDCT,
00058  * because the y0 and y4 inputs need not be divided by sqrt(N).
00059  *
00060  * We have to do addition and subtraction of the integer inputs, which
00061  * is no problem, and multiplication by fractional constants, which is
00062  * a problem to do in integer arithmetic.  We multiply all the constants
00063  * by CONST_SCALE and convert them to integer constants (thus retaining
00064  * CONST_BITS bits of precision in the constants).  After doing a
00065  * multiplication we have to divide the product by CONST_SCALE, with proper
00066  * rounding, to produce the correct output.  This division can be done
00067  * cheaply as a right shift of CONST_BITS bits.  We postpone shifting
00068  * as long as possible so that partial sums can be added together with
00069  * full fractional precision.
00070  *
00071  * The outputs of the first pass are scaled up by PASS1_BITS bits so that
00072  * they are represented to better-than-integral precision.  These outputs
00073  * require BITS_IN_JSAMPLE + PASS1_BITS + 3 bits; this fits in a 16-bit word
00074  * with the recommended scaling.  (To scale up 12-bit sample data further, an
00075  * intermediate int32 array would be needed.)
00076  *
00077  * To avoid overflow of the 32-bit intermediate results in pass 2, we must
00078  * have BITS_IN_JSAMPLE + CONST_BITS + PASS1_BITS <= 26.  Error analysis
00079  * shows that the values given below are the most effective.
00080  */
00081 
00082 #ifdef EIGHT_BIT_SAMPLES
00083 #define PASS1_BITS  2
00084 #else
00085 #define PASS1_BITS  1           /* lose a little precision to avoid overflow */
00086 #endif
00087 
00088 #define ONE     ((int32) 1)
00089 
00090 #define CONST_SCALE (ONE << CONST_BITS)
00091 
00092 /* Convert a positive real constant to an integer scaled by CONST_SCALE.
00093  * IMPORTANT: if your compiler doesn't do this arithmetic at compile time,
00094  * you will pay a significant penalty in run time.  In that case, figure
00095  * the correct integer constant values and insert them by hand.
00096  */
00097 
00098 /* Actually FIX is no longer used, we precomputed them all */
00099 #define FIX(x)  ((int32) ((x) * CONST_SCALE + 0.5)) 
00100 
00101 /* Descale and correctly round an int32 value that's scaled by N bits.
00102  * We assume RIGHT_SHIFT rounds towards minus infinity, so adding
00103  * the fudge factor is correct for either sign of X.
00104  */
00105 
00106 #define DESCALE(x,n)  RIGHT_SHIFT((x) + (ONE << ((n)-1)), n)
00107 
00108 /* Multiply an int32 variable by an int32 constant to yield an int32 result.
00109  * For 8-bit samples with the recommended scaling, all the variable
00110  * and constant values involved are no more than 16 bits wide, so a
00111  * 16x16->32 bit multiply can be used instead of a full 32x32 multiply;
00112  * this provides a useful speedup on many machines.
00113  * There is no way to specify a 16x16->32 multiply in portable C, but
00114  * some C compilers will do the right thing if you provide the correct
00115  * combination of casts.
00116  * NB: for 12-bit samples, a full 32-bit multiplication will be needed.
00117  */
00118 
00119 #ifdef EIGHT_BIT_SAMPLES
00120 #ifdef SHORTxSHORT_32           /* may work if 'int' is 32 bits */
00121 #define MULTIPLY(var,const)  (((INT16) (var)) * ((INT16) (const)))
00122 #endif
00123 #ifdef SHORTxLCONST_32          /* known to work with Microsoft C 6.0 */
00124 #define MULTIPLY(var,const)  (((INT16) (var)) * ((int32) (const)))
00125 #endif
00126 #endif
00127 
00128 #ifndef MULTIPLY                /* default definition */
00129 #define MULTIPLY(var,const)  ((var) * (const))
00130 #endif
00131 
00132 
00133 /* 
00134   Unlike our decoder where we approximate the FIXes, we need to use exact
00135 ones here or successive P-frames will drift too much with Reference frame coding 
00136 */
00137 #define FIX_0_211164243 1730
00138 #define FIX_0_275899380 2260
00139 #define FIX_0_298631336 2446
00140 #define FIX_0_390180644 3196
00141 #define FIX_0_509795579 4176
00142 #define FIX_0_541196100 4433
00143 #define FIX_0_601344887 4926
00144 #define FIX_0_765366865 6270
00145 #define FIX_0_785694958 6436
00146 #define FIX_0_899976223 7373
00147 #define FIX_1_061594337 8697
00148 #define FIX_1_111140466 9102
00149 #define FIX_1_175875602 9633
00150 #define FIX_1_306562965 10703
00151 #define FIX_1_387039845 11363
00152 #define FIX_1_451774981 11893
00153 #define FIX_1_501321110 12299
00154 #define FIX_1_662939225 13623
00155 #define FIX_1_847759065 15137
00156 #define FIX_1_961570560 16069
00157 #define FIX_2_053119869 16819
00158 #define FIX_2_172734803 17799
00159 #define FIX_2_562915447 20995
00160 #define FIX_3_072711026 25172
00161 
00162 /*
00163   Switch on reverse_dct choices
00164 */
00165 void reference_rev_dct _ANSI_ARGS_((int16 *block));
00166 void mpeg_jrevdct_quick _ANSI_ARGS_((int16 *block));
00167 void init_idctref _ANSI_ARGS_((void));
00168 
00169 extern boolean pureDCT;
00170 
00171 void
00172 mpeg_jrevdct(data)
00173     DCTBLOCK data;
00174 {
00175   if (pureDCT) reference_rev_dct(data);
00176   else mpeg_jrevdct_quick(data);
00177 }
00178 
00179 /*
00180  * Perform the inverse DCT on one block of coefficients.
00181  */
00182 
00183 void
00184 mpeg_jrevdct_quick(data)
00185     DCTBLOCK data;
00186 {
00187   int32 tmp0, tmp1, tmp2, tmp3;
00188   int32 tmp10, tmp11, tmp12, tmp13;
00189   int32 z1, z2, z3, z4, z5;
00190   int32 d0, d1, d2, d3, d4, d5, d6, d7;
00191   register DCTELEM *dataptr;
00192   int rowctr;
00193   SHIFT_TEMPS
00194    
00195   /* Pass 1: process rows. */
00196   /* Note results are scaled up by sqrt(8) compared to a true IDCT; */
00197   /* furthermore, we scale the results by 2**PASS1_BITS. */
00198 
00199   dataptr = data;
00200 
00201   for (rowctr = DCTSIZE-1; rowctr >= 0; rowctr--) {
00202     /* Due to quantization, we will usually find that many of the input
00203      * coefficients are zero, especially the AC terms.  We can exploit this
00204      * by short-circuiting the IDCT calculation for any row in which all
00205      * the AC terms are zero.  In that case each output is equal to the
00206      * DC coefficient (with scale factor as needed).
00207      * With typical images and quantization tables, half or more of the
00208      * row DCT calculations can be simplified this way.
00209      */
00210 
00211     register int *idataptr = (int*)dataptr;
00212     d0 = dataptr[0];
00213     d1 = dataptr[1];
00214     if ((d1 == 0) && (idataptr[1] | idataptr[2] | idataptr[3]) == 0) {
00215       /* AC terms all zero */
00216       if (d0) {
00217           /* Compute a 32 bit value to assign. */
00218           DCTELEM dcval = (DCTELEM) (d0 << PASS1_BITS);
00219           register int v = (dcval & 0xffff) | ((dcval << 16) & 0xffff0000);
00220           
00221           idataptr[0] = v;
00222           idataptr[1] = v;
00223           idataptr[2] = v;
00224           idataptr[3] = v;
00225       }
00226       
00227       dataptr += DCTSIZE;       /* advance pointer to next row */
00228       continue;
00229     }
00230     d2 = dataptr[2];
00231     d3 = dataptr[3];
00232     d4 = dataptr[4];
00233     d5 = dataptr[5];
00234     d6 = dataptr[6];
00235     d7 = dataptr[7];
00236 
00237     /* Even part: reverse the even part of the forward DCT. */
00238     /* The rotator is sqrt(2)*c(-6). */
00239 {
00240     if (d6) {
00241         if (d4) {
00242             if (d2) {
00243                 if (d0) {
00244                     /* d0 != 0, d2 != 0, d4 != 0, d6 != 0 */
00245                     z1 = MULTIPLY(d2 + d6, FIX_0_541196100);
00246                     tmp2 = z1 + MULTIPLY(-d6, FIX_1_847759065);
00247                     tmp3 = z1 + MULTIPLY(d2, FIX_0_765366865);
00248 
00249                     tmp0 = (d0 + d4) << CONST_BITS;
00250                     tmp1 = (d0 - d4) << CONST_BITS;
00251 
00252                     tmp10 = tmp0 + tmp3;
00253                     tmp13 = tmp0 - tmp3;
00254                     tmp11 = tmp1 + tmp2;
00255                     tmp12 = tmp1 - tmp2;
00256                 } else {
00257                     /* d0 == 0, d2 != 0, d4 != 0, d6 != 0 */
00258                     z1 = MULTIPLY(d2 + d6, FIX_0_541196100);
00259                     tmp2 = z1 + MULTIPLY(-d6, FIX_1_847759065);
00260                     tmp3 = z1 + MULTIPLY(d2, FIX_0_765366865);
00261 
00262                     tmp0 = d4 << CONST_BITS;
00263 
00264                     tmp10 = tmp0 + tmp3;
00265                     tmp13 = tmp0 - tmp3;
00266                     tmp11 = tmp2 - tmp0;
00267                     tmp12 = -(tmp0 + tmp2);
00268                 }
00269             } else {
00270                 if (d0) {
00271                     /* d0 != 0, d2 == 0, d4 != 0, d6 != 0 */
00272                     tmp2 = MULTIPLY(-d6, FIX_1_306562965);
00273                     tmp3 = MULTIPLY(d6, FIX_0_541196100);
00274 
00275                     tmp0 = (d0 + d4) << CONST_BITS;
00276                     tmp1 = (d0 - d4) << CONST_BITS;
00277 
00278                     tmp10 = tmp0 + tmp3;
00279                     tmp13 = tmp0 - tmp3;
00280                     tmp11 = tmp1 + tmp2;
00281                     tmp12 = tmp1 - tmp2;
00282                 } else {
00283                     /* d0 == 0, d2 == 0, d4 != 0, d6 != 0 */
00284                     tmp2 = MULTIPLY(-d6, FIX_1_306562965);
00285                     tmp3 = MULTIPLY(d6, FIX_0_541196100);
00286 
00287                     tmp0 = d4 << CONST_BITS;
00288 
00289                     tmp10 = tmp0 + tmp3;
00290                     tmp13 = tmp0 - tmp3;
00291                     tmp11 = tmp2 - tmp0;
00292                     tmp12 = -(tmp0 + tmp2);
00293                 }
00294             }
00295         } else {
00296             if (d2) {
00297                 if (d0) {
00298                     /* d0 != 0, d2 != 0, d4 == 0, d6 != 0 */
00299                     z1 = MULTIPLY(d2 + d6, FIX_0_541196100);
00300                     tmp2 = z1 + MULTIPLY(-d6, FIX_1_847759065);
00301                     tmp3 = z1 + MULTIPLY(d2, FIX_0_765366865);
00302 
00303                     tmp0 = d0 << CONST_BITS;
00304 
00305                     tmp10 = tmp0 + tmp3;
00306                     tmp13 = tmp0 - tmp3;
00307                     tmp11 = tmp0 + tmp2;
00308                     tmp12 = tmp0 - tmp2;
00309                 } else {
00310                     /* d0 == 0, d2 != 0, d4 == 0, d6 != 0 */
00311                     z1 = MULTIPLY(d2 + d6, FIX_0_541196100);
00312                     tmp2 = z1 + MULTIPLY(-d6, FIX_1_847759065);
00313                     tmp3 = z1 + MULTIPLY(d2, FIX_0_765366865);
00314 
00315                     tmp10 = tmp3;
00316                     tmp13 = -tmp3;
00317                     tmp11 = tmp2;
00318                     tmp12 = -tmp2;
00319                 }
00320             } else {
00321                 if (d0) {
00322                     /* d0 != 0, d2 == 0, d4 == 0, d6 != 0 */
00323                     tmp2 = MULTIPLY(-d6, FIX_1_306562965);
00324                     tmp3 = MULTIPLY(d6, FIX_0_541196100);
00325 
00326                     tmp0 = d0 << CONST_BITS;
00327 
00328                     tmp10 = tmp0 + tmp3;
00329                     tmp13 = tmp0 - tmp3;
00330                     tmp11 = tmp0 + tmp2;
00331                     tmp12 = tmp0 - tmp2;
00332                 } else {
00333                     /* d0 == 0, d2 == 0, d4 == 0, d6 != 0 */
00334                     tmp2 = MULTIPLY(-d6, FIX_1_306562965);
00335                     tmp3 = MULTIPLY(d6, FIX_0_541196100);
00336 
00337                     tmp10 = tmp3;
00338                     tmp13 = -tmp3;
00339                     tmp11 = tmp2;
00340                     tmp12 = -tmp2;
00341                 }
00342             }
00343         }
00344     } else {
00345         if (d4) {
00346             if (d2) {
00347                 if (d0) {
00348                     /* d0 != 0, d2 != 0, d4 != 0, d6 == 0 */
00349                     tmp2 = MULTIPLY(d2, FIX_0_541196100);
00350                     tmp3 = MULTIPLY(d2, FIX_1_306562965);
00351 
00352                     tmp0 = (d0 + d4) << CONST_BITS;
00353                     tmp1 = (d0 - d4) << CONST_BITS;
00354 
00355                     tmp10 = tmp0 + tmp3;
00356                     tmp13 = tmp0 - tmp3;
00357                     tmp11 = tmp1 + tmp2;
00358                     tmp12 = tmp1 - tmp2;
00359                 } else {
00360                     /* d0 == 0, d2 != 0, d4 != 0, d6 == 0 */
00361                     tmp2 = MULTIPLY(d2, FIX_0_541196100);
00362                     tmp3 = MULTIPLY(d2, FIX_1_306562965);
00363 
00364                     tmp0 = d4 << CONST_BITS;
00365 
00366                     tmp10 = tmp0 + tmp3;
00367                     tmp13 = tmp0 - tmp3;
00368                     tmp11 = tmp2 - tmp0;
00369                     tmp12 = -(tmp0 + tmp2);
00370                 }
00371             } else {
00372                 if (d0) {
00373                     /* d0 != 0, d2 == 0, d4 != 0, d6 == 0 */
00374                     tmp10 = tmp13 = (d0 + d4) << CONST_BITS;
00375                     tmp11 = tmp12 = (d0 - d4) << CONST_BITS;
00376                 } else {
00377                     /* d0 == 0, d2 == 0, d4 != 0, d6 == 0 */
00378                     tmp10 = tmp13 = d4 << CONST_BITS;
00379                     tmp11 = tmp12 = -tmp10;
00380                 }
00381             }
00382         } else {
00383             if (d2) {
00384                 if (d0) {
00385                     /* d0 != 0, d2 != 0, d4 == 0, d6 == 0 */
00386                     tmp2 = MULTIPLY(d2, FIX_0_541196100);
00387                     tmp3 = MULTIPLY(d2, FIX_1_306562965);
00388 
00389                     tmp0 = d0 << CONST_BITS;
00390 
00391                     tmp10 = tmp0 + tmp3;
00392                     tmp13 = tmp0 - tmp3;
00393                     tmp11 = tmp0 + tmp2;
00394                     tmp12 = tmp0 - tmp2;
00395                 } else {
00396                     /* d0 == 0, d2 != 0, d4 == 0, d6 == 0 */
00397                     tmp2 = MULTIPLY(d2, FIX_0_541196100);
00398                     tmp3 = MULTIPLY(d2, FIX_1_306562965);
00399 
00400                     tmp10 = tmp3;
00401                     tmp13 = -tmp3;
00402                     tmp11 = tmp2;
00403                     tmp12 = -tmp2;
00404                 }
00405             } else {
00406                 if (d0) {
00407                     /* d0 != 0, d2 == 0, d4 == 0, d6 == 0 */
00408                     tmp10 = tmp13 = tmp11 = tmp12 = d0 << CONST_BITS;
00409                 } else {
00410                     /* d0 == 0, d2 == 0, d4 == 0, d6 == 0 */
00411                     tmp10 = tmp13 = tmp11 = tmp12 = 0;
00412                 }
00413             }
00414         }
00415       }
00416 
00417     /* Odd part per figure 8; the matrix is unitary and hence its
00418      * transpose is its inverse.  i0..i3 are y7,y5,y3,y1 respectively.
00419      */
00420 
00421     if (d7) {
00422         if (d5) {
00423             if (d3) {
00424                 if (d1) {
00425                     /* d1 != 0, d3 != 0, d5 != 0, d7 != 0 */
00426                     z1 = d7 + d1;
00427                     z2 = d5 + d3;
00428                     z3 = d7 + d3;
00429                     z4 = d5 + d1;
00430                     z5 = MULTIPLY(z3 + z4, FIX_1_175875602);
00431                     
00432                     tmp0 = MULTIPLY(d7, FIX_0_298631336); 
00433                     tmp1 = MULTIPLY(d5, FIX_2_053119869);
00434                     tmp2 = MULTIPLY(d3, FIX_3_072711026);
00435                     tmp3 = MULTIPLY(d1, FIX_1_501321110);
00436                     z1 = MULTIPLY(-z1, FIX_0_899976223);
00437                     z2 = MULTIPLY(-z2, FIX_2_562915447);
00438                     z3 = MULTIPLY(-z3, FIX_1_961570560);
00439                     z4 = MULTIPLY(-z4, FIX_0_390180644);
00440                     
00441                     z3 += z5;
00442                     z4 += z5;
00443                     
00444                     tmp0 += z1 + z3;
00445                     tmp1 += z2 + z4;
00446                     tmp2 += z2 + z3;
00447                     tmp3 += z1 + z4;
00448                 } else {
00449                     /* d1 == 0, d3 != 0, d5 != 0, d7 != 0 */
00450                     z2 = d5 + d3;
00451                     z3 = d7 + d3;
00452                     z5 = MULTIPLY(z3 + d5, FIX_1_175875602);
00453                     
00454                     tmp0 = MULTIPLY(d7, FIX_0_298631336); 
00455                     tmp1 = MULTIPLY(d5, FIX_2_053119869);
00456                     tmp2 = MULTIPLY(d3, FIX_3_072711026);
00457                     z1 = MULTIPLY(-d7, FIX_0_899976223);
00458                     z2 = MULTIPLY(-z2, FIX_2_562915447);
00459                     z3 = MULTIPLY(-z3, FIX_1_961570560);
00460                     z4 = MULTIPLY(-d5, FIX_0_390180644);
00461                     
00462                     z3 += z5;
00463                     z4 += z5;
00464                     
00465                     tmp0 += z1 + z3;
00466                     tmp1 += z2 + z4;
00467                     tmp2 += z2 + z3;
00468                     tmp3 = z1 + z4;
00469                 }
00470             } else {
00471                 if (d1) {
00472                     /* d1 != 0, d3 == 0, d5 != 0, d7 != 0 */
00473                     z1 = d7 + d1;
00474                     z4 = d5 + d1;
00475                     z5 = MULTIPLY(d7 + z4, FIX_1_175875602);
00476                     
00477                     tmp0 = MULTIPLY(d7, FIX_0_298631336); 
00478                     tmp1 = MULTIPLY(d5, FIX_2_053119869);
00479                     tmp3 = MULTIPLY(d1, FIX_1_501321110);
00480                     z1 = MULTIPLY(-z1, FIX_0_899976223);
00481                     z2 = MULTIPLY(-d5, FIX_2_562915447);
00482                     z3 = MULTIPLY(-d7, FIX_1_961570560);
00483                     z4 = MULTIPLY(-z4, FIX_0_390180644);
00484                     
00485                     z3 += z5;
00486                     z4 += z5;
00487                     
00488                     tmp0 += z1 + z3;
00489                     tmp1 += z2 + z4;
00490                     tmp2 = z2 + z3;
00491                     tmp3 += z1 + z4;
00492                 } else {
00493                     /* d1 == 0, d3 == 0, d5 != 0, d7 != 0 */
00494                     tmp0 = MULTIPLY(-d7, FIX_0_601344887); 
00495                     z1 = MULTIPLY(-d7, FIX_0_899976223);
00496                     z3 = MULTIPLY(-d7, FIX_1_961570560);
00497                     tmp1 = MULTIPLY(-d5, FIX_0_509795579);
00498                     z2 = MULTIPLY(-d5, FIX_2_562915447);
00499                     z4 = MULTIPLY(-d5, FIX_0_390180644);
00500                     z5 = MULTIPLY(d5 + d7, FIX_1_175875602);
00501                     
00502                     z3 += z5;
00503                     z4 += z5;
00504                     
00505                     tmp0 += z3;
00506                     tmp1 += z4;
00507                     tmp2 = z2 + z3;
00508                     tmp3 = z1 + z4;
00509                 }
00510             }
00511         } else {
00512             if (d3) {
00513                 if (d1) {
00514                     /* d1 != 0, d3 != 0, d5 == 0, d7 != 0 */
00515                     z1 = d7 + d1;
00516                     z3 = d7 + d3;
00517                     z5 = MULTIPLY(z3 + d1, FIX_1_175875602);
00518                     
00519                     tmp0 = MULTIPLY(d7, FIX_0_298631336); 
00520                     tmp2 = MULTIPLY(d3, FIX_3_072711026);
00521                     tmp3 = MULTIPLY(d1, FIX_1_501321110);
00522                     z1 = MULTIPLY(-z1, FIX_0_899976223);
00523                     z2 = MULTIPLY(-d3, FIX_2_562915447);
00524                     z3 = MULTIPLY(-z3, FIX_1_961570560);
00525                     z4 = MULTIPLY(-d1, FIX_0_390180644);
00526                     
00527                     z3 += z5;
00528                     z4 += z5;
00529                     
00530                     tmp0 += z1 + z3;
00531                     tmp1 = z2 + z4;
00532                     tmp2 += z2 + z3;
00533                     tmp3 += z1 + z4;
00534                 } else {
00535                     /* d1 == 0, d3 != 0, d5 == 0, d7 != 0 */
00536                     z3 = d7 + d3;
00537                     
00538                     tmp0 = MULTIPLY(-d7, FIX_0_601344887); 
00539                     z1 = MULTIPLY(-d7, FIX_0_899976223);
00540                     tmp2 = MULTIPLY(d3, FIX_0_509795579);
00541                     z2 = MULTIPLY(-d3, FIX_2_562915447);
00542                     z5 = MULTIPLY(z3, FIX_1_175875602);
00543                     z3 = MULTIPLY(-z3, FIX_0_785694958);
00544                     
00545                     tmp0 += z3;
00546                     tmp1 = z2 + z5;
00547                     tmp2 += z3;
00548                     tmp3 = z1 + z5;
00549                 }
00550             } else {
00551                 if (d1) {
00552                     /* d1 != 0, d3 == 0, d5 == 0, d7 != 0 */
00553                     z1 = d7 + d1;
00554                     z5 = MULTIPLY(z1, FIX_1_175875602);
00555 
00556                     z1 = MULTIPLY(z1, FIX_0_275899380);
00557                     z3 = MULTIPLY(-d7, FIX_1_961570560);
00558                     tmp0 = MULTIPLY(-d7, FIX_1_662939225); 
00559                     z4 = MULTIPLY(-d1, FIX_0_390180644);
00560                     tmp3 = MULTIPLY(d1, FIX_1_111140466);
00561 
00562                     tmp0 += z1;
00563                     tmp1 = z4 + z5;
00564                     tmp2 = z3 + z5;
00565                     tmp3 += z1;
00566                 } else {
00567                     /* d1 == 0, d3 == 0, d5 == 0, d7 != 0 */
00568                     tmp0 = MULTIPLY(-d7, FIX_1_387039845);
00569                     tmp1 = MULTIPLY(d7, FIX_1_175875602);
00570                     tmp2 = MULTIPLY(-d7, FIX_0_785694958);
00571                     tmp3 = MULTIPLY(d7, FIX_0_275899380);
00572                 }
00573             }
00574         }
00575     } else {
00576         if (d5) {
00577             if (d3) {
00578                 if (d1) {
00579                     /* d1 != 0, d3 != 0, d5 != 0, d7 == 0 */
00580                     z2 = d5 + d3;
00581                     z4 = d5 + d1;
00582                     z5 = MULTIPLY(d3 + z4, FIX_1_175875602);
00583                     
00584                     tmp1 = MULTIPLY(d5, FIX_2_053119869);
00585                     tmp2 = MULTIPLY(d3, FIX_3_072711026);
00586                     tmp3 = MULTIPLY(d1, FIX_1_501321110);
00587                     z1 = MULTIPLY(-d1, FIX_0_899976223);
00588                     z2 = MULTIPLY(-z2, FIX_2_562915447);
00589                     z3 = MULTIPLY(-d3, FIX_1_961570560);
00590                     z4 = MULTIPLY(-z4, FIX_0_390180644);
00591                     
00592                     z3 += z5;
00593                     z4 += z5;
00594                     
00595                     tmp0 = z1 + z3;
00596                     tmp1 += z2 + z4;
00597                     tmp2 += z2 + z3;
00598                     tmp3 += z1 + z4;
00599                 } else {
00600                     /* d1 == 0, d3 != 0, d5 != 0, d7 == 0 */
00601                     z2 = d5 + d3;
00602                     
00603                     z5 = MULTIPLY(z2, FIX_1_175875602);
00604                     tmp1 = MULTIPLY(d5, FIX_1_662939225);
00605                     z4 = MULTIPLY(-d5, FIX_0_390180644);
00606                     z2 = MULTIPLY(-z2, FIX_1_387039845);
00607                     tmp2 = MULTIPLY(d3, FIX_1_111140466);
00608                     z3 = MULTIPLY(-d3, FIX_1_961570560);
00609                     
00610                     tmp0 = z3 + z5;
00611                     tmp1 += z2;
00612                     tmp2 += z2;
00613                     tmp3 = z4 + z5;
00614                 }
00615             } else {
00616                 if (d1) {
00617                     /* d1 != 0, d3 == 0, d5 != 0, d7 == 0 */
00618                     z4 = d5 + d1;
00619                     
00620                     z5 = MULTIPLY(z4, FIX_1_175875602);
00621                     z1 = MULTIPLY(-d1, FIX_0_899976223);
00622                     tmp3 = MULTIPLY(d1, FIX_0_601344887);
00623                     tmp1 = MULTIPLY(-d5, FIX_0_509795579);
00624                     z2 = MULTIPLY(-d5, FIX_2_562915447);
00625                     z4 = MULTIPLY(z4, FIX_0_785694958);
00626                     
00627                     tmp0 = z1 + z5;
00628                     tmp1 += z4;
00629                     tmp2 = z2 + z5;
00630                     tmp3 += z4;
00631                 } else {
00632                     /* d1 == 0, d3 == 0, d5 != 0, d7 == 0 */
00633                     tmp0 = MULTIPLY(d5, FIX_1_175875602);
00634                     tmp1 = MULTIPLY(d5, FIX_0_275899380);
00635                     tmp2 = MULTIPLY(-d5, FIX_1_387039845);
00636                     tmp3 = MULTIPLY(d5, FIX_0_785694958);
00637                 }
00638             }
00639         } else {
00640             if (d3) {
00641                 if (d1) {
00642                     /* d1 != 0, d3 != 0, d5 == 0, d7 == 0 */
00643                     z5 = d1 + d3;
00644                     tmp3 = MULTIPLY(d1, FIX_0_211164243);
00645                     tmp2 = MULTIPLY(-d3, FIX_1_451774981);
00646                     z1 = MULTIPLY(d1, FIX_1_061594337);
00647                     z2 = MULTIPLY(-d3, FIX_2_172734803);
00648                     z4 = MULTIPLY(z5, FIX_0_785694958);
00649                     z5 = MULTIPLY(z5, FIX_1_175875602);
00650                     
00651                     tmp0 = z1 - z4;
00652                     tmp1 = z2 + z4;
00653                     tmp2 += z5;
00654                     tmp3 += z5;
00655                 } else {
00656                     /* d1 == 0, d3 != 0, d5 == 0, d7 == 0 */
00657                     tmp0 = MULTIPLY(-d3, FIX_0_785694958);
00658                     tmp1 = MULTIPLY(-d3, FIX_1_387039845);
00659                     tmp2 = MULTIPLY(-d3, FIX_0_275899380);
00660                     tmp3 = MULTIPLY(d3, FIX_1_175875602);
00661                 }
00662             } else {
00663                 if (d1) {
00664                     /* d1 != 0, d3 == 0, d5 == 0, d7 == 0 */
00665                     tmp0 = MULTIPLY(d1, FIX_0_275899380);
00666                     tmp1 = MULTIPLY(d1, FIX_0_785694958);
00667                     tmp2 = MULTIPLY(d1, FIX_1_175875602);
00668                     tmp3 = MULTIPLY(d1, FIX_1_387039845);
00669                 } else {
00670                     /* d1 == 0, d3 == 0, d5 == 0, d7 == 0 */
00671                     tmp0 = tmp1 = tmp2 = tmp3 = 0;
00672                 }
00673             }
00674         }
00675     }
00676 }
00677     /* Final output stage: inputs are tmp10..tmp13, tmp0..tmp3 */
00678 
00679     dataptr[0] = (DCTELEM) DESCALE(tmp10 + tmp3, CONST_BITS-PASS1_BITS);
00680     dataptr[7] = (DCTELEM) DESCALE(tmp10 - tmp3, CONST_BITS-PASS1_BITS);
00681     dataptr[1] = (DCTELEM) DESCALE(tmp11 + tmp2, CONST_BITS-PASS1_BITS);
00682     dataptr[6] = (DCTELEM) DESCALE(tmp11 - tmp2, CONST_BITS-PASS1_BITS);
00683     dataptr[2] = (DCTELEM) DESCALE(tmp12 + tmp1, CONST_BITS-PASS1_BITS);
00684     dataptr[5] = (DCTELEM) DESCALE(tmp12 - tmp1, CONST_BITS-PASS1_BITS);
00685     dataptr[3] = (DCTELEM) DESCALE(tmp13 + tmp0, CONST_BITS-PASS1_BITS);
00686     dataptr[4] = (DCTELEM) DESCALE(tmp13 - tmp0, CONST_BITS-PASS1_BITS);
00687 
00688     dataptr += DCTSIZE;         /* advance pointer to next row */
00689   }
00690 
00691   /* Pass 2: process columns. */
00692   /* Note that we must descale the results by a factor of 8 == 2**3, */
00693   /* and also undo the PASS1_BITS scaling. */
00694 
00695   dataptr = data;
00696   for (rowctr = DCTSIZE-1; rowctr >= 0; rowctr--) {
00697     /* Columns of zeroes can be exploited in the same way as we did with rows.
00698      * However, the row calculation has created many nonzero AC terms, so the
00699      * simplification applies less often (typically 5% to 10% of the time).
00700      * On machines with very fast multiplication, it's possible that the
00701      * test takes more time than it's worth.  In that case this section
00702      * may be commented out.
00703      */
00704 
00705     d0 = dataptr[DCTSIZE*0];
00706     d1 = dataptr[DCTSIZE*1];
00707     d2 = dataptr[DCTSIZE*2];
00708     d3 = dataptr[DCTSIZE*3];
00709     d4 = dataptr[DCTSIZE*4];
00710     d5 = dataptr[DCTSIZE*5];
00711     d6 = dataptr[DCTSIZE*6];
00712     d7 = dataptr[DCTSIZE*7];
00713 
00714     /* Even part: reverse the even part of the forward DCT. */
00715     /* The rotator is sqrt(2)*c(-6). */
00716     if (d6) {
00717         if (d4) {
00718             if (d2) {
00719                 if (d0) {
00720                     /* d0 != 0, d2 != 0, d4 != 0, d6 != 0 */
00721                     z1 = MULTIPLY(d2 + d6, FIX_0_541196100);
00722                     tmp2 = z1 + MULTIPLY(-d6, FIX_1_847759065);
00723                     tmp3 = z1 + MULTIPLY(d2, FIX_0_765366865);
00724 
00725                     tmp0 = (d0 + d4) << CONST_BITS;
00726                     tmp1 = (d0 - d4) << CONST_BITS;
00727 
00728                     tmp10 = tmp0 + tmp3;
00729                     tmp13 = tmp0 - tmp3;
00730                     tmp11 = tmp1 + tmp2;
00731                     tmp12 = tmp1 - tmp2;
00732                 } else {
00733                     /* d0 == 0, d2 != 0, d4 != 0, d6 != 0 */
00734                     z1 = MULTIPLY(d2 + d6, FIX_0_541196100);
00735                     tmp2 = z1 + MULTIPLY(-d6, FIX_1_847759065);
00736                     tmp3 = z1 + MULTIPLY(d2, FIX_0_765366865);
00737 
00738                     tmp0 = d4 << CONST_BITS;
00739 
00740                     tmp10 = tmp0 + tmp3;
00741                     tmp13 = tmp0 - tmp3;
00742                     tmp11 = tmp2 - tmp0;
00743                     tmp12 = -(tmp0 + tmp2);
00744                 }
00745             } else {
00746                 if (d0) {
00747                     /* d0 != 0, d2 == 0, d4 != 0, d6 != 0 */
00748                     tmp2 = MULTIPLY(-d6, FIX_1_306562965);
00749                     tmp3 = MULTIPLY(d6, FIX_0_541196100);
00750 
00751                     tmp0 = (d0 + d4) << CONST_BITS;
00752                     tmp1 = (d0 - d4) << CONST_BITS;
00753 
00754                     tmp10 = tmp0 + tmp3;
00755                     tmp13 = tmp0 - tmp3;
00756                     tmp11 = tmp1 + tmp2;
00757                     tmp12 = tmp1 - tmp2;
00758                 } else {
00759                     /* d0 == 0, d2 == 0, d4 != 0, d6 != 0 */
00760                     tmp2 = MULTIPLY(-d6, FIX_1_306562965);
00761                     tmp3 = MULTIPLY(d6, FIX_0_541196100);
00762 
00763                     tmp0 = d4 << CONST_BITS;
00764 
00765                     tmp10 = tmp0 + tmp3;
00766                     tmp13 = tmp0 - tmp3;
00767                     tmp11 = tmp2 - tmp0;
00768                     tmp12 = -(tmp0 + tmp2);
00769                 }
00770             }
00771         } else {
00772             if (d2) {
00773                 if (d0) {
00774                     /* d0 != 0, d2 != 0, d4 == 0, d6 != 0 */
00775                     z1 = MULTIPLY(d2 + d6, FIX_0_541196100);
00776                     tmp2 = z1 + MULTIPLY(-d6, FIX_1_847759065);
00777                     tmp3 = z1 + MULTIPLY(d2, FIX_0_765366865);
00778 
00779                     tmp0 = d0 << CONST_BITS;
00780 
00781                     tmp10 = tmp0 + tmp3;
00782                     tmp13 = tmp0 - tmp3;
00783                     tmp11 = tmp0 + tmp2;
00784                     tmp12 = tmp0 - tmp2;
00785                 } else {
00786                     /* d0 == 0, d2 != 0, d4 == 0, d6 != 0 */
00787                     z1 = MULTIPLY(d2 + d6, FIX_0_541196100);
00788                     tmp2 = z1 + MULTIPLY(-d6, FIX_1_847759065);
00789                     tmp3 = z1 + MULTIPLY(d2, FIX_0_765366865);
00790 
00791                     tmp10 = tmp3;
00792                     tmp13 = -tmp3;
00793                     tmp11 = tmp2;
00794                     tmp12 = -tmp2;
00795                 }
00796             } else {
00797                 if (d0) {
00798                     /* d0 != 0, d2 == 0, d4 == 0, d6 != 0 */
00799                     tmp2 = MULTIPLY(-d6, FIX_1_306562965);
00800                     tmp3 = MULTIPLY(d6, FIX_0_541196100);
00801 
00802                     tmp0 = d0 << CONST_BITS;
00803 
00804                     tmp10 = tmp0 + tmp3;
00805                     tmp13 = tmp0 - tmp3;
00806                     tmp11 = tmp0 + tmp2;
00807                     tmp12 = tmp0 - tmp2;
00808                 } else {
00809                     /* d0 == 0, d2 == 0, d4 == 0, d6 != 0 */
00810                     tmp2 = MULTIPLY(-d6, FIX_1_306562965);
00811                     tmp3 = MULTIPLY(d6, FIX_0_541196100);
00812 
00813                     tmp10 = tmp3;
00814                     tmp13 = -tmp3;
00815                     tmp11 = tmp2;
00816                     tmp12 = -tmp2;
00817                 }
00818             }
00819         }
00820     } else {
00821         if (d4) {
00822             if (d2) {
00823                 if (d0) {
00824                     /* d0 != 0, d2 != 0, d4 != 0, d6 == 0 */
00825                     tmp2 = MULTIPLY(d2, FIX_0_541196100);
00826                     tmp3 = MULTIPLY(d2, FIX_1_306562965);
00827 
00828                     tmp0 = (d0 + d4) << CONST_BITS;
00829                     tmp1 = (d0 - d4) << CONST_BITS;
00830 
00831                     tmp10 = tmp0 + tmp3;
00832                     tmp13 = tmp0 - tmp3;
00833                     tmp11 = tmp1 + tmp2;
00834                     tmp12 = tmp1 - tmp2;
00835                 } else {
00836                     /* d0 == 0, d2 != 0, d4 != 0, d6 == 0 */
00837                     tmp2 = MULTIPLY(d2, FIX_0_541196100);
00838                     tmp3 = MULTIPLY(d2, FIX_1_306562965);
00839 
00840                     tmp0 = d4 << CONST_BITS;
00841 
00842                     tmp10 = tmp0 + tmp3;
00843                     tmp13 = tmp0 - tmp3;
00844                     tmp11 = tmp2 - tmp0;
00845                     tmp12 = -(tmp0 + tmp2);
00846                 }
00847             } else {
00848                 if (d0) {
00849                     /* d0 != 0, d2 == 0, d4 != 0, d6 == 0 */
00850                     tmp10 = tmp13 = (d0 + d4) << CONST_BITS;
00851                     tmp11 = tmp12 = (d0 - d4) << CONST_BITS;
00852                 } else {
00853                     /* d0 == 0, d2 == 0, d4 != 0, d6 == 0 */
00854                     tmp10 = tmp13 = d4 << CONST_BITS;
00855                     tmp11 = tmp12 = -tmp10;
00856                 }
00857             }
00858         } else {
00859             if (d2) {
00860                 if (d0) {
00861                     /* d0 != 0, d2 != 0, d4 == 0, d6 == 0 */
00862                     tmp2 = MULTIPLY(d2, FIX_0_541196100);
00863                     tmp3 = MULTIPLY(d2, FIX_1_306562965);
00864 
00865                     tmp0 = d0 << CONST_BITS;
00866 
00867                     tmp10 = tmp0 + tmp3;
00868                     tmp13 = tmp0 - tmp3;
00869                     tmp11 = tmp0 + tmp2;
00870                     tmp12 = tmp0 - tmp2;
00871                 } else {
00872                     /* d0 == 0, d2 != 0, d4 == 0, d6 == 0 */
00873                     tmp2 = MULTIPLY(d2, FIX_0_541196100);
00874                     tmp3 = MULTIPLY(d2, FIX_1_306562965);
00875 
00876                     tmp10 = tmp3;
00877                     tmp13 = -tmp3;
00878                     tmp11 = tmp2;
00879                     tmp12 = -tmp2;
00880                 }
00881             } else {
00882                 if (d0) {
00883                     /* d0 != 0, d2 == 0, d4 == 0, d6 == 0 */
00884                     tmp10 = tmp13 = tmp11 = tmp12 = d0 << CONST_BITS;
00885                 } else {
00886                     /* d0 == 0, d2 == 0, d4 == 0, d6 == 0 */
00887                     tmp10 = tmp13 = tmp11 = tmp12 = 0;
00888                 }
00889             }
00890         }
00891     }
00892 
00893     /* Odd part per figure 8; the matrix is unitary and hence its
00894      * transpose is its inverse.  i0..i3 are y7,y5,y3,y1 respectively.
00895      */
00896     if (d7) {
00897         if (d5) {
00898             if (d3) {
00899                 if (d1) {
00900                     /* d1 != 0, d3 != 0, d5 != 0, d7 != 0 */
00901                     z1 = d7 + d1;
00902                     z2 = d5 + d3;
00903                     z3 = d7 + d3;
00904                     z4 = d5 + d1;
00905                     z5 = MULTIPLY(z3 + z4, FIX_1_175875602);
00906                     
00907                     tmp0 = MULTIPLY(d7, FIX_0_298631336); 
00908                     tmp1 = MULTIPLY(d5, FIX_2_053119869);
00909                     tmp2 = MULTIPLY(d3, FIX_3_072711026);
00910                     tmp3 = MULTIPLY(d1, FIX_1_501321110);
00911                     z1 = MULTIPLY(-z1, FIX_0_899976223);
00912                     z2 = MULTIPLY(-z2, FIX_2_562915447);
00913                     z3 = MULTIPLY(-z3, FIX_1_961570560);
00914                     z4 = MULTIPLY(-z4, FIX_0_390180644);
00915                     
00916                     z3 += z5;
00917                     z4 += z5;
00918                     
00919                     tmp0 += z1 + z3;
00920                     tmp1 += z2 + z4;
00921                     tmp2 += z2 + z3;
00922                     tmp3 += z1 + z4;
00923                 } else {
00924                     /* d1 == 0, d3 != 0, d5 != 0, d7 != 0 */
00925                     z1 = d7;
00926                     z2 = d5 + d3;
00927                     z3 = d7 + d3;
00928                     z5 = MULTIPLY(z3 + d5, FIX_1_175875602);
00929                     
00930                     tmp0 = MULTIPLY(d7, FIX_0_298631336); 
00931                     tmp1 = MULTIPLY(d5, FIX_2_053119869);
00932                     tmp2 = MULTIPLY(d3, FIX_3_072711026);
00933                     z1 = MULTIPLY(-d7, FIX_0_899976223);
00934                     z2 = MULTIPLY(-z2, FIX_2_562915447);
00935                     z3 = MULTIPLY(-z3, FIX_1_961570560);
00936                     z4 = MULTIPLY(-d5, FIX_0_390180644);
00937                     
00938                     z3 += z5;
00939                     z4 += z5;
00940                     
00941                     tmp0 += z1 + z3;
00942                     tmp1 += z2 + z4;
00943                     tmp2 += z2 + z3;
00944                     tmp3 = z1 + z4;
00945                 }
00946             } else {
00947                 if (d1) {
00948                     /* d1 != 0, d3 == 0, d5 != 0, d7 != 0 */
00949                     z1 = d7 + d1;
00950                     z2 = d5;
00951                     z3 = d7;
00952                     z4 = d5 + d1;
00953                     z5 = MULTIPLY(z3 + z4, FIX_1_175875602);
00954                     
00955                     tmp0 = MULTIPLY(d7, FIX_0_298631336); 
00956                     tmp1 = MULTIPLY(d5, FIX_2_053119869);
00957                     tmp3 = MULTIPLY(d1, FIX_1_501321110);
00958                     z1 = MULTIPLY(-z1, FIX_0_899976223);
00959                     z2 = MULTIPLY(-d5, FIX_2_562915447);
00960                     z3 = MULTIPLY(-d7, FIX_1_961570560);
00961                     z4 = MULTIPLY(-z4, FIX_0_390180644);
00962                     
00963                     z3 += z5;
00964                     z4 += z5;
00965                     
00966                     tmp0 += z1 + z3;
00967                     tmp1 += z2 + z4;
00968                     tmp2 = z2 + z3;
00969                     tmp3 += z1 + z4;
00970                 } else {
00971                     /* d1 == 0, d3 == 0, d5 != 0, d7 != 0 */
00972                     tmp0 = MULTIPLY(-d7, FIX_0_601344887); 
00973                     z1 = MULTIPLY(-d7, FIX_0_899976223);
00974                     z3 = MULTIPLY(-d7, FIX_1_961570560);
00975                     tmp1 = MULTIPLY(-d5, FIX_0_509795579);
00976                     z2 = MULTIPLY(-d5, FIX_2_562915447);
00977                     z4 = MULTIPLY(-d5, FIX_0_390180644);
00978                     z5 = MULTIPLY(d5 + d7, FIX_1_175875602);
00979                     
00980                     z3 += z5;
00981                     z4 += z5;
00982                     
00983                     tmp0 += z3;
00984                     tmp1 += z4;
00985                     tmp2 = z2 + z3;
00986                     tmp3 = z1 + z4;
00987                 }
00988             }
00989         } else {
00990             if (d3) {
00991                 if (d1) {
00992                     /* d1 != 0, d3 != 0, d5 == 0, d7 != 0 */
00993                     z1 = d7 + d1;
00994                     z3 = d7 + d3;
00995                     z5 = MULTIPLY(z3 + d1, FIX_1_175875602);
00996                     
00997                     tmp0 = MULTIPLY(d7, FIX_0_298631336); 
00998                     tmp2 = MULTIPLY(d3, FIX_3_072711026);
00999                     tmp3 = MULTIPLY(d1, FIX_1_501321110);
01000                     z1 = MULTIPLY(-z1, FIX_0_899976223);
01001                     z2 = MULTIPLY(-d3, FIX_2_562915447);
01002                     z3 = MULTIPLY(-z3, FIX_1_961570560);
01003                     z4 = MULTIPLY(-d1, FIX_0_390180644);
01004                     
01005                     z3 += z5;
01006                     z4 += z5;
01007                     
01008                     tmp0 += z1 + z3;
01009                     tmp1 = z2 + z4;
01010                     tmp2 += z2 + z3;
01011                     tmp3 += z1 + z4;
01012                 } else {
01013                     /* d1 == 0, d3 != 0, d5 == 0, d7 != 0 */
01014                     z3 = d7 + d3;
01015                     
01016                     tmp0 = MULTIPLY(-d7, FIX_0_601344887); 
01017                     z1 = MULTIPLY(-d7, FIX_0_899976223);
01018                     tmp2 = MULTIPLY(d3, FIX_0_509795579);
01019                     z2 = MULTIPLY(-d3, FIX_2_562915447);
01020                     z5 = MULTIPLY(z3, FIX_1_175875602);
01021                     z3 = MULTIPLY(-z3, FIX_0_785694958);
01022                     
01023                     tmp0 += z3;
01024                     tmp1 = z2 + z5;
01025                     tmp2 += z3;
01026                     tmp3 = z1 + z5;
01027                 }
01028             } else {
01029                 if (d1) {
01030                     /* d1 != 0, d3 == 0, d5 == 0, d7 != 0 */
01031                     z1 = d7 + d1;
01032                     z5 = MULTIPLY(z1, FIX_1_175875602);
01033 
01034                     z1 = MULTIPLY(z1, FIX_0_275899380);
01035                     z3 = MULTIPLY(-d7, FIX_1_961570560);
01036                     tmp0 = MULTIPLY(-d7, FIX_1_662939225); 
01037                     z4 = MULTIPLY(-d1, FIX_0_390180644);
01038                     tmp3 = MULTIPLY(d1, FIX_1_111140466);
01039 
01040                     tmp0 += z1;
01041                     tmp1 = z4 + z5;
01042                     tmp2 = z3 + z5;
01043                     tmp3 += z1;
01044                 } else {
01045                     /* d1 == 0, d3 == 0, d5 == 0, d7 != 0 */
01046                     tmp0 = MULTIPLY(-d7, FIX_1_387039845);
01047                     tmp1 = MULTIPLY(d7, FIX_1_175875602);
01048                     tmp2 = MULTIPLY(-d7, FIX_0_785694958);
01049                     tmp3 = MULTIPLY(d7, FIX_0_275899380);
01050                 }
01051             }
01052         }
01053     } else {
01054         if (d5) {
01055             if (d3) {
01056                 if (d1) {
01057                     /* d1 != 0, d3 != 0, d5 != 0, d7 == 0 */
01058                     z2 = d5 + d3;
01059                     z4 = d5 + d1;
01060                     z5 = MULTIPLY(d3 + z4, FIX_1_175875602);
01061                     
01062                     tmp1 = MULTIPLY(d5, FIX_2_053119869);
01063                     tmp2 = MULTIPLY(d3, FIX_3_072711026);
01064                     tmp3 = MULTIPLY(d1, FIX_1_501321110);
01065                     z1 = MULTIPLY(-d1, FIX_0_899976223);
01066                     z2 = MULTIPLY(-z2, FIX_2_562915447);
01067                     z3 = MULTIPLY(-d3, FIX_1_961570560);
01068                     z4 = MULTIPLY(-z4, FIX_0_390180644);
01069                     
01070                     z3 += z5;
01071                     z4 += z5;
01072                     
01073                     tmp0 = z1 + z3;
01074                     tmp1 += z2 + z4;
01075                     tmp2 += z2 + z3;
01076                     tmp3 += z1 + z4;
01077                 } else {
01078                     /* d1 == 0, d3 != 0, d5 != 0, d7 == 0 */
01079                     z2 = d5 + d3;
01080                     
01081                     z5 = MULTIPLY(z2, FIX_1_175875602);
01082                     tmp1 = MULTIPLY(d5, FIX_1_662939225);
01083                     z4 = MULTIPLY(-d5, FIX_0_390180644);
01084                     z2 = MULTIPLY(-z2, FIX_1_387039845);
01085                     tmp2 = MULTIPLY(d3, FIX_1_111140466);
01086                     z3 = MULTIPLY(-d3, FIX_1_961570560);
01087                     
01088                     tmp0 = z3 + z5;
01089                     tmp1 += z2;
01090                     tmp2 += z2;
01091                     tmp3 = z4 + z5;
01092                 }
01093             } else {
01094                 if (d1) {
01095                     /* d1 != 0, d3 == 0, d5 != 0, d7 == 0 */
01096                     z4 = d5 + d1;
01097                     
01098                     z5 = MULTIPLY(z4, FIX_1_175875602);
01099                     z1 = MULTIPLY(-d1, FIX_0_899976223);
01100                     tmp3 = MULTIPLY(d1, FIX_0_601344887);
01101                     tmp1 = MULTIPLY(-d5, FIX_0_509795579);
01102                     z2 = MULTIPLY(-d5, FIX_2_562915447);
01103                     z4 = MULTIPLY(z4, FIX_0_785694958);
01104                     
01105                     tmp0 = z1 + z5;
01106                     tmp1 += z4;
01107                     tmp2 = z2 + z5;
01108                     tmp3 += z4;
01109                 } else {
01110                     /* d1 == 0, d3 == 0, d5 != 0, d7 == 0 */
01111                     tmp0 = MULTIPLY(d5, FIX_1_175875602);
01112                     tmp1 = MULTIPLY(d5, FIX_0_275899380);
01113                     tmp2 = MULTIPLY(-d5, FIX_1_387039845);
01114                     tmp3 = MULTIPLY(d5, FIX_0_785694958);
01115                 }
01116             }
01117         } else {
01118             if (d3) {
01119                 if (d1) {
01120                     /* d1 != 0, d3 != 0, d5 == 0, d7 == 0 */
01121                     z5 = d1 + d3;
01122                     tmp3 = MULTIPLY(d1, FIX_0_211164243);
01123                     tmp2 = MULTIPLY(-d3, FIX_1_451774981);
01124                     z1 = MULTIPLY(d1, FIX_1_061594337);
01125                     z2 = MULTIPLY(-d3, FIX_2_172734803);
01126                     z4 = MULTIPLY(z5, FIX_0_785694958);
01127                     z5 = MULTIPLY(z5, FIX_1_175875602);
01128                     
01129                     tmp0 = z1 - z4;
01130                     tmp1 = z2 + z4;
01131                     tmp2 += z5;
01132                     tmp3 += z5;
01133                 } else {
01134                     /* d1 == 0, d3 != 0, d5 == 0, d7 == 0 */
01135                     tmp0 = MULTIPLY(-d3, FIX_0_785694958);
01136                     tmp1 = MULTIPLY(-d3, FIX_1_387039845);
01137                     tmp2 = MULTIPLY(-d3, FIX_0_275899380);
01138                     tmp3 = MULTIPLY(d3, FIX_1_175875602);
01139                 }
01140             } else {
01141                 if (d1) {
01142                     /* d1 != 0, d3 == 0, d5 == 0, d7 == 0 */
01143                     tmp0 = MULTIPLY(d1, FIX_0_275899380);
01144                     tmp1 = MULTIPLY(d1, FIX_0_785694958);
01145                     tmp2 = MULTIPLY(d1, FIX_1_175875602);
01146                     tmp3 = MULTIPLY(d1, FIX_1_387039845);
01147                 } else {
01148                     /* d1 == 0, d3 == 0, d5 == 0, d7 == 0 */
01149                     tmp0 = tmp1 = tmp2 = tmp3 = 0;
01150                 }
01151             }
01152         }
01153     }
01154 
01155     /* Final output stage: inputs are tmp10..tmp13, tmp0..tmp3 */
01156 
01157     dataptr[DCTSIZE*0] = (DCTELEM) DESCALE(tmp10 + tmp3,
01158                                            CONST_BITS+PASS1_BITS+3);
01159     dataptr[DCTSIZE*7] = (DCTELEM) DESCALE(tmp10 - tmp3,
01160                                            CONST_BITS+PASS1_BITS+3);
01161     dataptr[DCTSIZE*1] = (DCTELEM) DESCALE(tmp11 + tmp2,
01162                                            CONST_BITS+PASS1_BITS+3);
01163     dataptr[DCTSIZE*6] = (DCTELEM) DESCALE(tmp11 - tmp2,
01164                                            CONST_BITS+PASS1_BITS+3);
01165     dataptr[DCTSIZE*2] = (DCTELEM) DESCALE(tmp12 + tmp1,
01166                                            CONST_BITS+PASS1_BITS+3);
01167     dataptr[DCTSIZE*5] = (DCTELEM) DESCALE(tmp12 - tmp1,
01168                                            CONST_BITS+PASS1_BITS+3);
01169     dataptr[DCTSIZE*3] = (DCTELEM) DESCALE(tmp13 + tmp0,
01170                                            CONST_BITS+PASS1_BITS+3);
01171     dataptr[DCTSIZE*4] = (DCTELEM) DESCALE(tmp13 - tmp0,
01172                                            CONST_BITS+PASS1_BITS+3);
01173     
01174     dataptr++;                  /* advance pointer to next column */
01175   }
01176 }
01177 
01178 
01179 /* here is the reference one, in case of problems with the normal one */
01180 
01181 /* idctref.c, Inverse Discrete Fourier Transform, double precision          */
01182 
01183 /* Copyright (C) 1994, MPEG Software Simulation Group. All Rights Reserved. */
01184 
01185 /*
01186  * Disclaimer of Warranty
01187  *
01188  * These software programs are available to the user without any license fee or
01189  * royalty on an "as is" basis.  The MPEG Software Simulation Group disclaims
01190  * any and all warranties, whether express, implied, or statuary, including any
01191  * implied warranties or merchantability or of fitness for a particular
01192  * purpose.  In no event shall the copyright-holder be liable for any
01193  * incidental, punitive, or consequential damages of any kind whatsoever
01194  * arising from the use of these programs.
01195  *
01196  * This disclaimer of warranty extends to the user of these programs and user's
01197  * customers, employees, agents, transferees, successors, and assigns.
01198  *
01199  * The MPEG Software Simulation Group does not represent or warrant that the
01200  * programs furnished hereunder are free of infringement of any third-party
01201  * patents.
01202  *
01203  * Commercial implementations of MPEG-1 and MPEG-2 video, including shareware,
01204  * are subject to royalty fees to patent holders.  Many of these patents are
01205  * general enough such that they are unavoidable regardless of implementation
01206  * design.
01207  *
01208  */
01209 
01210 /*  Perform IEEE 1180 reference (64-bit floating point, separable 8x1
01211  *  direct matrix multiply) Inverse Discrete Cosine Transform
01212 */
01213 
01214 
01215 /* Here we use math.h to generate constants.  Compiler results may
01216    vary a little */
01217 
01218 #ifndef PI
01219 #ifdef M_PI
01220 #define PI M_PI
01221 #else
01222 #define PI 3.14159265358979323846
01223 #endif
01224 #endif
01225 
01226 /* cosine transform matrix for 8x1 IDCT */
01227 static double itrans_coef[8][8];
01228 
01229 /* initialize DCT coefficient matrix */
01230 
01231 void init_idctref()
01232 {
01233   int freq, time;
01234   double scale;
01235 
01236   for (freq=0; freq < 8; freq++)
01237   {
01238     scale = (freq == 0) ? sqrt(0.125) : 0.5;
01239     for (time=0; time<8; time++)
01240       itrans_coef[freq][time] = scale*cos((PI/8.0)*freq*(time + 0.5));
01241   }
01242 }
01243 
01244 /* perform IDCT matrix multiply for 8x8 coefficient block */
01245 
01246 void reference_rev_dct(block)
01247 int16 *block;
01248 {
01249   int i, j, k, v;
01250   double partial_product;
01251   double tmp[64];
01252 
01253   for (i=0; i<8; i++)
01254     for (j=0; j<8; j++)
01255     {
01256       partial_product = 0.0;
01257 
01258       for (k=0; k<8; k++)
01259         partial_product+= itrans_coef[k][j]*block[8*i+k];
01260 
01261       tmp[8*i+j] = partial_product;
01262     }
01263 
01264   /* Transpose operation is integrated into address mapping by switching 
01265      loop order of i and j */
01266 
01267   for (j=0; j<8; j++)
01268     for (i=0; i<8; i++)
01269     {
01270       partial_product = 0.0;
01271 
01272       for (k=0; k<8; k++)
01273         partial_product+= itrans_coef[k][i]*tmp[8*k+j];
01274 
01275       v = floor(partial_product+0.5);
01276       block[8*i+j] = (v<-256) ? -256 : ((v>255) ? 255 : v);
01277     }
01278 }