The Ecere SDK - sdk/blob - deps/jpeg-9a/jfdctint.c

   1 /*
   2  * jfdctint.c
   3  *
   4  * Copyright (C) 1991-1996, Thomas G. Lane.
   5  * Modification developed 2003-2013 by Guido Vollbeding.
   6  * This file is part of the Independent JPEG Group's software.
   7  * For conditions of distribution and use, see the accompanying README file.
   8  *
   9  * This file contains a slow-but-accurate integer implementation of the
  10  * forward DCT (Discrete Cosine Transform).
  11  *
  12  * A 2-D DCT can be done by 1-D DCT on each row followed by 1-D DCT
  13  * on each column.  Direct algorithms are also available, but they are
  14  * much more complex and seem not to be any faster when reduced to code.
  15  *
  16  * This implementation is based on an algorithm described in
  17  *   C. Loeffler, A. Ligtenberg and G. Moschytz, "Practical Fast 1-D DCT
  18  *   Algorithms with 11 Multiplications", Proc. Int'l. Conf. on Acoustics,
  19  *   Speech, and Signal Processing 1989 (ICASSP '89), pp. 988-991.
  20  * The primary algorithm described there uses 11 multiplies and 29 adds.
  21  * We use their alternate method with 12 multiplies and 32 adds.
  22  * The advantage of this method is that no data path contains more than one
  23  * multiplication; this allows a very simple and accurate implementation in
  24  * scaled fixed-point arithmetic, with a minimal number of shifts.
  25  *
  26  * We also provide FDCT routines with various input sample block sizes for
  27  * direct resolution reduction or enlargement and for direct resolving the
  28  * common 2x1 and 1x2 subsampling cases without additional resampling: NxN
  29  * (N=1...16), 2NxN, and Nx2N (N=1...8) pixels for one 8x8 output DCT block.
  30  *
  31  * For N<8 we fill the remaining block coefficients with zero.
  32  * For N>8 we apply a partial N-point FDCT on the input samples, computing
  33  * just the lower 8 frequency coefficients and discarding the rest.
  34  *
  35  * We must scale the output coefficients of the N-point FDCT appropriately
  36  * to the standard 8-point FDCT level by 8/N per 1-D pass.  This scaling
  37  * is folded into the constant multipliers (pass 2) and/or final/initial
  38  * shifting.
  39  *
  40  * CAUTION: We rely on the FIX() macro except for the N=1,2,4,8 cases
  41  * since there would be too many additional constants to pre-calculate.
  42  */
  43
  44 #define JPEG_INTERNALS
  45 #include "jinclude.h"
  46 #include "jpeglib.h"
  47 #include "jdct.h"               /* Private declarations for DCT subsystem */
  48
  49 #ifdef DCT_ISLOW_SUPPORTED
  50
  51
  52 /*
  53  * This module is specialized to the case DCTSIZE = 8.
  54  */
  55
  56 #if DCTSIZE != 8
  57   Sorry, this code only copes with 8x8 DCT blocks. /* deliberate syntax err */
  58 #endif
  59
  60
  61 /*
  62  * The poop on this scaling stuff is as follows:
  63  *
  64  * Each 1-D DCT step produces outputs which are a factor of sqrt(N)
  65  * larger than the true DCT outputs.  The final outputs are therefore
  66  * a factor of N larger than desired; since N=8 this can be cured by
  67  * a simple right shift at the end of the algorithm.  The advantage of
  68  * this arrangement is that we save two multiplications per 1-D DCT,
  69  * because the y0 and y4 outputs need not be divided by sqrt(N).
  70  * In the IJG code, this factor of 8 is removed by the quantization step
  71  * (in jcdctmgr.c), NOT in this module.
  72  *
  73  * We have to do addition and subtraction of the integer inputs, which
  74  * is no problem, and multiplication by fractional constants, which is
  75  * a problem to do in integer arithmetic.  We multiply all the constants
  76  * by CONST_SCALE and convert them to integer constants (thus retaining
  77  * CONST_BITS bits of precision in the constants).  After doing a
  78  * multiplication we have to divide the product by CONST_SCALE, with proper
  79  * rounding, to produce the correct output.  This division can be done
  80  * cheaply as a right shift of CONST_BITS bits.  We postpone shifting
  81  * as long as possible so that partial sums can be added together with
  82  * full fractional precision.
  83  *
  84  * The outputs of the first pass are scaled up by PASS1_BITS bits so that
  85  * they are represented to better-than-integral precision.  These outputs
  86  * require BITS_IN_JSAMPLE + PASS1_BITS + 3 bits; this fits in a 16-bit word
  87  * with the recommended scaling.  (For 12-bit sample data, the intermediate
  88  * array is INT32 anyway.)
  89  *
  90  * To avoid overflow of the 32-bit intermediate results in pass 2, we must
  91  * have BITS_IN_JSAMPLE + CONST_BITS + PASS1_BITS <= 26.  Error analysis
  92  * shows that the values given below are the most effective.
  93  */
  94
  95 #if BITS_IN_JSAMPLE == 8
  96 #define CONST_BITS  13
  97 #define PASS1_BITS  2
  98 #else
  99 #define CONST_BITS  13
 100 #define PASS1_BITS  1           /* lose a little precision to avoid overflow */
 101 #endif
 102
 103 /* Some C compilers fail to reduce "FIX(constant)" at compile time, thus
 104  * causing a lot of useless floating-point operations at run time.
 105  * To get around this we use the following pre-calculated constants.
 106  * If you change CONST_BITS you may want to add appropriate values.
 107  * (With a reasonable C compiler, you can just rely on the FIX() macro...)
 108  */
 109
 110 #if CONST_BITS == 13
 111 #define FIX_0_298631336  ((INT32)  2446)        /* FIX(0.298631336) */
 112 #define FIX_0_390180644  ((INT32)  3196)        /* FIX(0.390180644) */
 113 #define FIX_0_541196100  ((INT32)  4433)        /* FIX(0.541196100) */
 114 #define FIX_0_765366865  ((INT32)  6270)        /* FIX(0.765366865) */
 115 #define FIX_0_899976223  ((INT32)  7373)        /* FIX(0.899976223) */
 116 #define FIX_1_175875602  ((INT32)  9633)        /* FIX(1.175875602) */
 117 #define FIX_1_501321110  ((INT32)  12299)       /* FIX(1.501321110) */
 118 #define FIX_1_847759065  ((INT32)  15137)       /* FIX(1.847759065) */
 119 #define FIX_1_961570560  ((INT32)  16069)       /* FIX(1.961570560) */
 120 #define FIX_2_053119869  ((INT32)  16819)       /* FIX(2.053119869) */
 121 #define FIX_2_562915447  ((INT32)  20995)       /* FIX(2.562915447) */
 122 #define FIX_3_072711026  ((INT32)  25172)       /* FIX(3.072711026) */
 123 #else
 124 #define FIX_0_298631336  FIX(0.298631336)
 125 #define FIX_0_390180644  FIX(0.390180644)
 126 #define FIX_0_541196100  FIX(0.541196100)
 127 #define FIX_0_765366865  FIX(0.765366865)
 128 #define FIX_0_899976223  FIX(0.899976223)
 129 #define FIX_1_175875602  FIX(1.175875602)
 130 #define FIX_1_501321110  FIX(1.501321110)
 131 #define FIX_1_847759065  FIX(1.847759065)
 132 #define FIX_1_961570560  FIX(1.961570560)
 133 #define FIX_2_053119869  FIX(2.053119869)
 134 #define FIX_2_562915447  FIX(2.562915447)
 135 #define FIX_3_072711026  FIX(3.072711026)
 136 #endif
 137
 138
 139 /* Multiply an INT32 variable by an INT32 constant to yield an INT32 result.
 140  * For 8-bit samples with the recommended scaling, all the variable
 141  * and constant values involved are no more than 16 bits wide, so a
 142  * 16x16->32 bit multiply can be used instead of a full 32x32 multiply.
 143  * For 12-bit samples, a full 32-bit multiplication will be needed.
 144  */
 145
 146 #if BITS_IN_JSAMPLE == 8
 147 #define MULTIPLY(var,const)  MULTIPLY16C16(var,const)
 148 #else
 149 #define MULTIPLY(var,const)  ((var) * (const))
 150 #endif
 151
 152
 153 /*
 154  * Perform the forward DCT on one block of samples.
 155  */
 156
 157 GLOBAL(void)
 158 jpeg_fdct_islow (DCTELEM * data, JSAMPARRAY sample_data, JDIMENSION start_col)
 159 {
 160   INT32 tmp0, tmp1, tmp2, tmp3;
 161   INT32 tmp10, tmp11, tmp12, tmp13;
 162   INT32 z1;
 163   DCTELEM *dataptr;
 164   JSAMPROW elemptr;
 165   int ctr;
 166   SHIFT_TEMPS
 167
 168   /* Pass 1: process rows.
 169    * Note results are scaled up by sqrt(8) compared to a true DCT;
 170    * furthermore, we scale the results by 2**PASS1_BITS.
 171    * cK represents sqrt(2) * cos(K*pi/16).
 172    */
 173
 174   dataptr = data;
 175   for (ctr = 0; ctr < DCTSIZE; ctr++) {
 176     elemptr = sample_data[ctr] + start_col;
 177
 178     /* Even part per LL&M figure 1 --- note that published figure is faulty;
 179      * rotator "c1" should be "c6".
 180      */
 181
 182     tmp0 = GETJSAMPLE(elemptr[0]) + GETJSAMPLE(elemptr[7]);
 183     tmp1 = GETJSAMPLE(elemptr[1]) + GETJSAMPLE(elemptr[6]);
 184     tmp2 = GETJSAMPLE(elemptr[2]) + GETJSAMPLE(elemptr[5]);
 185     tmp3 = GETJSAMPLE(elemptr[3]) + GETJSAMPLE(elemptr[4]);
 186
 187     tmp10 = tmp0 + tmp3;
 188     tmp12 = tmp0 - tmp3;
 189     tmp11 = tmp1 + tmp2;
 190     tmp13 = tmp1 - tmp2;
 191
 192     tmp0 = GETJSAMPLE(elemptr[0]) - GETJSAMPLE(elemptr[7]);
 193     tmp1 = GETJSAMPLE(elemptr[1]) - GETJSAMPLE(elemptr[6]);
 194     tmp2 = GETJSAMPLE(elemptr[2]) - GETJSAMPLE(elemptr[5]);
 195     tmp3 = GETJSAMPLE(elemptr[3]) - GETJSAMPLE(elemptr[4]);
 196
 197     /* Apply unsigned->signed conversion */
 198     dataptr[0] = (DCTELEM) ((tmp10 + tmp11 - 8 * CENTERJSAMPLE) << PASS1_BITS);
 199     dataptr[4] = (DCTELEM) ((tmp10 - tmp11) << PASS1_BITS);
 200
 201     z1 = MULTIPLY(tmp12 + tmp13, FIX_0_541196100);       /* c6 */
 202     /* Add fudge factor here for final descale. */
 203     z1 += ONE << (CONST_BITS-PASS1_BITS-1);
 204
 205     dataptr[2] = (DCTELEM)
 206       RIGHT_SHIFT(z1 + MULTIPLY(tmp12, FIX_0_765366865), /* c2-c6 */
 207                   CONST_BITS-PASS1_BITS);
 208     dataptr[6] = (DCTELEM)
 209       RIGHT_SHIFT(z1 - MULTIPLY(tmp13, FIX_1_847759065), /* c2+c6 */
 210                   CONST_BITS-PASS1_BITS);
 211
 212     /* Odd part per figure 8 --- note paper omits factor of sqrt(2).
 213      * i0..i3 in the paper are tmp0..tmp3 here.
 214      */
 215
 216     tmp12 = tmp0 + tmp2;
 217     tmp13 = tmp1 + tmp3;
 218
 219     z1 = MULTIPLY(tmp12 + tmp13, FIX_1_175875602);       /*  c3 */
 220     /* Add fudge factor here for final descale. */
 221     z1 += ONE << (CONST_BITS-PASS1_BITS-1);
 222
 223     tmp12 = MULTIPLY(tmp12, - FIX_0_390180644);          /* -c3+c5 */
 224     tmp13 = MULTIPLY(tmp13, - FIX_1_961570560);          /* -c3-c5 */
 225     tmp12 += z1;
 226     tmp13 += z1;
 227
 228     z1 = MULTIPLY(tmp0 + tmp3, - FIX_0_899976223);       /* -c3+c7 */
 229     tmp0 = MULTIPLY(tmp0, FIX_1_501321110);              /*  c1+c3-c5-c7 */
 230     tmp3 = MULTIPLY(tmp3, FIX_0_298631336);              /* -c1+c3+c5-c7 */
 231     tmp0 += z1 + tmp12;
 232     tmp3 += z1 + tmp13;
 233
 234     z1 = MULTIPLY(tmp1 + tmp2, - FIX_2_562915447);       /* -c1-c3 */
 235     tmp1 = MULTIPLY(tmp1, FIX_3_072711026);              /*  c1+c3+c5-c7 */
 236     tmp2 = MULTIPLY(tmp2, FIX_2_053119869);              /*  c1+c3-c5+c7 */
 237     tmp1 += z1 + tmp13;
 238     tmp2 += z1 + tmp12;
 239
 240     dataptr[1] = (DCTELEM) RIGHT_SHIFT(tmp0, CONST_BITS-PASS1_BITS);
 241     dataptr[3] = (DCTELEM) RIGHT_SHIFT(tmp1, CONST_BITS-PASS1_BITS);
 242     dataptr[5] = (DCTELEM) RIGHT_SHIFT(tmp2, CONST_BITS-PASS1_BITS);
 243     dataptr[7] = (DCTELEM) RIGHT_SHIFT(tmp3, CONST_BITS-PASS1_BITS);
 244
 245     dataptr += DCTSIZE;         /* advance pointer to next row */
 246   }
 247
 248   /* Pass 2: process columns.
 249    * We remove the PASS1_BITS scaling, but leave the results scaled up
 250    * by an overall factor of 8.
 251    * cK represents sqrt(2) * cos(K*pi/16).
 252    */
 253
 254   dataptr = data;
 255   for (ctr = DCTSIZE-1; ctr >= 0; ctr--) {
 256     /* Even part per LL&M figure 1 --- note that published figure is faulty;
 257      * rotator "c1" should be "c6".
 258      */
 259
 260     tmp0 = dataptr[DCTSIZE*0] + dataptr[DCTSIZE*7];
 261     tmp1 = dataptr[DCTSIZE*1] + dataptr[DCTSIZE*6];
 262     tmp2 = dataptr[DCTSIZE*2] + dataptr[DCTSIZE*5];
 263     tmp3 = dataptr[DCTSIZE*3] + dataptr[DCTSIZE*4];
 264
 265     /* Add fudge factor here for final descale. */
 266     tmp10 = tmp0 + tmp3 + (ONE << (PASS1_BITS-1));
 267     tmp12 = tmp0 - tmp3;
 268     tmp11 = tmp1 + tmp2;
 269     tmp13 = tmp1 - tmp2;
 270
 271     tmp0 = dataptr[DCTSIZE*0] - dataptr[DCTSIZE*7];
 272     tmp1 = dataptr[DCTSIZE*1] - dataptr[DCTSIZE*6];
 273     tmp2 = dataptr[DCTSIZE*2] - dataptr[DCTSIZE*5];
 274     tmp3 = dataptr[DCTSIZE*3] - dataptr[DCTSIZE*4];
 275
 276     dataptr[DCTSIZE*0] = (DCTELEM) RIGHT_SHIFT(tmp10 + tmp11, PASS1_BITS);
 277     dataptr[DCTSIZE*4] = (DCTELEM) RIGHT_SHIFT(tmp10 - tmp11, PASS1_BITS);
 278
 279     z1 = MULTIPLY(tmp12 + tmp13, FIX_0_541196100);       /* c6 */
 280     /* Add fudge factor here for final descale. */
 281     z1 += ONE << (CONST_BITS+PASS1_BITS-1);
 282
 283     dataptr[DCTSIZE*2] = (DCTELEM)
 284       RIGHT_SHIFT(z1 + MULTIPLY(tmp12, FIX_0_765366865), /* c2-c6 */
 285                   CONST_BITS+PASS1_BITS);
 286     dataptr[DCTSIZE*6] = (DCTELEM)
 287       RIGHT_SHIFT(z1 - MULTIPLY(tmp13, FIX_1_847759065), /* c2+c6 */
 288                   CONST_BITS+PASS1_BITS);
 289
 290     /* Odd part per figure 8 --- note paper omits factor of sqrt(2).
 291      * i0..i3 in the paper are tmp0..tmp3 here.
 292      */
 293
 294     tmp12 = tmp0 + tmp2;
 295     tmp13 = tmp1 + tmp3;
 296
 297     z1 = MULTIPLY(tmp12 + tmp13, FIX_1_175875602);       /*  c3 */
 298     /* Add fudge factor here for final descale. */
 299     z1 += ONE << (CONST_BITS+PASS1_BITS-1);
 300
 301     tmp12 = MULTIPLY(tmp12, - FIX_0_390180644);          /* -c3+c5 */
 302     tmp13 = MULTIPLY(tmp13, - FIX_1_961570560);          /* -c3-c5 */
 303     tmp12 += z1;
 304     tmp13 += z1;
 305
 306     z1 = MULTIPLY(tmp0 + tmp3, - FIX_0_899976223);       /* -c3+c7 */
 307     tmp0 = MULTIPLY(tmp0, FIX_1_501321110);              /*  c1+c3-c5-c7 */
 308     tmp3 = MULTIPLY(tmp3, FIX_0_298631336);              /* -c1+c3+c5-c7 */
 309     tmp0 += z1 + tmp12;
 310     tmp3 += z1 + tmp13;
 311
 312     z1 = MULTIPLY(tmp1 + tmp2, - FIX_2_562915447);       /* -c1-c3 */
 313     tmp1 = MULTIPLY(tmp1, FIX_3_072711026);              /*  c1+c3+c5-c7 */
 314     tmp2 = MULTIPLY(tmp2, FIX_2_053119869);              /*  c1+c3-c5+c7 */
 315     tmp1 += z1 + tmp13;
 316     tmp2 += z1 + tmp12;
 317
 318     dataptr[DCTSIZE*1] = (DCTELEM) RIGHT_SHIFT(tmp0, CONST_BITS+PASS1_BITS);
 319     dataptr[DCTSIZE*3] = (DCTELEM) RIGHT_SHIFT(tmp1, CONST_BITS+PASS1_BITS);
 320     dataptr[DCTSIZE*5] = (DCTELEM) RIGHT_SHIFT(tmp2, CONST_BITS+PASS1_BITS);
 321     dataptr[DCTSIZE*7] = (DCTELEM) RIGHT_SHIFT(tmp3, CONST_BITS+PASS1_BITS);
 322
 323     dataptr++;                  /* advance pointer to next column */
 324   }
 325 }
 326
 327 #ifdef DCT_SCALING_SUPPORTED
 328
 329
 330 /*
 331  * Perform the forward DCT on a 7x7 sample block.
 332  */
 333
 334 GLOBAL(void)
 335 jpeg_fdct_7x7 (DCTELEM * data, JSAMPARRAY sample_data, JDIMENSION start_col)
 336 {
 337   INT32 tmp0, tmp1, tmp2, tmp3;
 338   INT32 tmp10, tmp11, tmp12;
 339   INT32 z1, z2, z3;
 340   DCTELEM *dataptr;
 341   JSAMPROW elemptr;
 342   int ctr;
 343   SHIFT_TEMPS
 344
 345   /* Pre-zero output coefficient block. */
 346   MEMZERO(data, SIZEOF(DCTELEM) * DCTSIZE2);
 347
 348   /* Pass 1: process rows.
 349    * Note results are scaled up by sqrt(8) compared to a true DCT;
 350    * furthermore, we scale the results by 2**PASS1_BITS.
 351    * cK represents sqrt(2) * cos(K*pi/14).
 352    */
 353
 354   dataptr = data;
 355   for (ctr = 0; ctr < 7; ctr++) {
 356     elemptr = sample_data[ctr] + start_col;
 357
 358     /* Even part */
 359
 360     tmp0 = GETJSAMPLE(elemptr[0]) + GETJSAMPLE(elemptr[6]);
 361     tmp1 = GETJSAMPLE(elemptr[1]) + GETJSAMPLE(elemptr[5]);
 362     tmp2 = GETJSAMPLE(elemptr[2]) + GETJSAMPLE(elemptr[4]);
 363     tmp3 = GETJSAMPLE(elemptr[3]);
 364
 365     tmp10 = GETJSAMPLE(elemptr[0]) - GETJSAMPLE(elemptr[6]);
 366     tmp11 = GETJSAMPLE(elemptr[1]) - GETJSAMPLE(elemptr[5]);
 367     tmp12 = GETJSAMPLE(elemptr[2]) - GETJSAMPLE(elemptr[4]);
 368
 369     z1 = tmp0 + tmp2;
 370     /* Apply unsigned->signed conversion */
 371     dataptr[0] = (DCTELEM)
 372       ((z1 + tmp1 + tmp3 - 7 * CENTERJSAMPLE) << PASS1_BITS);
 373     tmp3 += tmp3;
 374     z1 -= tmp3;
 375     z1 -= tmp3;
 376     z1 = MULTIPLY(z1, FIX(0.353553391));                /* (c2+c6-c4)/2 */
 377     z2 = MULTIPLY(tmp0 - tmp2, FIX(0.920609002));       /* (c2+c4-c6)/2 */
 378     z3 = MULTIPLY(tmp1 - tmp2, FIX(0.314692123));       /* c6 */
 379     dataptr[2] = (DCTELEM) DESCALE(z1 + z2 + z3, CONST_BITS-PASS1_BITS);
 380     z1 -= z2;
 381     z2 = MULTIPLY(tmp0 - tmp1, FIX(0.881747734));       /* c4 */
 382     dataptr[4] = (DCTELEM)
 383       DESCALE(z2 + z3 - MULTIPLY(tmp1 - tmp3, FIX(0.707106781)), /* c2+c6-c4 */
 384               CONST_BITS-PASS1_BITS);
 385     dataptr[6] = (DCTELEM) DESCALE(z1 + z2, CONST_BITS-PASS1_BITS);
 386
 387     /* Odd part */
 388
 389     tmp1 = MULTIPLY(tmp10 + tmp11, FIX(0.935414347));   /* (c3+c1-c5)/2 */
 390     tmp2 = MULTIPLY(tmp10 - tmp11, FIX(0.170262339));   /* (c3+c5-c1)/2 */
 391     tmp0 = tmp1 - tmp2;
 392     tmp1 += tmp2;
 393     tmp2 = MULTIPLY(tmp11 + tmp12, - FIX(1.378756276)); /* -c1 */
 394     tmp1 += tmp2;
 395     tmp3 = MULTIPLY(tmp10 + tmp12, FIX(0.613604268));   /* c5 */
 396     tmp0 += tmp3;
 397     tmp2 += tmp3 + MULTIPLY(tmp12, FIX(1.870828693));   /* c3+c1-c5 */
 398
 399     dataptr[1] = (DCTELEM) DESCALE(tmp0, CONST_BITS-PASS1_BITS);
 400     dataptr[3] = (DCTELEM) DESCALE(tmp1, CONST_BITS-PASS1_BITS);
 401     dataptr[5] = (DCTELEM) DESCALE(tmp2, CONST_BITS-PASS1_BITS);
 402
 403     dataptr += DCTSIZE;         /* advance pointer to next row */
 404   }
 405
 406   /* Pass 2: process columns.
 407    * We remove the PASS1_BITS scaling, but leave the results scaled up
 408    * by an overall factor of 8.
 409    * We must also scale the output by (8/7)**2 = 64/49, which we fold
 410    * into the constant multipliers:
 411    * cK now represents sqrt(2) * cos(K*pi/14) * 64/49.
 412    */
 413
 414   dataptr = data;
 415   for (ctr = 0; ctr < 7; ctr++) {
 416     /* Even part */
 417
 418     tmp0 = dataptr[DCTSIZE*0] + dataptr[DCTSIZE*6];
 419     tmp1 = dataptr[DCTSIZE*1] + dataptr[DCTSIZE*5];
 420     tmp2 = dataptr[DCTSIZE*2] + dataptr[DCTSIZE*4];
 421     tmp3 = dataptr[DCTSIZE*3];
 422
 423     tmp10 = dataptr[DCTSIZE*0] - dataptr[DCTSIZE*6];
 424     tmp11 = dataptr[DCTSIZE*1] - dataptr[DCTSIZE*5];
 425     tmp12 = dataptr[DCTSIZE*2] - dataptr[DCTSIZE*4];
 426
 427     z1 = tmp0 + tmp2;
 428     dataptr[DCTSIZE*0] = (DCTELEM)
 429       DESCALE(MULTIPLY(z1 + tmp1 + tmp3, FIX(1.306122449)), /* 64/49 */
 430               CONST_BITS+PASS1_BITS);
 431     tmp3 += tmp3;
 432     z1 -= tmp3;
 433     z1 -= tmp3;
 434     z1 = MULTIPLY(z1, FIX(0.461784020));                /* (c2+c6-c4)/2 */
 435     z2 = MULTIPLY(tmp0 - tmp2, FIX(1.202428084));       /* (c2+c4-c6)/2 */
 436     z3 = MULTIPLY(tmp1 - tmp2, FIX(0.411026446));       /* c6 */
 437     dataptr[DCTSIZE*2] = (DCTELEM) DESCALE(z1 + z2 + z3, CONST_BITS+PASS1_BITS);
 438     z1 -= z2;
 439     z2 = MULTIPLY(tmp0 - tmp1, FIX(1.151670509));       /* c4 */
 440     dataptr[DCTSIZE*4] = (DCTELEM)
 441       DESCALE(z2 + z3 - MULTIPLY(tmp1 - tmp3, FIX(0.923568041)), /* c2+c6-c4 */
 442               CONST_BITS+PASS1_BITS);
 443     dataptr[DCTSIZE*6] = (DCTELEM) DESCALE(z1 + z2, CONST_BITS+PASS1_BITS);
 444
 445     /* Odd part */
 446
 447     tmp1 = MULTIPLY(tmp10 + tmp11, FIX(1.221765677));   /* (c3+c1-c5)/2 */
 448     tmp2 = MULTIPLY(tmp10 - tmp11, FIX(0.222383464));   /* (c3+c5-c1)/2 */
 449     tmp0 = tmp1 - tmp2;
 450     tmp1 += tmp2;
 451     tmp2 = MULTIPLY(tmp11 + tmp12, - FIX(1.800824523)); /* -c1 */
 452     tmp1 += tmp2;
 453     tmp3 = MULTIPLY(tmp10 + tmp12, FIX(0.801442310));   /* c5 */
 454     tmp0 += tmp3;
 455     tmp2 += tmp3 + MULTIPLY(tmp12, FIX(2.443531355));   /* c3+c1-c5 */
 456
 457     dataptr[DCTSIZE*1] = (DCTELEM) DESCALE(tmp0, CONST_BITS+PASS1_BITS);
 458     dataptr[DCTSIZE*3] = (DCTELEM) DESCALE(tmp1, CONST_BITS+PASS1_BITS);
 459     dataptr[DCTSIZE*5] = (DCTELEM) DESCALE(tmp2, CONST_BITS+PASS1_BITS);
 460
 461     dataptr++;                  /* advance pointer to next column */
 462   }
 463 }
 464
 465
 466 /*
 467  * Perform the forward DCT on a 6x6 sample block.
 468  */
 469
 470 GLOBAL(void)
 471 jpeg_fdct_6x6 (DCTELEM * data, JSAMPARRAY sample_data, JDIMENSION start_col)
 472 {
 473   INT32 tmp0, tmp1, tmp2;
 474   INT32 tmp10, tmp11, tmp12;
 475   DCTELEM *dataptr;
 476   JSAMPROW elemptr;
 477   int ctr;
 478   SHIFT_TEMPS
 479
 480   /* Pre-zero output coefficient block. */
 481   MEMZERO(data, SIZEOF(DCTELEM) * DCTSIZE2);
 482
 483   /* Pass 1: process rows.
 484    * Note results are scaled up by sqrt(8) compared to a true DCT;
 485    * furthermore, we scale the results by 2**PASS1_BITS.
 486    * cK represents sqrt(2) * cos(K*pi/12).
 487    */
 488
 489   dataptr = data;
 490   for (ctr = 0; ctr < 6; ctr++) {
 491     elemptr = sample_data[ctr] + start_col;
 492
 493     /* Even part */
 494
 495     tmp0 = GETJSAMPLE(elemptr[0]) + GETJSAMPLE(elemptr[5]);
 496     tmp11 = GETJSAMPLE(elemptr[1]) + GETJSAMPLE(elemptr[4]);
 497     tmp2 = GETJSAMPLE(elemptr[2]) + GETJSAMPLE(elemptr[3]);
 498
 499     tmp10 = tmp0 + tmp2;
 500     tmp12 = tmp0 - tmp2;
 501
 502     tmp0 = GETJSAMPLE(elemptr[0]) - GETJSAMPLE(elemptr[5]);
 503     tmp1 = GETJSAMPLE(elemptr[1]) - GETJSAMPLE(elemptr[4]);
 504     tmp2 = GETJSAMPLE(elemptr[2]) - GETJSAMPLE(elemptr[3]);
 505
 506     /* Apply unsigned->signed conversion */
 507     dataptr[0] = (DCTELEM)
 508       ((tmp10 + tmp11 - 6 * CENTERJSAMPLE) << PASS1_BITS);
 509     dataptr[2] = (DCTELEM)
 510       DESCALE(MULTIPLY(tmp12, FIX(1.224744871)),                 /* c2 */
 511               CONST_BITS-PASS1_BITS);
 512     dataptr[4] = (DCTELEM)
 513       DESCALE(MULTIPLY(tmp10 - tmp11 - tmp11, FIX(0.707106781)), /* c4 */
 514               CONST_BITS-PASS1_BITS);
 515
 516     /* Odd part */
 517
 518     tmp10 = DESCALE(MULTIPLY(tmp0 + tmp2, FIX(0.366025404)),     /* c5 */
 519                     CONST_BITS-PASS1_BITS);
 520
 521     dataptr[1] = (DCTELEM) (tmp10 + ((tmp0 + tmp1) << PASS1_BITS));
 522     dataptr[3] = (DCTELEM) ((tmp0 - tmp1 - tmp2) << PASS1_BITS);
 523     dataptr[5] = (DCTELEM) (tmp10 + ((tmp2 - tmp1) << PASS1_BITS));
 524
 525     dataptr += DCTSIZE;         /* advance pointer to next row */
 526   }
 527
 528   /* Pass 2: process columns.
 529    * We remove the PASS1_BITS scaling, but leave the results scaled up
 530    * by an overall factor of 8.
 531    * We must also scale the output by (8/6)**2 = 16/9, which we fold
 532    * into the constant multipliers:
 533    * cK now represents sqrt(2) * cos(K*pi/12) * 16/9.
 534    */
 535
 536   dataptr = data;
 537   for (ctr = 0; ctr < 6; ctr++) {
 538     /* Even part */
 539
 540     tmp0 = dataptr[DCTSIZE*0] + dataptr[DCTSIZE*5];
 541     tmp11 = dataptr[DCTSIZE*1] + dataptr[DCTSIZE*4];
 542     tmp2 = dataptr[DCTSIZE*2] + dataptr[DCTSIZE*3];
 543
 544     tmp10 = tmp0 + tmp2;
 545     tmp12 = tmp0 - tmp2;
 546
 547     tmp0 = dataptr[DCTSIZE*0] - dataptr[DCTSIZE*5];
 548     tmp1 = dataptr[DCTSIZE*1] - dataptr[DCTSIZE*4];
 549     tmp2 = dataptr[DCTSIZE*2] - dataptr[DCTSIZE*3];
 550
 551     dataptr[DCTSIZE*0] = (DCTELEM)
 552       DESCALE(MULTIPLY(tmp10 + tmp11, FIX(1.777777778)),         /* 16/9 */
 553               CONST_BITS+PASS1_BITS);
 554     dataptr[DCTSIZE*2] = (DCTELEM)
 555       DESCALE(MULTIPLY(tmp12, FIX(2.177324216)),                 /* c2 */
 556               CONST_BITS+PASS1_BITS);
 557     dataptr[DCTSIZE*4] = (DCTELEM)
 558       DESCALE(MULTIPLY(tmp10 - tmp11 - tmp11, FIX(1.257078722)), /* c4 */
 559               CONST_BITS+PASS1_BITS);
 560
 561     /* Odd part */
 562
 563     tmp10 = MULTIPLY(tmp0 + tmp2, FIX(0.650711829));             /* c5 */
 564
 565     dataptr[DCTSIZE*1] = (DCTELEM)
 566       DESCALE(tmp10 + MULTIPLY(tmp0 + tmp1, FIX(1.777777778)),   /* 16/9 */
 567               CONST_BITS+PASS1_BITS);
 568     dataptr[DCTSIZE*3] = (DCTELEM)
 569       DESCALE(MULTIPLY(tmp0 - tmp1 - tmp2, FIX(1.777777778)),    /* 16/9 */
 570               CONST_BITS+PASS1_BITS);
 571     dataptr[DCTSIZE*5] = (DCTELEM)
 572       DESCALE(tmp10 + MULTIPLY(tmp2 - tmp1, FIX(1.777777778)),   /* 16/9 */
 573               CONST_BITS+PASS1_BITS);
 574
 575     dataptr++;                  /* advance pointer to next column */
 576   }
 577 }
 578
 579
 580 /*
 581  * Perform the forward DCT on a 5x5 sample block.
 582  */
 583
 584 GLOBAL(void)
 585 jpeg_fdct_5x5 (DCTELEM * data, JSAMPARRAY sample_data, JDIMENSION start_col)
 586 {
 587   INT32 tmp0, tmp1, tmp2;
 588   INT32 tmp10, tmp11;
 589   DCTELEM *dataptr;
 590   JSAMPROW elemptr;
 591   int ctr;
 592   SHIFT_TEMPS
 593
 594   /* Pre-zero output coefficient block. */
 595   MEMZERO(data, SIZEOF(DCTELEM) * DCTSIZE2);
 596
 597   /* Pass 1: process rows.
 598    * Note results are scaled up by sqrt(8) compared to a true DCT;
 599    * furthermore, we scale the results by 2**PASS1_BITS.
 600    * We scale the results further by 2 as part of output adaption
 601    * scaling for different DCT size.
 602    * cK represents sqrt(2) * cos(K*pi/10).
 603    */
 604
 605   dataptr = data;
 606   for (ctr = 0; ctr < 5; ctr++) {
 607     elemptr = sample_data[ctr] + start_col;
 608
 609     /* Even part */
 610
 611     tmp0 = GETJSAMPLE(elemptr[0]) + GETJSAMPLE(elemptr[4]);
 612     tmp1 = GETJSAMPLE(elemptr[1]) + GETJSAMPLE(elemptr[3]);
 613     tmp2 = GETJSAMPLE(elemptr[2]);
 614
 615     tmp10 = tmp0 + tmp1;
 616     tmp11 = tmp0 - tmp1;
 617
 618     tmp0 = GETJSAMPLE(elemptr[0]) - GETJSAMPLE(elemptr[4]);
 619     tmp1 = GETJSAMPLE(elemptr[1]) - GETJSAMPLE(elemptr[3]);
 620
 621     /* Apply unsigned->signed conversion */
 622     dataptr[0] = (DCTELEM)
 623       ((tmp10 + tmp2 - 5 * CENTERJSAMPLE) << (PASS1_BITS+1));
 624     tmp11 = MULTIPLY(tmp11, FIX(0.790569415));          /* (c2+c4)/2 */
 625     tmp10 -= tmp2 << 2;
 626     tmp10 = MULTIPLY(tmp10, FIX(0.353553391));          /* (c2-c4)/2 */
 627     dataptr[2] = (DCTELEM) DESCALE(tmp11 + tmp10, CONST_BITS-PASS1_BITS-1);
 628     dataptr[4] = (DCTELEM) DESCALE(tmp11 - tmp10, CONST_BITS-PASS1_BITS-1);
 629
 630     /* Odd part */
 631
 632     tmp10 = MULTIPLY(tmp0 + tmp1, FIX(0.831253876));    /* c3 */
 633
 634     dataptr[1] = (DCTELEM)
 635       DESCALE(tmp10 + MULTIPLY(tmp0, FIX(0.513743148)), /* c1-c3 */
 636               CONST_BITS-PASS1_BITS-1);
 637     dataptr[3] = (DCTELEM)
 638       DESCALE(tmp10 - MULTIPLY(tmp1, FIX(2.176250899)), /* c1+c3 */
 639               CONST_BITS-PASS1_BITS-1);
 640
 641     dataptr += DCTSIZE;         /* advance pointer to next row */
 642   }
 643
 644   /* Pass 2: process columns.
 645    * We remove the PASS1_BITS scaling, but leave the results scaled up
 646    * by an overall factor of 8.
 647    * We must also scale the output by (8/5)**2 = 64/25, which we partially
 648    * fold into the constant multipliers (other part was done in pass 1):
 649    * cK now represents sqrt(2) * cos(K*pi/10) * 32/25.
 650    */
 651
 652   dataptr = data;
 653   for (ctr = 0; ctr < 5; ctr++) {
 654     /* Even part */
 655
 656     tmp0 = dataptr[DCTSIZE*0] + dataptr[DCTSIZE*4];
 657     tmp1 = dataptr[DCTSIZE*1] + dataptr[DCTSIZE*3];
 658     tmp2 = dataptr[DCTSIZE*2];
 659
 660     tmp10 = tmp0 + tmp1;
 661     tmp11 = tmp0 - tmp1;
 662
 663     tmp0 = dataptr[DCTSIZE*0] - dataptr[DCTSIZE*4];
 664     tmp1 = dataptr[DCTSIZE*1] - dataptr[DCTSIZE*3];
 665
 666     dataptr[DCTSIZE*0] = (DCTELEM)
 667       DESCALE(MULTIPLY(tmp10 + tmp2, FIX(1.28)),        /* 32/25 */
 668               CONST_BITS+PASS1_BITS);
 669     tmp11 = MULTIPLY(tmp11, FIX(1.011928851));          /* (c2+c4)/2 */
 670     tmp10 -= tmp2 << 2;
 671     tmp10 = MULTIPLY(tmp10, FIX(0.452548340));          /* (c2-c4)/2 */
 672     dataptr[DCTSIZE*2] = (DCTELEM) DESCALE(tmp11 + tmp10, CONST_BITS+PASS1_BITS);
 673     dataptr[DCTSIZE*4] = (DCTELEM) DESCALE(tmp11 - tmp10, CONST_BITS+PASS1_BITS);
 674
 675     /* Odd part */
 676
 677     tmp10 = MULTIPLY(tmp0 + tmp1, FIX(1.064004961));    /* c3 */
 678
 679     dataptr[DCTSIZE*1] = (DCTELEM)
 680       DESCALE(tmp10 + MULTIPLY(tmp0, FIX(0.657591230)), /* c1-c3 */
 681               CONST_BITS+PASS1_BITS);
 682     dataptr[DCTSIZE*3] = (DCTELEM)
 683       DESCALE(tmp10 - MULTIPLY(tmp1, FIX(2.785601151)), /* c1+c3 */
 684               CONST_BITS+PASS1_BITS);
 685
 686     dataptr++;                  /* advance pointer to next column */
 687   }
 688 }
 689
 690
 691 /*
 692  * Perform the forward DCT on a 4x4 sample block.
 693  */
 694
 695 GLOBAL(void)
 696 jpeg_fdct_4x4 (DCTELEM * data, JSAMPARRAY sample_data, JDIMENSION start_col)
 697 {
 698   INT32 tmp0, tmp1;
 699   INT32 tmp10, tmp11;
 700   DCTELEM *dataptr;
 701   JSAMPROW elemptr;
 702   int ctr;
 703   SHIFT_TEMPS
 704
 705   /* Pre-zero output coefficient block. */
 706   MEMZERO(data, SIZEOF(DCTELEM) * DCTSIZE2);
 707
 708   /* Pass 1: process rows.
 709    * Note results are scaled up by sqrt(8) compared to a true DCT;
 710    * furthermore, we scale the results by 2**PASS1_BITS.
 711    * We must also scale the output by (8/4)**2 = 2**2, which we add here.
 712    * cK represents sqrt(2) * cos(K*pi/16) [refers to 8-point FDCT].
 713    */
 714
 715   dataptr = data;
 716   for (ctr = 0; ctr < 4; ctr++) {
 717     elemptr = sample_data[ctr] + start_col;
 718
 719     /* Even part */
 720
 721     tmp0 = GETJSAMPLE(elemptr[0]) + GETJSAMPLE(elemptr[3]);
 722     tmp1 = GETJSAMPLE(elemptr[1]) + GETJSAMPLE(elemptr[2]);
 723
 724     tmp10 = GETJSAMPLE(elemptr[0]) - GETJSAMPLE(elemptr[3]);
 725     tmp11 = GETJSAMPLE(elemptr[1]) - GETJSAMPLE(elemptr[2]);
 726
 727     /* Apply unsigned->signed conversion */
 728     dataptr[0] = (DCTELEM)
 729       ((tmp0 + tmp1 - 4 * CENTERJSAMPLE) << (PASS1_BITS+2));
 730     dataptr[2] = (DCTELEM) ((tmp0 - tmp1) << (PASS1_BITS+2));
 731
 732     /* Odd part */
 733
 734     tmp0 = MULTIPLY(tmp10 + tmp11, FIX_0_541196100);       /* c6 */
 735     /* Add fudge factor here for final descale. */
 736     tmp0 += ONE << (CONST_BITS-PASS1_BITS-3);
 737
 738     dataptr[1] = (DCTELEM)
 739       RIGHT_SHIFT(tmp0 + MULTIPLY(tmp10, FIX_0_765366865), /* c2-c6 */
 740                   CONST_BITS-PASS1_BITS-2);
 741     dataptr[3] = (DCTELEM)
 742       RIGHT_SHIFT(tmp0 - MULTIPLY(tmp11, FIX_1_847759065), /* c2+c6 */
 743                   CONST_BITS-PASS1_BITS-2);
 744
 745     dataptr += DCTSIZE;         /* advance pointer to next row */
 746   }
 747
 748   /* Pass 2: process columns.
 749    * We remove the PASS1_BITS scaling, but leave the results scaled up
 750    * by an overall factor of 8.
 751    * cK represents sqrt(2) * cos(K*pi/16) [refers to 8-point FDCT].
 752    */
 753
 754   dataptr = data;
 755   for (ctr = 0; ctr < 4; ctr++) {
 756     /* Even part */
 757
 758     /* Add fudge factor here for final descale. */
 759     tmp0 = dataptr[DCTSIZE*0] + dataptr[DCTSIZE*3] + (ONE << (PASS1_BITS-1));
 760     tmp1 = dataptr[DCTSIZE*1] + dataptr[DCTSIZE*2];
 761
 762     tmp10 = dataptr[DCTSIZE*0] - dataptr[DCTSIZE*3];
 763     tmp11 = dataptr[DCTSIZE*1] - dataptr[DCTSIZE*2];
 764
 765     dataptr[DCTSIZE*0] = (DCTELEM) RIGHT_SHIFT(tmp0 + tmp1, PASS1_BITS);
 766     dataptr[DCTSIZE*2] = (DCTELEM) RIGHT_SHIFT(tmp0 - tmp1, PASS1_BITS);
 767
 768     /* Odd part */
 769
 770     tmp0 = MULTIPLY(tmp10 + tmp11, FIX_0_541196100);       /* c6 */
 771     /* Add fudge factor here for final descale. */
 772     tmp0 += ONE << (CONST_BITS+PASS1_BITS-1);
 773
 774     dataptr[DCTSIZE*1] = (DCTELEM)
 775       RIGHT_SHIFT(tmp0 + MULTIPLY(tmp10, FIX_0_765366865), /* c2-c6 */
 776                   CONST_BITS+PASS1_BITS);
 777     dataptr[DCTSIZE*3] = (DCTELEM)
 778       RIGHT_SHIFT(tmp0 - MULTIPLY(tmp11, FIX_1_847759065), /* c2+c6 */
 779                   CONST_BITS+PASS1_BITS);
 780
 781     dataptr++;                  /* advance pointer to next column */
 782   }
 783 }
 784
 785
 786 /*
 787  * Perform the forward DCT on a 3x3 sample block.
 788  */
 789
 790 GLOBAL(void)
 791 jpeg_fdct_3x3 (DCTELEM * data, JSAMPARRAY sample_data, JDIMENSION start_col)
 792 {
 793   INT32 tmp0, tmp1, tmp2;
 794   DCTELEM *dataptr;
 795   JSAMPROW elemptr;
 796   int ctr;
 797   SHIFT_TEMPS
 798
 799   /* Pre-zero output coefficient block. */
 800   MEMZERO(data, SIZEOF(DCTELEM) * DCTSIZE2);
 801
 802   /* Pass 1: process rows.
 803    * Note results are scaled up by sqrt(8) compared to a true DCT;
 804    * furthermore, we scale the results by 2**PASS1_BITS.
 805    * We scale the results further by 2**2 as part of output adaption
 806    * scaling for different DCT size.
 807    * cK represents sqrt(2) * cos(K*pi/6).
 808    */
 809
 810   dataptr = data;
 811   for (ctr = 0; ctr < 3; ctr++) {
 812     elemptr = sample_data[ctr] + start_col;
 813
 814     /* Even part */
 815
 816     tmp0 = GETJSAMPLE(elemptr[0]) + GETJSAMPLE(elemptr[2]);
 817     tmp1 = GETJSAMPLE(elemptr[1]);
 818
 819     tmp2 = GETJSAMPLE(elemptr[0]) - GETJSAMPLE(elemptr[2]);
 820
 821     /* Apply unsigned->signed conversion */
 822     dataptr[0] = (DCTELEM)
 823       ((tmp0 + tmp1 - 3 * CENTERJSAMPLE) << (PASS1_BITS+2));
 824     dataptr[2] = (DCTELEM)
 825       DESCALE(MULTIPLY(tmp0 - tmp1 - tmp1, FIX(0.707106781)), /* c2 */
 826               CONST_BITS-PASS1_BITS-2);
 827
 828     /* Odd part */
 829
 830     dataptr[1] = (DCTELEM)
 831       DESCALE(MULTIPLY(tmp2, FIX(1.224744871)),               /* c1 */
 832               CONST_BITS-PASS1_BITS-2);
 833
 834     dataptr += DCTSIZE;         /* advance pointer to next row */
 835   }
 836
 837   /* Pass 2: process columns.
 838    * We remove the PASS1_BITS scaling, but leave the results scaled up
 839    * by an overall factor of 8.
 840    * We must also scale the output by (8/3)**2 = 64/9, which we partially
 841    * fold into the constant multipliers (other part was done in pass 1):
 842    * cK now represents sqrt(2) * cos(K*pi/6) * 16/9.
 843    */
 844
 845   dataptr = data;
 846   for (ctr = 0; ctr < 3; ctr++) {
 847     /* Even part */
 848
 849     tmp0 = dataptr[DCTSIZE*0] + dataptr[DCTSIZE*2];
 850     tmp1 = dataptr[DCTSIZE*1];
 851
 852     tmp2 = dataptr[DCTSIZE*0] - dataptr[DCTSIZE*2];
 853
 854     dataptr[DCTSIZE*0] = (DCTELEM)
 855       DESCALE(MULTIPLY(tmp0 + tmp1, FIX(1.777777778)),        /* 16/9 */
 856               CONST_BITS+PASS1_BITS);
 857     dataptr[DCTSIZE*2] = (DCTELEM)
 858       DESCALE(MULTIPLY(tmp0 - tmp1 - tmp1, FIX(1.257078722)), /* c2 */
 859               CONST_BITS+PASS1_BITS);
 860
 861     /* Odd part */
 862
 863     dataptr[DCTSIZE*1] = (DCTELEM)
 864       DESCALE(MULTIPLY(tmp2, FIX(2.177324216)),               /* c1 */
 865               CONST_BITS+PASS1_BITS);
 866
 867     dataptr++;                  /* advance pointer to next column */
 868   }
 869 }
 870
 871
 872 /*
 873  * Perform the forward DCT on a 2x2 sample block.
 874  */
 875
 876 GLOBAL(void)
 877 jpeg_fdct_2x2 (DCTELEM * data, JSAMPARRAY sample_data, JDIMENSION start_col)
 878 {
 879   INT32 tmp0, tmp1, tmp2, tmp3;
 880   JSAMPROW elemptr;
 881
 882   /* Pre-zero output coefficient block. */
 883   MEMZERO(data, SIZEOF(DCTELEM) * DCTSIZE2);
 884
 885   /* Pass 1: process rows.
 886    * Note results are scaled up by sqrt(8) compared to a true DCT.
 887    */
 888
 889   /* Row 0 */
 890   elemptr = sample_data[0] + start_col;
 891
 892   tmp0 = GETJSAMPLE(elemptr[0]) + GETJSAMPLE(elemptr[1]);
 893   tmp1 = GETJSAMPLE(elemptr[0]) - GETJSAMPLE(elemptr[1]);
 894
 895   /* Row 1 */
 896   elemptr = sample_data[1] + start_col;
 897
 898   tmp2 = GETJSAMPLE(elemptr[0]) + GETJSAMPLE(elemptr[1]);
 899   tmp3 = GETJSAMPLE(elemptr[0]) - GETJSAMPLE(elemptr[1]);
 900
 901   /* Pass 2: process columns.
 902    * We leave the results scaled up by an overall factor of 8.
 903    * We must also scale the output by (8/2)**2 = 2**4.
 904    */
 905
 906   /* Column 0 */
 907   /* Apply unsigned->signed conversion */
 908   data[DCTSIZE*0] = (DCTELEM) ((tmp0 + tmp2 - 4 * CENTERJSAMPLE) << 4);
 909   data[DCTSIZE*1] = (DCTELEM) ((tmp0 - tmp2) << 4);
 910
 911   /* Column 1 */
 912   data[DCTSIZE*0+1] = (DCTELEM) ((tmp1 + tmp3) << 4);
 913   data[DCTSIZE*1+1] = (DCTELEM) ((tmp1 - tmp3) << 4);
 914 }
 915
 916
 917 /*
 918  * Perform the forward DCT on a 1x1 sample block.
 919  */
 920
 921 GLOBAL(void)
 922 jpeg_fdct_1x1 (DCTELEM * data, JSAMPARRAY sample_data, JDIMENSION start_col)
 923 {
 924   /* Pre-zero output coefficient block. */
 925   MEMZERO(data, SIZEOF(DCTELEM) * DCTSIZE2);
 926
 927   /* We leave the result scaled up by an overall factor of 8. */
 928   /* We must also scale the output by (8/1)**2 = 2**6. */
 929   /* Apply unsigned->signed conversion */
 930   data[0] = (DCTELEM)
 931     ((GETJSAMPLE(sample_data[0][start_col]) - CENTERJSAMPLE) << 6);
 932 }
 933
 934
 935 /*
 936  * Perform the forward DCT on a 9x9 sample block.
 937  */
 938
 939 GLOBAL(void)
 940 jpeg_fdct_9x9 (DCTELEM * data, JSAMPARRAY sample_data, JDIMENSION start_col)
 941 {
 942   INT32 tmp0, tmp1, tmp2, tmp3, tmp4;
 943   INT32 tmp10, tmp11, tmp12, tmp13;
 944   INT32 z1, z2;
 945   DCTELEM workspace[8];
 946   DCTELEM *dataptr;
 947   DCTELEM *wsptr;
 948   JSAMPROW elemptr;
 949   int ctr;
 950   SHIFT_TEMPS
 951
 952   /* Pass 1: process rows.
 953    * Note results are scaled up by sqrt(8) compared to a true DCT;
 954    * we scale the results further by 2 as part of output adaption
 955    * scaling for different DCT size.
 956    * cK represents sqrt(2) * cos(K*pi/18).
 957    */
 958
 959   dataptr = data;
 960   ctr = 0;
 961   for (;;) {
 962     elemptr = sample_data[ctr] + start_col;
 963
 964     /* Even part */
 965
 966     tmp0 = GETJSAMPLE(elemptr[0]) + GETJSAMPLE(elemptr[8]);
 967     tmp1 = GETJSAMPLE(elemptr[1]) + GETJSAMPLE(elemptr[7]);
 968     tmp2 = GETJSAMPLE(elemptr[2]) + GETJSAMPLE(elemptr[6]);
 969     tmp3 = GETJSAMPLE(elemptr[3]) + GETJSAMPLE(elemptr[5]);
 970     tmp4 = GETJSAMPLE(elemptr[4]);
 971
 972     tmp10 = GETJSAMPLE(elemptr[0]) - GETJSAMPLE(elemptr[8]);
 973     tmp11 = GETJSAMPLE(elemptr[1]) - GETJSAMPLE(elemptr[7]);
 974     tmp12 = GETJSAMPLE(elemptr[2]) - GETJSAMPLE(elemptr[6]);
 975     tmp13 = GETJSAMPLE(elemptr[3]) - GETJSAMPLE(elemptr[5]);
 976
 977     z1 = tmp0 + tmp2 + tmp3;
 978     z2 = tmp1 + tmp4;
 979     /* Apply unsigned->signed conversion */
 980     dataptr[0] = (DCTELEM) ((z1 + z2 - 9 * CENTERJSAMPLE) << 1);
 981     dataptr[6] = (DCTELEM)
 982       DESCALE(MULTIPLY(z1 - z2 - z2, FIX(0.707106781)),  /* c6 */
 983               CONST_BITS-1);
 984     z1 = MULTIPLY(tmp0 - tmp2, FIX(1.328926049));        /* c2 */
 985     z2 = MULTIPLY(tmp1 - tmp4 - tmp4, FIX(0.707106781)); /* c6 */
 986     dataptr[2] = (DCTELEM)
 987       DESCALE(MULTIPLY(tmp2 - tmp3, FIX(1.083350441))    /* c4 */
 988               + z1 + z2, CONST_BITS-1);
 989     dataptr[4] = (DCTELEM)
 990       DESCALE(MULTIPLY(tmp3 - tmp0, FIX(0.245575608))    /* c8 */
 991               + z1 - z2, CONST_BITS-1);
 992
 993     /* Odd part */
 994
 995     dataptr[3] = (DCTELEM)
 996       DESCALE(MULTIPLY(tmp10 - tmp12 - tmp13, FIX(1.224744871)), /* c3 */
 997               CONST_BITS-1);
 998
 999     tmp11 = MULTIPLY(tmp11, FIX(1.224744871));        /* c3 */
1000     tmp0 = MULTIPLY(tmp10 + tmp12, FIX(0.909038955)); /* c5 */
1001     tmp1 = MULTIPLY(tmp10 + tmp13, FIX(0.483689525)); /* c7 */
1002
1003     dataptr[1] = (DCTELEM) DESCALE(tmp11 + tmp0 + tmp1, CONST_BITS-1);
1004
1005     tmp2 = MULTIPLY(tmp12 - tmp13, FIX(1.392728481)); /* c1 */
1006
1007     dataptr[5] = (DCTELEM) DESCALE(tmp0 - tmp11 - tmp2, CONST_BITS-1);
1008     dataptr[7] = (DCTELEM) DESCALE(tmp1 - tmp11 + tmp2, CONST_BITS-1);
1009
1010     ctr++;
1011
1012     if (ctr != DCTSIZE) {
1013       if (ctr == 9)
1014         break;                  /* Done. */
1015       dataptr += DCTSIZE;       /* advance pointer to next row */
1016     } else
1017       dataptr = workspace;      /* switch pointer to extended workspace */
1018   }
1019
1020   /* Pass 2: process columns.
1021    * We leave the results scaled up by an overall factor of 8.
1022    * We must also scale the output by (8/9)**2 = 64/81, which we partially
1023    * fold into the constant multipliers and final/initial shifting:
1024    * cK now represents sqrt(2) * cos(K*pi/18) * 128/81.
1025    */
1026
1027   dataptr = data;
1028   wsptr = workspace;
1029   for (ctr = DCTSIZE-1; ctr >= 0; ctr--) {
1030     /* Even part */
1031
1032     tmp0 = dataptr[DCTSIZE*0] + wsptr[DCTSIZE*0];
1033     tmp1 = dataptr[DCTSIZE*1] + dataptr[DCTSIZE*7];
1034     tmp2 = dataptr[DCTSIZE*2] + dataptr[DCTSIZE*6];
1035     tmp3 = dataptr[DCTSIZE*3] + dataptr[DCTSIZE*5];
1036     tmp4 = dataptr[DCTSIZE*4];
1037
1038     tmp10 = dataptr[DCTSIZE*0] - wsptr[DCTSIZE*0];
1039     tmp11 = dataptr[DCTSIZE*1] - dataptr[DCTSIZE*7];
1040     tmp12 = dataptr[DCTSIZE*2] - dataptr[DCTSIZE*6];
1041     tmp13 = dataptr[DCTSIZE*3] - dataptr[DCTSIZE*5];
1042
1043     z1 = tmp0 + tmp2 + tmp3;
1044     z2 = tmp1 + tmp4;
1045     dataptr[DCTSIZE*0] = (DCTELEM)
1046       DESCALE(MULTIPLY(z1 + z2, FIX(1.580246914)),       /* 128/81 */
1047               CONST_BITS+2);
1048     dataptr[DCTSIZE*6] = (DCTELEM)
1049       DESCALE(MULTIPLY(z1 - z2 - z2, FIX(1.117403309)),  /* c6 */
1050               CONST_BITS+2);
1051     z1 = MULTIPLY(tmp0 - tmp2, FIX(2.100031287));        /* c2 */
1052     z2 = MULTIPLY(tmp1 - tmp4 - tmp4, FIX(1.117403309)); /* c6 */
1053     dataptr[DCTSIZE*2] = (DCTELEM)
1054       DESCALE(MULTIPLY(tmp2 - tmp3, FIX(1.711961190))    /* c4 */
1055               + z1 + z2, CONST_BITS+2);
1056     dataptr[DCTSIZE*4] = (DCTELEM)
1057       DESCALE(MULTIPLY(tmp3 - tmp0, FIX(0.388070096))    /* c8 */
1058               + z1 - z2, CONST_BITS+2);
1059
1060     /* Odd part */
1061
1062     dataptr[DCTSIZE*3] = (DCTELEM)
1063       DESCALE(MULTIPLY(tmp10 - tmp12 - tmp13, FIX(1.935399303)), /* c3 */
1064               CONST_BITS+2);
1065
1066     tmp11 = MULTIPLY(tmp11, FIX(1.935399303));        /* c3 */
1067     tmp0 = MULTIPLY(tmp10 + tmp12, FIX(1.436506004)); /* c5 */
1068     tmp1 = MULTIPLY(tmp10 + tmp13, FIX(0.764348879)); /* c7 */
1069
1070     dataptr[DCTSIZE*1] = (DCTELEM)
1071       DESCALE(tmp11 + tmp0 + tmp1, CONST_BITS+2);
1072
1073     tmp2 = MULTIPLY(tmp12 - tmp13, FIX(2.200854883)); /* c1 */
1074
1075     dataptr[DCTSIZE*5] = (DCTELEM)
1076       DESCALE(tmp0 - tmp11 - tmp2, CONST_BITS+2);
1077     dataptr[DCTSIZE*7] = (DCTELEM)
1078       DESCALE(tmp1 - tmp11 + tmp2, CONST_BITS+2);
1079
1080     dataptr++;                  /* advance pointer to next column */
1081     wsptr++;                    /* advance pointer to next column */
1082   }
1083 }
1084
1085
1086 /*
1087  * Perform the forward DCT on a 10x10 sample block.
1088  */
1089
1090 GLOBAL(void)
1091 jpeg_fdct_10x10 (DCTELEM * data, JSAMPARRAY sample_data, JDIMENSION start_col)
1092 {
1093   INT32 tmp0, tmp1, tmp2, tmp3, tmp4;
1094   INT32 tmp10, tmp11, tmp12, tmp13, tmp14;
1095   DCTELEM workspace[8*2];
1096   DCTELEM *dataptr;
1097   DCTELEM *wsptr;
1098   JSAMPROW elemptr;
1099   int ctr;
1100   SHIFT_TEMPS
1101
1102   /* Pass 1: process rows.
1103    * Note results are scaled up by sqrt(8) compared to a true DCT;
1104    * we scale the results further by 2 as part of output adaption
1105    * scaling for different DCT size.
1106    * cK represents sqrt(2) * cos(K*pi/20).
1107    */
1108
1109   dataptr = data;
1110   ctr = 0;
1111   for (;;) {
1112     elemptr = sample_data[ctr] + start_col;
1113
1114     /* Even part */
1115
1116     tmp0 = GETJSAMPLE(elemptr[0]) + GETJSAMPLE(elemptr[9]);
1117     tmp1 = GETJSAMPLE(elemptr[1]) + GETJSAMPLE(elemptr[8]);
1118     tmp12 = GETJSAMPLE(elemptr[2]) + GETJSAMPLE(elemptr[7]);
1119     tmp3 = GETJSAMPLE(elemptr[3]) + GETJSAMPLE(elemptr[6]);
1120     tmp4 = GETJSAMPLE(elemptr[4]) + GETJSAMPLE(elemptr[5]);
1121
1122     tmp10 = tmp0 + tmp4;
1123     tmp13 = tmp0 - tmp4;
1124     tmp11 = tmp1 + tmp3;
1125     tmp14 = tmp1 - tmp3;
1126
1127     tmp0 = GETJSAMPLE(elemptr[0]) - GETJSAMPLE(elemptr[9]);
1128     tmp1 = GETJSAMPLE(elemptr[1]) - GETJSAMPLE(elemptr[8]);
1129     tmp2 = GETJSAMPLE(elemptr[2]) - GETJSAMPLE(elemptr[7]);
1130     tmp3 = GETJSAMPLE(elemptr[3]) - GETJSAMPLE(elemptr[6]);
1131     tmp4 = GETJSAMPLE(elemptr[4]) - GETJSAMPLE(elemptr[5]);
1132
1133     /* Apply unsigned->signed conversion */
1134     dataptr[0] = (DCTELEM)
1135       ((tmp10 + tmp11 + tmp12 - 10 * CENTERJSAMPLE) << 1);
1136     tmp12 += tmp12;
1137     dataptr[4] = (DCTELEM)
1138       DESCALE(MULTIPLY(tmp10 - tmp12, FIX(1.144122806)) - /* c4 */
1139               MULTIPLY(tmp11 - tmp12, FIX(0.437016024)),  /* c8 */
1140               CONST_BITS-1);
1141     tmp10 = MULTIPLY(tmp13 + tmp14, FIX(0.831253876));    /* c6 */
1142     dataptr[2] = (DCTELEM)
1143       DESCALE(tmp10 + MULTIPLY(tmp13, FIX(0.513743148)),  /* c2-c6 */
1144               CONST_BITS-1);
1145     dataptr[6] = (DCTELEM)
1146       DESCALE(tmp10 - MULTIPLY(tmp14, FIX(2.176250899)),  /* c2+c6 */
1147               CONST_BITS-1);
1148
1149     /* Odd part */
1150
1151     tmp10 = tmp0 + tmp4;
1152     tmp11 = tmp1 - tmp3;
1153     dataptr[5] = (DCTELEM) ((tmp10 - tmp11 - tmp2) << 1);
1154     tmp2 <<= CONST_BITS;
1155     dataptr[1] = (DCTELEM)
1156       DESCALE(MULTIPLY(tmp0, FIX(1.396802247)) +          /* c1 */
1157               MULTIPLY(tmp1, FIX(1.260073511)) + tmp2 +   /* c3 */
1158               MULTIPLY(tmp3, FIX(0.642039522)) +          /* c7 */
1159               MULTIPLY(tmp4, FIX(0.221231742)),           /* c9 */
1160               CONST_BITS-1);
1161     tmp12 = MULTIPLY(tmp0 - tmp4, FIX(0.951056516)) -     /* (c3+c7)/2 */
1162             MULTIPLY(tmp1 + tmp3, FIX(0.587785252));      /* (c1-c9)/2 */
1163     tmp13 = MULTIPLY(tmp10 + tmp11, FIX(0.309016994)) +   /* (c3-c7)/2 */
1164             (tmp11 << (CONST_BITS - 1)) - tmp2;
1165     dataptr[3] = (DCTELEM) DESCALE(tmp12 + tmp13, CONST_BITS-1);
1166     dataptr[7] = (DCTELEM) DESCALE(tmp12 - tmp13, CONST_BITS-1);
1167
1168     ctr++;
1169
1170     if (ctr != DCTSIZE) {
1171       if (ctr == 10)
1172         break;                  /* Done. */
1173       dataptr += DCTSIZE;       /* advance pointer to next row */
1174     } else
1175       dataptr = workspace;      /* switch pointer to extended workspace */
1176   }
1177
1178   /* Pass 2: process columns.
1179    * We leave the results scaled up by an overall factor of 8.
1180    * We must also scale the output by (8/10)**2 = 16/25, which we partially
1181    * fold into the constant multipliers and final/initial shifting:
1182    * cK now represents sqrt(2) * cos(K*pi/20) * 32/25.
1183    */
1184
1185   dataptr = data;
1186   wsptr = workspace;
1187   for (ctr = DCTSIZE-1; ctr >= 0; ctr--) {
1188     /* Even part */
1189
1190     tmp0 = dataptr[DCTSIZE*0] + wsptr[DCTSIZE*1];
1191     tmp1 = dataptr[DCTSIZE*1] + wsptr[DCTSIZE*0];
1192     tmp12 = dataptr[DCTSIZE*2] + dataptr[DCTSIZE*7];
1193     tmp3 = dataptr[DCTSIZE*3] + dataptr[DCTSIZE*6];
1194     tmp4 = dataptr[DCTSIZE*4] + dataptr[DCTSIZE*5];
1195
1196     tmp10 = tmp0 + tmp4;
1197     tmp13 = tmp0 - tmp4;
1198     tmp11 = tmp1 + tmp3;
1199     tmp14 = tmp1 - tmp3;
1200
1201     tmp0 = dataptr[DCTSIZE*0] - wsptr[DCTSIZE*1];
1202     tmp1 = dataptr[DCTSIZE*1] - wsptr[DCTSIZE*0];
1203     tmp2 = dataptr[DCTSIZE*2] - dataptr[DCTSIZE*7];
1204     tmp3 = dataptr[DCTSIZE*3] - dataptr[DCTSIZE*6];
1205     tmp4 = dataptr[DCTSIZE*4] - dataptr[DCTSIZE*5];
1206
1207     dataptr[DCTSIZE*0] = (DCTELEM)
1208       DESCALE(MULTIPLY(tmp10 + tmp11 + tmp12, FIX(1.28)), /* 32/25 */
1209               CONST_BITS+2);
1210     tmp12 += tmp12;
1211     dataptr[DCTSIZE*4] = (DCTELEM)
1212       DESCALE(MULTIPLY(tmp10 - tmp12, FIX(1.464477191)) - /* c4 */
1213               MULTIPLY(tmp11 - tmp12, FIX(0.559380511)),  /* c8 */
1214               CONST_BITS+2);
1215     tmp10 = MULTIPLY(tmp13 + tmp14, FIX(1.064004961));    /* c6 */
1216     dataptr[DCTSIZE*2] = (DCTELEM)
1217       DESCALE(tmp10 + MULTIPLY(tmp13, FIX(0.657591230)),  /* c2-c6 */
1218               CONST_BITS+2);
1219     dataptr[DCTSIZE*6] = (DCTELEM)
1220       DESCALE(tmp10 - MULTIPLY(tmp14, FIX(2.785601151)),  /* c2+c6 */
1221               CONST_BITS+2);
1222
1223     /* Odd part */
1224
1225     tmp10 = tmp0 + tmp4;
1226     tmp11 = tmp1 - tmp3;
1227     dataptr[DCTSIZE*5] = (DCTELEM)
1228       DESCALE(MULTIPLY(tmp10 - tmp11 - tmp2, FIX(1.28)),  /* 32/25 */
1229               CONST_BITS+2);
1230     tmp2 = MULTIPLY(tmp2, FIX(1.28));                     /* 32/25 */
1231     dataptr[DCTSIZE*1] = (DCTELEM)
1232       DESCALE(MULTIPLY(tmp0, FIX(1.787906876)) +          /* c1 */
1233               MULTIPLY(tmp1, FIX(1.612894094)) + tmp2 +   /* c3 */
1234               MULTIPLY(tmp3, FIX(0.821810588)) +          /* c7 */
1235               MULTIPLY(tmp4, FIX(0.283176630)),           /* c9 */
1236               CONST_BITS+2);
1237     tmp12 = MULTIPLY(tmp0 - tmp4, FIX(1.217352341)) -     /* (c3+c7)/2 */
1238             MULTIPLY(tmp1 + tmp3, FIX(0.752365123));      /* (c1-c9)/2 */
1239     tmp13 = MULTIPLY(tmp10 + tmp11, FIX(0.395541753)) +   /* (c3-c7)/2 */
1240             MULTIPLY(tmp11, FIX(0.64)) - tmp2;            /* 16/25 */
1241     dataptr[DCTSIZE*3] = (DCTELEM) DESCALE(tmp12 + tmp13, CONST_BITS+2);
1242     dataptr[DCTSIZE*7] = (DCTELEM) DESCALE(tmp12 - tmp13, CONST_BITS+2);
1243
1244     dataptr++;                  /* advance pointer to next column */
1245     wsptr++;                    /* advance pointer to next column */
1246   }
1247 }
1248
1249
1250 /*
1251  * Perform the forward DCT on an 11x11 sample block.
1252  */
1253
1254 GLOBAL(void)
1255 jpeg_fdct_11x11 (DCTELEM * data, JSAMPARRAY sample_data, JDIMENSION start_col)
1256 {
1257   INT32 tmp0, tmp1, tmp2, tmp3, tmp4, tmp5;
1258   INT32 tmp10, tmp11, tmp12, tmp13, tmp14;
1259   INT32 z1, z2, z3;
1260   DCTELEM workspace[8*3];
1261   DCTELEM *dataptr;
1262   DCTELEM *wsptr;
1263   JSAMPROW elemptr;
1264   int ctr;
1265   SHIFT_TEMPS
1266
1267   /* Pass 1: process rows.
1268    * Note results are scaled up by sqrt(8) compared to a true DCT;
1269    * we scale the results further by 2 as part of output adaption
1270    * scaling for different DCT size.
1271    * cK represents sqrt(2) * cos(K*pi/22).
1272    */
1273
1274   dataptr = data;
1275   ctr = 0;
1276   for (;;) {
1277     elemptr = sample_data[ctr] + start_col;
1278
1279     /* Even part */
1280
1281     tmp0 = GETJSAMPLE(elemptr[0]) + GETJSAMPLE(elemptr[10]);
1282     tmp1 = GETJSAMPLE(elemptr[1]) + GETJSAMPLE(elemptr[9]);
1283     tmp2 = GETJSAMPLE(elemptr[2]) + GETJSAMPLE(elemptr[8]);
1284     tmp3 = GETJSAMPLE(elemptr[3]) + GETJSAMPLE(elemptr[7]);
1285     tmp4 = GETJSAMPLE(elemptr[4]) + GETJSAMPLE(elemptr[6]);
1286     tmp5 = GETJSAMPLE(elemptr[5]);
1287
1288     tmp10 = GETJSAMPLE(elemptr[0]) - GETJSAMPLE(elemptr[10]);
1289     tmp11 = GETJSAMPLE(elemptr[1]) - GETJSAMPLE(elemptr[9]);
1290     tmp12 = GETJSAMPLE(elemptr[2]) - GETJSAMPLE(elemptr[8]);
1291     tmp13 = GETJSAMPLE(elemptr[3]) - GETJSAMPLE(elemptr[7]);
1292     tmp14 = GETJSAMPLE(elemptr[4]) - GETJSAMPLE(elemptr[6]);
1293
1294     /* Apply unsigned->signed conversion */
1295     dataptr[0] = (DCTELEM)
1296       ((tmp0 + tmp1 + tmp2 + tmp3 + tmp4 + tmp5 - 11 * CENTERJSAMPLE) << 1);
1297     tmp5 += tmp5;
1298     tmp0 -= tmp5;
1299     tmp1 -= tmp5;
1300     tmp2 -= tmp5;
1301     tmp3 -= tmp5;
1302     tmp4 -= tmp5;
1303     z1 = MULTIPLY(tmp0 + tmp3, FIX(1.356927976)) +       /* c2 */
1304          MULTIPLY(tmp2 + tmp4, FIX(0.201263574));        /* c10 */
1305     z2 = MULTIPLY(tmp1 - tmp3, FIX(0.926112931));        /* c6 */
1306     z3 = MULTIPLY(tmp0 - tmp1, FIX(1.189712156));        /* c4 */
1307     dataptr[2] = (DCTELEM)
1308       DESCALE(z1 + z2 - MULTIPLY(tmp3, FIX(1.018300590)) /* c2+c8-c6 */
1309               - MULTIPLY(tmp4, FIX(1.390975730)),        /* c4+c10 */
1310               CONST_BITS-1);
1311     dataptr[4] = (DCTELEM)
1312       DESCALE(z2 + z3 + MULTIPLY(tmp1, FIX(0.062335650)) /* c4-c6-c10 */
1313               - MULTIPLY(tmp2, FIX(1.356927976))         /* c2 */
1314               + MULTIPLY(tmp4, FIX(0.587485545)),        /* c8 */
1315               CONST_BITS-1);
1316     dataptr[6] = (DCTELEM)
1317       DESCALE(z1 + z3 - MULTIPLY(tmp0, FIX(1.620527200)) /* c2+c4-c6 */
1318               - MULTIPLY(tmp2, FIX(0.788749120)),        /* c8+c10 */
1319               CONST_BITS-1);
1320
1321     /* Odd part */
1322
1323     tmp1 = MULTIPLY(tmp10 + tmp11, FIX(1.286413905));    /* c3 */
1324     tmp2 = MULTIPLY(tmp10 + tmp12, FIX(1.068791298));    /* c5 */
1325     tmp3 = MULTIPLY(tmp10 + tmp13, FIX(0.764581576));    /* c7 */
1326     tmp0 = tmp1 + tmp2 + tmp3 - MULTIPLY(tmp10, FIX(1.719967871)) /* c7+c5+c3-c1 */
1327            + MULTIPLY(tmp14, FIX(0.398430003));          /* c9 */
1328     tmp4 = MULTIPLY(tmp11 + tmp12, - FIX(0.764581576));  /* -c7 */
1329     tmp5 = MULTIPLY(tmp11 + tmp13, - FIX(1.399818907));  /* -c1 */
1330     tmp1 += tmp4 + tmp5 + MULTIPLY(tmp11, FIX(1.276416582)) /* c9+c7+c1-c3 */
1331             - MULTIPLY(tmp14, FIX(1.068791298));         /* c5 */
1332     tmp10 = MULTIPLY(tmp12 + tmp13, FIX(0.398430003));   /* c9 */
1333     tmp2 += tmp4 + tmp10 - MULTIPLY(tmp12, FIX(1.989053629)) /* c9+c5+c3-c7 */
1334             + MULTIPLY(tmp14, FIX(1.399818907));         /* c1 */
1335     tmp3 += tmp5 + tmp10 + MULTIPLY(tmp13, FIX(1.305598626)) /* c1+c5-c9-c7 */
1336             - MULTIPLY(tmp14, FIX(1.286413905));         /* c3 */
1337
1338     dataptr[1] = (DCTELEM) DESCALE(tmp0, CONST_BITS-1);
1339     dataptr[3] = (DCTELEM) DESCALE(tmp1, CONST_BITS-1);
1340     dataptr[5] = (DCTELEM) DESCALE(tmp2, CONST_BITS-1);
1341     dataptr[7] = (DCTELEM) DESCALE(tmp3, CONST_BITS-1);
1342
1343     ctr++;
1344
1345     if (ctr != DCTSIZE) {
1346       if (ctr == 11)
1347         break;                  /* Done. */
1348       dataptr += DCTSIZE;       /* advance pointer to next row */
1349     } else
1350       dataptr = workspace;      /* switch pointer to extended workspace */
1351   }
1352
1353   /* Pass 2: process columns.
1354    * We leave the results scaled up by an overall factor of 8.
1355    * We must also scale the output by (8/11)**2 = 64/121, which we partially
1356    * fold into the constant multipliers and final/initial shifting:
1357    * cK now represents sqrt(2) * cos(K*pi/22) * 128/121.
1358    */
1359
1360   dataptr = data;
1361   wsptr = workspace;
1362   for (ctr = DCTSIZE-1; ctr >= 0; ctr--) {
1363     /* Even part */
1364
1365     tmp0 = dataptr[DCTSIZE*0] + wsptr[DCTSIZE*2];
1366     tmp1 = dataptr[DCTSIZE*1] + wsptr[DCTSIZE*1];
1367     tmp2 = dataptr[DCTSIZE*2] + wsptr[DCTSIZE*0];
1368     tmp3 = dataptr[DCTSIZE*3] + dataptr[DCTSIZE*7];
1369     tmp4 = dataptr[DCTSIZE*4] + dataptr[DCTSIZE*6];
1370     tmp5 = dataptr[DCTSIZE*5];
1371
1372     tmp10 = dataptr[DCTSIZE*0] - wsptr[DCTSIZE*2];
1373     tmp11 = dataptr[DCTSIZE*1] - wsptr[DCTSIZE*1];
1374     tmp12 = dataptr[DCTSIZE*2] - wsptr[DCTSIZE*0];
1375     tmp13 = dataptr[DCTSIZE*3] - dataptr[DCTSIZE*7];
1376     tmp14 = dataptr[DCTSIZE*4] - dataptr[DCTSIZE*6];
1377
1378     dataptr[DCTSIZE*0] = (DCTELEM)
1379       DESCALE(MULTIPLY(tmp0 + tmp1 + tmp2 + tmp3 + tmp4 + tmp5,
1380                        FIX(1.057851240)),                /* 128/121 */
1381               CONST_BITS+2);
1382     tmp5 += tmp5;
1383     tmp0 -= tmp5;
1384     tmp1 -= tmp5;
1385     tmp2 -= tmp5;
1386     tmp3 -= tmp5;
1387     tmp4 -= tmp5;
1388     z1 = MULTIPLY(tmp0 + tmp3, FIX(1.435427942)) +       /* c2 */
1389          MULTIPLY(tmp2 + tmp4, FIX(0.212906922));        /* c10 */
1390     z2 = MULTIPLY(tmp1 - tmp3, FIX(0.979689713));        /* c6 */
1391     z3 = MULTIPLY(tmp0 - tmp1, FIX(1.258538479));        /* c4 */
1392     dataptr[DCTSIZE*2] = (DCTELEM)
1393       DESCALE(z1 + z2 - MULTIPLY(tmp3, FIX(1.077210542)) /* c2+c8-c6 */
1394               - MULTIPLY(tmp4, FIX(1.471445400)),        /* c4+c10 */
1395               CONST_BITS+2);
1396     dataptr[DCTSIZE*4] = (DCTELEM)
1397       DESCALE(z2 + z3 + MULTIPLY(tmp1, FIX(0.065941844)) /* c4-c6-c10 */
1398               - MULTIPLY(tmp2, FIX(1.435427942))         /* c2 */
1399               + MULTIPLY(tmp4, FIX(0.621472312)),        /* c8 */
1400               CONST_BITS+2);
1401     dataptr[DCTSIZE*6] = (DCTELEM)
1402       DESCALE(z1 + z3 - MULTIPLY(tmp0, FIX(1.714276708)) /* c2+c4-c6 */
1403               - MULTIPLY(tmp2, FIX(0.834379234)),        /* c8+c10 */
1404               CONST_BITS+2);
1405
1406     /* Odd part */
1407
1408     tmp1 = MULTIPLY(tmp10 + tmp11, FIX(1.360834544));    /* c3 */
1409     tmp2 = MULTIPLY(tmp10 + tmp12, FIX(1.130622199));    /* c5 */
1410     tmp3 = MULTIPLY(tmp10 + tmp13, FIX(0.808813568));    /* c7 */
1411     tmp0 = tmp1 + tmp2 + tmp3 - MULTIPLY(tmp10, FIX(1.819470145)) /* c7+c5+c3-c1 */
1412            + MULTIPLY(tmp14, FIX(0.421479672));          /* c9 */
1413     tmp4 = MULTIPLY(tmp11 + tmp12, - FIX(0.808813568));  /* -c7 */
1414     tmp5 = MULTIPLY(tmp11 + tmp13, - FIX(1.480800167));  /* -c1 */
1415     tmp1 += tmp4 + tmp5 + MULTIPLY(tmp11, FIX(1.350258864)) /* c9+c7+c1-c3 */
1416             - MULTIPLY(tmp14, FIX(1.130622199));         /* c5 */
1417     tmp10 = MULTIPLY(tmp12 + tmp13, FIX(0.421479672));   /* c9 */
1418     tmp2 += tmp4 + tmp10 - MULTIPLY(tmp12, FIX(2.104122847)) /* c9+c5+c3-c7 */
1419             + MULTIPLY(tmp14, FIX(1.480800167));         /* c1 */
1420     tmp3 += tmp5 + tmp10 + MULTIPLY(tmp13, FIX(1.381129125)) /* c1+c5-c9-c7 */
1421             - MULTIPLY(tmp14, FIX(1.360834544));         /* c3 */
1422
1423     dataptr[DCTSIZE*1] = (DCTELEM) DESCALE(tmp0, CONST_BITS+2);
1424     dataptr[DCTSIZE*3] = (DCTELEM) DESCALE(tmp1, CONST_BITS+2);
1425     dataptr[DCTSIZE*5] = (DCTELEM) DESCALE(tmp2, CONST_BITS+2);
1426     dataptr[DCTSIZE*7] = (DCTELEM) DESCALE(tmp3, CONST_BITS+2);
1427
1428     dataptr++;                  /* advance pointer to next column */
1429     wsptr++;                    /* advance pointer to next column */
1430   }
1431 }
1432
1433
1434 /*
1435  * Perform the forward DCT on a 12x12 sample block.
1436  */
1437
1438 GLOBAL(void)
1439 jpeg_fdct_12x12 (DCTELEM * data, JSAMPARRAY sample_data, JDIMENSION start_col)
1440 {
1441   INT32 tmp0, tmp1, tmp2, tmp3, tmp4, tmp5;
1442   INT32 tmp10, tmp11, tmp12, tmp13, tmp14, tmp15;
1443   DCTELEM workspace[8*4];
1444   DCTELEM *dataptr;
1445   DCTELEM *wsptr;
1446   JSAMPROW elemptr;
1447   int ctr;
1448   SHIFT_TEMPS
1449
1450   /* Pass 1: process rows.
1451    * Note results are scaled up by sqrt(8) compared to a true DCT.
1452    * cK represents sqrt(2) * cos(K*pi/24).
1453    */
1454
1455   dataptr = data;
1456   ctr = 0;
1457   for (;;) {
1458     elemptr = sample_data[ctr] + start_col;
1459
1460     /* Even part */
1461
1462     tmp0 = GETJSAMPLE(elemptr[0]) + GETJSAMPLE(elemptr[11]);
1463     tmp1 = GETJSAMPLE(elemptr[1]) + GETJSAMPLE(elemptr[10]);
1464     tmp2 = GETJSAMPLE(elemptr[2]) + GETJSAMPLE(elemptr[9]);
1465     tmp3 = GETJSAMPLE(elemptr[3]) + GETJSAMPLE(elemptr[8]);
1466     tmp4 = GETJSAMPLE(elemptr[4]) + GETJSAMPLE(elemptr[7]);
1467     tmp5 = GETJSAMPLE(elemptr[5]) + GETJSAMPLE(elemptr[6]);
1468
1469     tmp10 = tmp0 + tmp5;
1470     tmp13 = tmp0 - tmp5;
1471     tmp11 = tmp1 + tmp4;
1472     tmp14 = tmp1 - tmp4;
1473     tmp12 = tmp2 + tmp3;
1474     tmp15 = tmp2 - tmp3;
1475
1476     tmp0 = GETJSAMPLE(elemptr[0]) - GETJSAMPLE(elemptr[11]);
1477     tmp1 = GETJSAMPLE(elemptr[1]) - GETJSAMPLE(elemptr[10]);
1478     tmp2 = GETJSAMPLE(elemptr[2]) - GETJSAMPLE(elemptr[9]);
1479     tmp3 = GETJSAMPLE(elemptr[3]) - GETJSAMPLE(elemptr[8]);
1480     tmp4 = GETJSAMPLE(elemptr[4]) - GETJSAMPLE(elemptr[7]);
1481     tmp5 = GETJSAMPLE(elemptr[5]) - GETJSAMPLE(elemptr[6]);
1482
1483     /* Apply unsigned->signed conversion */
1484     dataptr[0] = (DCTELEM) (tmp10 + tmp11 + tmp12 - 12 * CENTERJSAMPLE);
1485     dataptr[6] = (DCTELEM) (tmp13 - tmp14 - tmp15);
1486     dataptr[4] = (DCTELEM)
1487       DESCALE(MULTIPLY(tmp10 - tmp12, FIX(1.224744871)), /* c4 */
1488               CONST_BITS);
1489     dataptr[2] = (DCTELEM)
1490       DESCALE(tmp14 - tmp15 + MULTIPLY(tmp13 + tmp15, FIX(1.366025404)), /* c2 */
1491               CONST_BITS);
1492
1493     /* Odd part */
1494
1495     tmp10 = MULTIPLY(tmp1 + tmp4, FIX_0_541196100);    /* c9 */
1496     tmp14 = tmp10 + MULTIPLY(tmp1, FIX_0_765366865);   /* c3-c9 */
1497     tmp15 = tmp10 - MULTIPLY(tmp4, FIX_1_847759065);   /* c3+c9 */
1498     tmp12 = MULTIPLY(tmp0 + tmp2, FIX(1.121971054));   /* c5 */
1499     tmp13 = MULTIPLY(tmp0 + tmp3, FIX(0.860918669));   /* c7 */
1500     tmp10 = tmp12 + tmp13 + tmp14 - MULTIPLY(tmp0, FIX(0.580774953)) /* c5+c7-c1 */
1501             + MULTIPLY(tmp5, FIX(0.184591911));        /* c11 */
1502     tmp11 = MULTIPLY(tmp2 + tmp3, - FIX(0.184591911)); /* -c11 */
1503     tmp12 += tmp11 - tmp15 - MULTIPLY(tmp2, FIX(2.339493912)) /* c1+c5-c11 */
1504             + MULTIPLY(tmp5, FIX(0.860918669));        /* c7 */
1505     tmp13 += tmp11 - tmp14 + MULTIPLY(tmp3, FIX(0.725788011)) /* c1+c11-c7 */
1506             - MULTIPLY(tmp5, FIX(1.121971054));        /* c5 */
1507     tmp11 = tmp15 + MULTIPLY(tmp0 - tmp3, FIX(1.306562965)) /* c3 */
1508             - MULTIPLY(tmp2 + tmp5, FIX_0_541196100);  /* c9 */
1509
1510     dataptr[1] = (DCTELEM) DESCALE(tmp10, CONST_BITS);
1511     dataptr[3] = (DCTELEM) DESCALE(tmp11, CONST_BITS);
1512     dataptr[5] = (DCTELEM) DESCALE(tmp12, CONST_BITS);
1513     dataptr[7] = (DCTELEM) DESCALE(tmp13, CONST_BITS);
1514
1515     ctr++;
1516
1517     if (ctr != DCTSIZE) {
1518       if (ctr == 12)
1519         break;                  /* Done. */
1520       dataptr += DCTSIZE;       /* advance pointer to next row */
1521     } else
1522       dataptr = workspace;      /* switch pointer to extended workspace */
1523   }
1524
1525   /* Pass 2: process columns.
1526    * We leave the results scaled up by an overall factor of 8.
1527    * We must also scale the output by (8/12)**2 = 4/9, which we partially
1528    * fold into the constant multipliers and final shifting:
1529    * cK now represents sqrt(2) * cos(K*pi/24) * 8/9.
1530    */
1531
1532   dataptr = data;
1533   wsptr = workspace;
1534   for (ctr = DCTSIZE-1; ctr >= 0; ctr--) {
1535     /* Even part */
1536
1537     tmp0 = dataptr[DCTSIZE*0] + wsptr[DCTSIZE*3];
1538     tmp1 = dataptr[DCTSIZE*1] + wsptr[DCTSIZE*2];
1539     tmp2 = dataptr[DCTSIZE*2] + wsptr[DCTSIZE*1];
1540     tmp3 = dataptr[DCTSIZE*3] + wsptr[DCTSIZE*0];
1541     tmp4 = dataptr[DCTSIZE*4] + dataptr[DCTSIZE*7];
1542     tmp5 = dataptr[DCTSIZE*5] + dataptr[DCTSIZE*6];
1543
1544     tmp10 = tmp0 + tmp5;
1545     tmp13 = tmp0 - tmp5;
1546     tmp11 = tmp1 + tmp4;
1547     tmp14 = tmp1 - tmp4;
1548     tmp12 = tmp2 + tmp3;
1549     tmp15 = tmp2 - tmp3;
1550
1551     tmp0 = dataptr[DCTSIZE*0] - wsptr[DCTSIZE*3];
1552     tmp1 = dataptr[DCTSIZE*1] - wsptr[DCTSIZE*2];
1553     tmp2 = dataptr[DCTSIZE*2] - wsptr[DCTSIZE*1];
1554     tmp3 = dataptr[DCTSIZE*3] - wsptr[DCTSIZE*0];
1555     tmp4 = dataptr[DCTSIZE*4] - dataptr[DCTSIZE*7];
1556     tmp5 = dataptr[DCTSIZE*5] - dataptr[DCTSIZE*6];
1557
1558     dataptr[DCTSIZE*0] = (DCTELEM)
1559       DESCALE(MULTIPLY(tmp10 + tmp11 + tmp12, FIX(0.888888889)), /* 8/9 */
1560               CONST_BITS+1);
1561     dataptr[DCTSIZE*6] = (DCTELEM)
1562       DESCALE(MULTIPLY(tmp13 - tmp14 - tmp15, FIX(0.888888889)), /* 8/9 */
1563               CONST_BITS+1);
1564     dataptr[DCTSIZE*4] = (DCTELEM)
1565       DESCALE(MULTIPLY(tmp10 - tmp12, FIX(1.088662108)),         /* c4 */
1566               CONST_BITS+1);
1567     dataptr[DCTSIZE*2] = (DCTELEM)
1568       DESCALE(MULTIPLY(tmp14 - tmp15, FIX(0.888888889)) +        /* 8/9 */
1569               MULTIPLY(tmp13 + tmp15, FIX(1.214244803)),         /* c2 */
1570               CONST_BITS+1);
1571
1572     /* Odd part */
1573
1574     tmp10 = MULTIPLY(tmp1 + tmp4, FIX(0.481063200));   /* c9 */
1575     tmp14 = tmp10 + MULTIPLY(tmp1, FIX(0.680326102));  /* c3-c9 */
1576     tmp15 = tmp10 - MULTIPLY(tmp4, FIX(1.642452502));  /* c3+c9 */
1577     tmp12 = MULTIPLY(tmp0 + tmp2, FIX(0.997307603));   /* c5 */
1578     tmp13 = MULTIPLY(tmp0 + tmp3, FIX(0.765261039));   /* c7 */
1579     tmp10 = tmp12 + tmp13 + tmp14 - MULTIPLY(tmp0, FIX(0.516244403)) /* c5+c7-c1 */
1580             + MULTIPLY(tmp5, FIX(0.164081699));        /* c11 */
1581     tmp11 = MULTIPLY(tmp2 + tmp3, - FIX(0.164081699)); /* -c11 */
1582     tmp12 += tmp11 - tmp15 - MULTIPLY(tmp2, FIX(2.079550144)) /* c1+c5-c11 */
1583             + MULTIPLY(tmp5, FIX(0.765261039));        /* c7 */
1584     tmp13 += tmp11 - tmp14 + MULTIPLY(tmp3, FIX(0.645144899)) /* c1+c11-c7 */
1585             - MULTIPLY(tmp5, FIX(0.997307603));        /* c5 */
1586     tmp11 = tmp15 + MULTIPLY(tmp0 - tmp3, FIX(1.161389302)) /* c3 */
1587             - MULTIPLY(tmp2 + tmp5, FIX(0.481063200)); /* c9 */
1588
1589     dataptr[DCTSIZE*1] = (DCTELEM) DESCALE(tmp10, CONST_BITS+1);
1590     dataptr[DCTSIZE*3] = (DCTELEM) DESCALE(tmp11, CONST_BITS+1);
1591     dataptr[DCTSIZE*5] = (DCTELEM) DESCALE(tmp12, CONST_BITS+1);
1592     dataptr[DCTSIZE*7] = (DCTELEM) DESCALE(tmp13, CONST_BITS+1);
1593
1594     dataptr++;                  /* advance pointer to next column */
1595     wsptr++;                    /* advance pointer to next column */
1596   }
1597 }
1598
1599
1600 /*
1601  * Perform the forward DCT on a 13x13 sample block.
1602  */
1603
1604 GLOBAL(void)
1605 jpeg_fdct_13x13 (DCTELEM * data, JSAMPARRAY sample_data, JDIMENSION start_col)
1606 {
1607   INT32 tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6;
1608   INT32 tmp10, tmp11, tmp12, tmp13, tmp14, tmp15;
1609   INT32 z1, z2;
1610   DCTELEM workspace[8*5];
1611   DCTELEM *dataptr;
1612   DCTELEM *wsptr;
1613   JSAMPROW elemptr;
1614   int ctr;
1615   SHIFT_TEMPS
1616
1617   /* Pass 1: process rows.
1618    * Note results are scaled up by sqrt(8) compared to a true DCT.
1619    * cK represents sqrt(2) * cos(K*pi/26).
1620    */
1621
1622   dataptr = data;
1623   ctr = 0;
1624   for (;;) {
1625     elemptr = sample_data[ctr] + start_col;
1626
1627     /* Even part */
1628
1629     tmp0 = GETJSAMPLE(elemptr[0]) + GETJSAMPLE(elemptr[12]);
1630     tmp1 = GETJSAMPLE(elemptr[1]) + GETJSAMPLE(elemptr[11]);
1631     tmp2 = GETJSAMPLE(elemptr[2]) + GETJSAMPLE(elemptr[10]);
1632     tmp3 = GETJSAMPLE(elemptr[3]) + GETJSAMPLE(elemptr[9]);
1633     tmp4 = GETJSAMPLE(elemptr[4]) + GETJSAMPLE(elemptr[8]);
1634     tmp5 = GETJSAMPLE(elemptr[5]) + GETJSAMPLE(elemptr[7]);
1635     tmp6 = GETJSAMPLE(elemptr[6]);
1636
1637     tmp10 = GETJSAMPLE(elemptr[0]) - GETJSAMPLE(elemptr[12]);
1638     tmp11 = GETJSAMPLE(elemptr[1]) - GETJSAMPLE(elemptr[11]);
1639     tmp12 = GETJSAMPLE(elemptr[2]) - GETJSAMPLE(elemptr[10]);
1640     tmp13 = GETJSAMPLE(elemptr[3]) - GETJSAMPLE(elemptr[9]);
1641     tmp14 = GETJSAMPLE(elemptr[4]) - GETJSAMPLE(elemptr[8]);
1642     tmp15 = GETJSAMPLE(elemptr[5]) - GETJSAMPLE(elemptr[7]);
1643
1644     /* Apply unsigned->signed conversion */
1645     dataptr[0] = (DCTELEM)
1646       (tmp0 + tmp1 + tmp2 + tmp3 + tmp4 + tmp5 + tmp6 - 13 * CENTERJSAMPLE);
1647     tmp6 += tmp6;
1648     tmp0 -= tmp6;
1649     tmp1 -= tmp6;
1650     tmp2 -= tmp6;
1651     tmp3 -= tmp6;
1652     tmp4 -= tmp6;
1653     tmp5 -= tmp6;
1654     dataptr[2] = (DCTELEM)
1655       DESCALE(MULTIPLY(tmp0, FIX(1.373119086)) +   /* c2 */
1656               MULTIPLY(tmp1, FIX(1.058554052)) +   /* c6 */
1657               MULTIPLY(tmp2, FIX(0.501487041)) -   /* c10 */
1658               MULTIPLY(tmp3, FIX(0.170464608)) -   /* c12 */
1659               MULTIPLY(tmp4, FIX(0.803364869)) -   /* c8 */
1660               MULTIPLY(tmp5, FIX(1.252223920)),    /* c4 */
1661               CONST_BITS);
1662     z1 = MULTIPLY(tmp0 - tmp2, FIX(1.155388986)) - /* (c4+c6)/2 */
1663          MULTIPLY(tmp3 - tmp4, FIX(0.435816023)) - /* (c2-c10)/2 */
1664          MULTIPLY(tmp1 - tmp5, FIX(0.316450131));  /* (c8-c12)/2 */
1665     z2 = MULTIPLY(tmp0 + tmp2, FIX(0.096834934)) - /* (c4-c6)/2 */
1666          MULTIPLY(tmp3 + tmp4, FIX(0.937303064)) + /* (c2+c10)/2 */
1667          MULTIPLY(tmp1 + tmp5, FIX(0.486914739));  /* (c8+c12)/2 */
1668
1669     dataptr[4] = (DCTELEM) DESCALE(z1 + z2, CONST_BITS);
1670     dataptr[6] = (DCTELEM) DESCALE(z1 - z2, CONST_BITS);
1671
1672     /* Odd part */
1673
1674     tmp1 = MULTIPLY(tmp10 + tmp11, FIX(1.322312651));   /* c3 */
1675     tmp2 = MULTIPLY(tmp10 + tmp12, FIX(1.163874945));   /* c5 */
1676     tmp3 = MULTIPLY(tmp10 + tmp13, FIX(0.937797057)) +  /* c7 */
1677            MULTIPLY(tmp14 + tmp15, FIX(0.338443458));   /* c11 */
1678     tmp0 = tmp1 + tmp2 + tmp3 -
1679            MULTIPLY(tmp10, FIX(2.020082300)) +          /* c3+c5+c7-c1 */
1680            MULTIPLY(tmp14, FIX(0.318774355));           /* c9-c11 */
1681     tmp4 = MULTIPLY(tmp14 - tmp15, FIX(0.937797057)) -  /* c7 */
1682            MULTIPLY(tmp11 + tmp12, FIX(0.338443458));   /* c11 */
1683     tmp5 = MULTIPLY(tmp11 + tmp13, - FIX(1.163874945)); /* -c5 */
1684     tmp1 += tmp4 + tmp5 +
1685             MULTIPLY(tmp11, FIX(0.837223564)) -         /* c5+c9+c11-c3 */
1686             MULTIPLY(tmp14, FIX(2.341699410));          /* c1+c7 */
1687     tmp6 = MULTIPLY(tmp12 + tmp13, - FIX(0.657217813)); /* -c9 */
1688     tmp2 += tmp4 + tmp6 -
1689             MULTIPLY(tmp12, FIX(1.572116027)) +         /* c1+c5-c9-c11 */
1690             MULTIPLY(tmp15, FIX(2.260109708));          /* c3+c7 */
1691     tmp3 += tmp5 + tmp6 +
1692             MULTIPLY(tmp13, FIX(2.205608352)) -         /* c3+c5+c9-c7 */
1693             MULTIPLY(tmp15, FIX(1.742345811));          /* c1+c11 */
1694
1695     dataptr[1] = (DCTELEM) DESCALE(tmp0, CONST_BITS);
1696     dataptr[3] = (DCTELEM) DESCALE(tmp1, CONST_BITS);
1697     dataptr[5] = (DCTELEM) DESCALE(tmp2, CONST_BITS);
1698     dataptr[7] = (DCTELEM) DESCALE(tmp3, CONST_BITS);
1699
1700     ctr++;
1701
1702     if (ctr != DCTSIZE) {
1703       if (ctr == 13)
1704         break;                  /* Done. */
1705       dataptr += DCTSIZE;       /* advance pointer to next row */
1706     } else
1707       dataptr = workspace;      /* switch pointer to extended workspace */
1708   }
1709
1710   /* Pass 2: process columns.
1711    * We leave the results scaled up by an overall factor of 8.
1712    * We must also scale the output by (8/13)**2 = 64/169, which we partially
1713    * fold into the constant multipliers and final shifting:
1714    * cK now represents sqrt(2) * cos(K*pi/26) * 128/169.
1715    */
1716
1717   dataptr = data;
1718   wsptr = workspace;
1719   for (ctr = DCTSIZE-1; ctr >= 0; ctr--) {
1720     /* Even part */
1721
1722     tmp0 = dataptr[DCTSIZE*0] + wsptr[DCTSIZE*4];
1723     tmp1 = dataptr[DCTSIZE*1] + wsptr[DCTSIZE*3];
1724     tmp2 = dataptr[DCTSIZE*2] + wsptr[DCTSIZE*2];
1725     tmp3 = dataptr[DCTSIZE*3] + wsptr[DCTSIZE*1];
1726     tmp4 = dataptr[DCTSIZE*4] + wsptr[DCTSIZE*0];
1727     tmp5 = dataptr[DCTSIZE*5] + dataptr[DCTSIZE*7];
1728     tmp6 = dataptr[DCTSIZE*6];
1729
1730     tmp10 = dataptr[DCTSIZE*0] - wsptr[DCTSIZE*4];
1731     tmp11 = dataptr[DCTSIZE*1] - wsptr[DCTSIZE*3];
1732     tmp12 = dataptr[DCTSIZE*2] - wsptr[DCTSIZE*2];
1733     tmp13 = dataptr[DCTSIZE*3] - wsptr[DCTSIZE*1];
1734     tmp14 = dataptr[DCTSIZE*4] - wsptr[DCTSIZE*0];
1735     tmp15 = dataptr[DCTSIZE*5] - dataptr[DCTSIZE*7];
1736
1737     dataptr[DCTSIZE*0] = (DCTELEM)
1738       DESCALE(MULTIPLY(tmp0 + tmp1 + tmp2 + tmp3 + tmp4 + tmp5 + tmp6,
1739                        FIX(0.757396450)),          /* 128/169 */
1740               CONST_BITS+1);
1741     tmp6 += tmp6;
1742     tmp0 -= tmp6;
1743     tmp1 -= tmp6;
1744     tmp2 -= tmp6;
1745     tmp3 -= tmp6;
1746     tmp4 -= tmp6;
1747     tmp5 -= tmp6;
1748     dataptr[DCTSIZE*2] = (DCTELEM)
1749       DESCALE(MULTIPLY(tmp0, FIX(1.039995521)) +   /* c2 */
1750               MULTIPLY(tmp1, FIX(0.801745081)) +   /* c6 */
1751               MULTIPLY(tmp2, FIX(0.379824504)) -   /* c10 */
1752               MULTIPLY(tmp3, FIX(0.129109289)) -   /* c12 */
1753               MULTIPLY(tmp4, FIX(0.608465700)) -   /* c8 */
1754               MULTIPLY(tmp5, FIX(0.948429952)),    /* c4 */
1755               CONST_BITS+1);
1756     z1 = MULTIPLY(tmp0 - tmp2, FIX(0.875087516)) - /* (c4+c6)/2 */
1757          MULTIPLY(tmp3 - tmp4, FIX(0.330085509)) - /* (c2-c10)/2 */
1758          MULTIPLY(tmp1 - tmp5, FIX(0.239678205));  /* (c8-c12)/2 */
1759     z2 = MULTIPLY(tmp0 + tmp2, FIX(0.073342435)) - /* (c4-c6)/2 */
1760          MULTIPLY(tmp3 + tmp4, FIX(0.709910013)) + /* (c2+c10)/2 */
1761          MULTIPLY(tmp1 + tmp5, FIX(0.368787494));  /* (c8+c12)/2 */
1762
1763     dataptr[DCTSIZE*4] = (DCTELEM) DESCALE(z1 + z2, CONST_BITS+1);
1764     dataptr[DCTSIZE*6] = (DCTELEM) DESCALE(z1 - z2, CONST_BITS+1);
1765
1766     /* Odd part */
1767
1768     tmp1 = MULTIPLY(tmp10 + tmp11, FIX(1.001514908));   /* c3 */
1769     tmp2 = MULTIPLY(tmp10 + tmp12, FIX(0.881514751));   /* c5 */
1770     tmp3 = MULTIPLY(tmp10 + tmp13, FIX(0.710284161)) +  /* c7 */
1771            MULTIPLY(tmp14 + tmp15, FIX(0.256335874));   /* c11 */
1772     tmp0 = tmp1 + tmp2 + tmp3 -
1773            MULTIPLY(tmp10, FIX(1.530003162)) +          /* c3+c5+c7-c1 */
1774            MULTIPLY(tmp14, FIX(0.241438564));           /* c9-c11 */
1775     tmp4 = MULTIPLY(tmp14 - tmp15, FIX(0.710284161)) -  /* c7 */
1776            MULTIPLY(tmp11 + tmp12, FIX(0.256335874));   /* c11 */
1777     tmp5 = MULTIPLY(tmp11 + tmp13, - FIX(0.881514751)); /* -c5 */
1778     tmp1 += tmp4 + tmp5 +
1779             MULTIPLY(tmp11, FIX(0.634110155)) -         /* c5+c9+c11-c3 */
1780             MULTIPLY(tmp14, FIX(1.773594819));          /* c1+c7 */
1781     tmp6 = MULTIPLY(tmp12 + tmp13, - FIX(0.497774438)); /* -c9 */
1782     tmp2 += tmp4 + tmp6 -
1783             MULTIPLY(tmp12, FIX(1.190715098)) +         /* c1+c5-c9-c11 */
1784             MULTIPLY(tmp15, FIX(1.711799069));          /* c3+c7 */
1785     tmp3 += tmp5 + tmp6 +
1786             MULTIPLY(tmp13, FIX(1.670519935)) -         /* c3+c5+c9-c7 */
1787             MULTIPLY(tmp15, FIX(1.319646532));          /* c1+c11 */
1788
1789     dataptr[DCTSIZE*1] = (DCTELEM) DESCALE(tmp0, CONST_BITS+1);
1790     dataptr[DCTSIZE*3] = (DCTELEM) DESCALE(tmp1, CONST_BITS+1);
1791     dataptr[DCTSIZE*5] = (DCTELEM) DESCALE(tmp2, CONST_BITS+1);
1792     dataptr[DCTSIZE*7] = (DCTELEM) DESCALE(tmp3, CONST_BITS+1);
1793
1794     dataptr++;                  /* advance pointer to next column */
1795     wsptr++;                    /* advance pointer to next column */
1796   }
1797 }
1798
1799
1800 /*
1801  * Perform the forward DCT on a 14x14 sample block.
1802  */
1803
1804 GLOBAL(void)
1805 jpeg_fdct_14x14 (DCTELEM * data, JSAMPARRAY sample_data, JDIMENSION start_col)
1806 {
1807   INT32 tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6;
1808   INT32 tmp10, tmp11, tmp12, tmp13, tmp14, tmp15, tmp16;
1809   DCTELEM workspace[8*6];
1810   DCTELEM *dataptr;
1811   DCTELEM *wsptr;
1812   JSAMPROW elemptr;
1813   int ctr;
1814   SHIFT_TEMPS
1815
1816   /* Pass 1: process rows.
1817    * Note results are scaled up by sqrt(8) compared to a true DCT.
1818    * cK represents sqrt(2) * cos(K*pi/28).
1819    */
1820
1821   dataptr = data;
1822   ctr = 0;
1823   for (;;) {
1824     elemptr = sample_data[ctr] + start_col;
1825
1826     /* Even part */
1827
1828     tmp0 = GETJSAMPLE(elemptr[0]) + GETJSAMPLE(elemptr[13]);
1829     tmp1 = GETJSAMPLE(elemptr[1]) + GETJSAMPLE(elemptr[12]);
1830     tmp2 = GETJSAMPLE(elemptr[2]) + GETJSAMPLE(elemptr[11]);
1831     tmp13 = GETJSAMPLE(elemptr[3]) + GETJSAMPLE(elemptr[10]);
1832     tmp4 = GETJSAMPLE(elemptr[4]) + GETJSAMPLE(elemptr[9]);
1833     tmp5 = GETJSAMPLE(elemptr[5]) + GETJSAMPLE(elemptr[8]);
1834     tmp6 = GETJSAMPLE(elemptr[6]) + GETJSAMPLE(elemptr[7]);
1835
1836     tmp10 = tmp0 + tmp6;
1837     tmp14 = tmp0 - tmp6;
1838     tmp11 = tmp1 + tmp5;
1839     tmp15 = tmp1 - tmp5;
1840     tmp12 = tmp2 + tmp4;
1841     tmp16 = tmp2 - tmp4;
1842
1843     tmp0 = GETJSAMPLE(elemptr[0]) - GETJSAMPLE(elemptr[13]);
1844     tmp1 = GETJSAMPLE(elemptr[1]) - GETJSAMPLE(elemptr[12]);
1845     tmp2 = GETJSAMPLE(elemptr[2]) - GETJSAMPLE(elemptr[11]);
1846     tmp3 = GETJSAMPLE(elemptr[3]) - GETJSAMPLE(elemptr[10]);
1847     tmp4 = GETJSAMPLE(elemptr[4]) - GETJSAMPLE(elemptr[9]);
1848     tmp5 = GETJSAMPLE(elemptr[5]) - GETJSAMPLE(elemptr[8]);
1849     tmp6 = GETJSAMPLE(elemptr[6]) - GETJSAMPLE(elemptr[7]);
1850
1851     /* Apply unsigned->signed conversion */
1852     dataptr[0] = (DCTELEM)
1853       (tmp10 + tmp11 + tmp12 + tmp13 - 14 * CENTERJSAMPLE);
1854     tmp13 += tmp13;
1855     dataptr[4] = (DCTELEM)
1856       DESCALE(MULTIPLY(tmp10 - tmp13, FIX(1.274162392)) + /* c4 */
1857               MULTIPLY(tmp11 - tmp13, FIX(0.314692123)) - /* c12 */
1858               MULTIPLY(tmp12 - tmp13, FIX(0.881747734)),  /* c8 */
1859               CONST_BITS);
1860
1861     tmp10 = MULTIPLY(tmp14 + tmp15, FIX(1.105676686));    /* c6 */
1862
1863     dataptr[2] = (DCTELEM)
1864       DESCALE(tmp10 + MULTIPLY(tmp14, FIX(0.273079590))   /* c2-c6 */
1865               + MULTIPLY(tmp16, FIX(0.613604268)),        /* c10 */
1866               CONST_BITS);
1867     dataptr[6] = (DCTELEM)
1868       DESCALE(tmp10 - MULTIPLY(tmp15, FIX(1.719280954))   /* c6+c10 */
1869               - MULTIPLY(tmp16, FIX(1.378756276)),        /* c2 */
1870               CONST_BITS);
1871
1872     /* Odd part */
1873
1874     tmp10 = tmp1 + tmp2;
1875     tmp11 = tmp5 - tmp4;
1876     dataptr[7] = (DCTELEM) (tmp0 - tmp10 + tmp3 - tmp11 - tmp6);
1877     tmp3 <<= CONST_BITS;
1878     tmp10 = MULTIPLY(tmp10, - FIX(0.158341681));          /* -c13 */
1879     tmp11 = MULTIPLY(tmp11, FIX(1.405321284));            /* c1 */
1880     tmp10 += tmp11 - tmp3;
1881     tmp11 = MULTIPLY(tmp0 + tmp2, FIX(1.197448846)) +     /* c5 */
1882             MULTIPLY(tmp4 + tmp6, FIX(0.752406978));      /* c9 */
1883     dataptr[5] = (DCTELEM)
1884       DESCALE(tmp10 + tmp11 - MULTIPLY(tmp2, FIX(2.373959773)) /* c3+c5-c13 */
1885               + MULTIPLY(tmp4, FIX(1.119999435)),         /* c1+c11-c9 */
1886               CONST_BITS);
1887     tmp12 = MULTIPLY(tmp0 + tmp1, FIX(1.334852607)) +     /* c3 */
1888             MULTIPLY(tmp5 - tmp6, FIX(0.467085129));      /* c11 */
1889     dataptr[3] = (DCTELEM)
1890       DESCALE(tmp10 + tmp12 - MULTIPLY(tmp1, FIX(0.424103948)) /* c3-c9-c13 */
1891               - MULTIPLY(tmp5, FIX(3.069855259)),         /* c1+c5+c11 */
1892               CONST_BITS);
1893     dataptr[1] = (DCTELEM)
1894       DESCALE(tmp11 + tmp12 + tmp3 + tmp6 -
1895               MULTIPLY(tmp0 + tmp6, FIX(1.126980169)),    /* c3+c5-c1 */
1896               CONST_BITS);
1897
1898     ctr++;
1899
1900     if (ctr != DCTSIZE) {
1901       if (ctr == 14)
1902         break;                  /* Done. */
1903       dataptr += DCTSIZE;       /* advance pointer to next row */
1904     } else
1905       dataptr = workspace;      /* switch pointer to extended workspace */
1906   }
1907
1908   /* Pass 2: process columns.
1909    * We leave the results scaled up by an overall factor of 8.
1910    * We must also scale the output by (8/14)**2 = 16/49, which we partially
1911    * fold into the constant multipliers and final shifting:
1912    * cK now represents sqrt(2) * cos(K*pi/28) * 32/49.
1913    */
1914
1915   dataptr = data;
1916   wsptr = workspace;
1917   for (ctr = DCTSIZE-1; ctr >= 0; ctr--) {
1918     /* Even part */
1919
1920     tmp0 = dataptr[DCTSIZE*0] + wsptr[DCTSIZE*5];
1921     tmp1 = dataptr[DCTSIZE*1] + wsptr[DCTSIZE*4];
1922     tmp2 = dataptr[DCTSIZE*2] + wsptr[DCTSIZE*3];
1923     tmp13 = dataptr[DCTSIZE*3] + wsptr[DCTSIZE*2];
1924     tmp4 = dataptr[DCTSIZE*4] + wsptr[DCTSIZE*1];
1925     tmp5 = dataptr[DCTSIZE*5] + wsptr[DCTSIZE*0];
1926     tmp6 = dataptr[DCTSIZE*6] + dataptr[DCTSIZE*7];
1927
1928     tmp10 = tmp0 + tmp6;
1929     tmp14 = tmp0 - tmp6;
1930     tmp11 = tmp1 + tmp5;
1931     tmp15 = tmp1 - tmp5;
1932     tmp12 = tmp2 + tmp4;
1933     tmp16 = tmp2 - tmp4;
1934
1935     tmp0 = dataptr[DCTSIZE*0] - wsptr[DCTSIZE*5];
1936     tmp1 = dataptr[DCTSIZE*1] - wsptr[DCTSIZE*4];
1937     tmp2 = dataptr[DCTSIZE*2] - wsptr[DCTSIZE*3];
1938     tmp3 = dataptr[DCTSIZE*3] - wsptr[DCTSIZE*2];
1939     tmp4 = dataptr[DCTSIZE*4] - wsptr[DCTSIZE*1];
1940     tmp5 = dataptr[DCTSIZE*5] - wsptr[DCTSIZE*0];
1941     tmp6 = dataptr[DCTSIZE*6] - dataptr[DCTSIZE*7];
1942
1943     dataptr[DCTSIZE*0] = (DCTELEM)
1944       DESCALE(MULTIPLY(tmp10 + tmp11 + tmp12 + tmp13,
1945                        FIX(0.653061224)),                 /* 32/49 */
1946               CONST_BITS+1);
1947     tmp13 += tmp13;
1948     dataptr[DCTSIZE*4] = (DCTELEM)
1949       DESCALE(MULTIPLY(tmp10 - tmp13, FIX(0.832106052)) + /* c4 */
1950               MULTIPLY(tmp11 - tmp13, FIX(0.205513223)) - /* c12 */
1951               MULTIPLY(tmp12 - tmp13, FIX(0.575835255)),  /* c8 */
1952               CONST_BITS+1);
1953
1954     tmp10 = MULTIPLY(tmp14 + tmp15, FIX(0.722074570));    /* c6 */
1955
1956     dataptr[DCTSIZE*2] = (DCTELEM)
1957       DESCALE(tmp10 + MULTIPLY(tmp14, FIX(0.178337691))   /* c2-c6 */
1958               + MULTIPLY(tmp16, FIX(0.400721155)),        /* c10 */
1959               CONST_BITS+1);
1960     dataptr[DCTSIZE*6] = (DCTELEM)
1961       DESCALE(tmp10 - MULTIPLY(tmp15, FIX(1.122795725))   /* c6+c10 */
1962               - MULTIPLY(tmp16, FIX(0.900412262)),        /* c2 */
1963               CONST_BITS+1);
1964
1965     /* Odd part */
1966
1967     tmp10 = tmp1 + tmp2;
1968     tmp11 = tmp5 - tmp4;
1969     dataptr[DCTSIZE*7] = (DCTELEM)
1970       DESCALE(MULTIPLY(tmp0 - tmp10 + tmp3 - tmp11 - tmp6,
1971                        FIX(0.653061224)),                 /* 32/49 */
1972               CONST_BITS+1);
1973     tmp3  = MULTIPLY(tmp3 , FIX(0.653061224));            /* 32/49 */
1974     tmp10 = MULTIPLY(tmp10, - FIX(0.103406812));          /* -c13 */
1975     tmp11 = MULTIPLY(tmp11, FIX(0.917760839));            /* c1 */
1976     tmp10 += tmp11 - tmp3;
1977     tmp11 = MULTIPLY(tmp0 + tmp2, FIX(0.782007410)) +     /* c5 */
1978             MULTIPLY(tmp4 + tmp6, FIX(0.491367823));      /* c9 */
1979     dataptr[DCTSIZE*5] = (DCTELEM)
1980       DESCALE(tmp10 + tmp11 - MULTIPLY(tmp2, FIX(1.550341076)) /* c3+c5-c13 */
1981               + MULTIPLY(tmp4, FIX(0.731428202)),         /* c1+c11-c9 */
1982               CONST_BITS+1);
1983     tmp12 = MULTIPLY(tmp0 + tmp1, FIX(0.871740478)) +     /* c3 */
1984             MULTIPLY(tmp5 - tmp6, FIX(0.305035186));      /* c11 */
1985     dataptr[DCTSIZE*3] = (DCTELEM)
1986       DESCALE(tmp10 + tmp12 - MULTIPLY(tmp1, FIX(0.276965844)) /* c3-c9-c13 */
1987               - MULTIPLY(tmp5, FIX(2.004803435)),         /* c1+c5+c11 */
1988               CONST_BITS+1);
1989     dataptr[DCTSIZE*1] = (DCTELEM)
1990       DESCALE(tmp11 + tmp12 + tmp3
1991               - MULTIPLY(tmp0, FIX(0.735987049))          /* c3+c5-c1 */
1992               - MULTIPLY(tmp6, FIX(0.082925825)),         /* c9-c11-c13 */
1993               CONST_BITS+1);
1994
1995     dataptr++;                  /* advance pointer to next column */
1996     wsptr++;                    /* advance pointer to next column */
1997   }
1998 }
1999
2000
2001 /*
2002  * Perform the forward DCT on a 15x15 sample block.
2003  */
2004
2005 GLOBAL(void)
2006 jpeg_fdct_15x15 (DCTELEM * data, JSAMPARRAY sample_data, JDIMENSION start_col)
2007 {
2008   INT32 tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7;
2009   INT32 tmp10, tmp11, tmp12, tmp13, tmp14, tmp15, tmp16;
2010   INT32 z1, z2, z3;
2011   DCTELEM workspace[8*7];
2012   DCTELEM *dataptr;
2013   DCTELEM *wsptr;
2014   JSAMPROW elemptr;
2015   int ctr;
2016   SHIFT_TEMPS
2017
2018   /* Pass 1: process rows.
2019    * Note results are scaled up by sqrt(8) compared to a true DCT.
2020    * cK represents sqrt(2) * cos(K*pi/30).
2021    */
2022
2023   dataptr = data;
2024   ctr = 0;
2025   for (;;) {
2026     elemptr = sample_data[ctr] + start_col;
2027
2028     /* Even part */
2029
2030     tmp0 = GETJSAMPLE(elemptr[0]) + GETJSAMPLE(elemptr[14]);
2031     tmp1 = GETJSAMPLE(elemptr[1]) + GETJSAMPLE(elemptr[13]);
2032     tmp2 = GETJSAMPLE(elemptr[2]) + GETJSAMPLE(elemptr[12]);
2033     tmp3 = GETJSAMPLE(elemptr[3]) + GETJSAMPLE(elemptr[11]);
2034     tmp4 = GETJSAMPLE(elemptr[4]) + GETJSAMPLE(elemptr[10]);
2035     tmp5 = GETJSAMPLE(elemptr[5]) + GETJSAMPLE(elemptr[9]);
2036     tmp6 = GETJSAMPLE(elemptr[6]) + GETJSAMPLE(elemptr[8]);
2037     tmp7 = GETJSAMPLE(elemptr[7]);
2038
2039     tmp10 = GETJSAMPLE(elemptr[0]) - GETJSAMPLE(elemptr[14]);
2040     tmp11 = GETJSAMPLE(elemptr[1]) - GETJSAMPLE(elemptr[13]);
2041     tmp12 = GETJSAMPLE(elemptr[2]) - GETJSAMPLE(elemptr[12]);
2042     tmp13 = GETJSAMPLE(elemptr[3]) - GETJSAMPLE(elemptr[11]);
2043     tmp14 = GETJSAMPLE(elemptr[4]) - GETJSAMPLE(elemptr[10]);
2044     tmp15 = GETJSAMPLE(elemptr[5]) - GETJSAMPLE(elemptr[9]);
2045     tmp16 = GETJSAMPLE(elemptr[6]) - GETJSAMPLE(elemptr[8]);
2046
2047     z1 = tmp0 + tmp4 + tmp5;
2048     z2 = tmp1 + tmp3 + tmp6;
2049     z3 = tmp2 + tmp7;
2050     /* Apply unsigned->signed conversion */
2051     dataptr[0] = (DCTELEM) (z1 + z2 + z3 - 15 * CENTERJSAMPLE);
2052     z3 += z3;
2053     dataptr[6] = (DCTELEM)
2054       DESCALE(MULTIPLY(z1 - z3, FIX(1.144122806)) - /* c6 */
2055               MULTIPLY(z2 - z3, FIX(0.437016024)),  /* c12 */
2056               CONST_BITS);
2057     tmp2 += ((tmp1 + tmp4) >> 1) - tmp7 - tmp7;
2058     z1 = MULTIPLY(tmp3 - tmp2, FIX(1.531135173)) -  /* c2+c14 */
2059          MULTIPLY(tmp6 - tmp2, FIX(2.238241955));   /* c4+c8 */
2060     z2 = MULTIPLY(tmp5 - tmp2, FIX(0.798468008)) -  /* c8-c14 */
2061          MULTIPLY(tmp0 - tmp2, FIX(0.091361227));   /* c2-c4 */
2062     z3 = MULTIPLY(tmp0 - tmp3, FIX(1.383309603)) +  /* c2 */
2063          MULTIPLY(tmp6 - tmp5, FIX(0.946293579)) +  /* c8 */
2064          MULTIPLY(tmp1 - tmp4, FIX(0.790569415));   /* (c6+c12)/2 */
2065
2066     dataptr[2] = (DCTELEM) DESCALE(z1 + z3, CONST_BITS);
2067     dataptr[4] = (DCTELEM) DESCALE(z2 + z3, CONST_BITS);
2068
2069     /* Odd part */
2070
2071     tmp2 = MULTIPLY(tmp10 - tmp12 - tmp13 + tmp15 + tmp16,
2072                     FIX(1.224744871));                         /* c5 */
2073     tmp1 = MULTIPLY(tmp10 - tmp14 - tmp15, FIX(1.344997024)) + /* c3 */
2074            MULTIPLY(tmp11 - tmp13 - tmp16, FIX(0.831253876));  /* c9 */
2075     tmp12 = MULTIPLY(tmp12, FIX(1.224744871));                 /* c5 */
2076     tmp4 = MULTIPLY(tmp10 - tmp16, FIX(1.406466353)) +         /* c1 */
2077            MULTIPLY(tmp11 + tmp14, FIX(1.344997024)) +         /* c3 */
2078            MULTIPLY(tmp13 + tmp15, FIX(0.575212477));          /* c11 */
2079     tmp0 = MULTIPLY(tmp13, FIX(0.475753014)) -                 /* c7-c11 */
2080            MULTIPLY(tmp14, FIX(0.513743148)) +                 /* c3-c9 */
2081            MULTIPLY(tmp16, FIX(1.700497885)) + tmp4 + tmp12;   /* c1+c13 */
2082     tmp3 = MULTIPLY(tmp10, - FIX(0.355500862)) -               /* -(c1-c7) */
2083            MULTIPLY(tmp11, FIX(2.176250899)) -                 /* c3+c9 */
2084            MULTIPLY(tmp15, FIX(0.869244010)) + tmp4 - tmp12;   /* c11+c13 */
2085
2086     dataptr[1] = (DCTELEM) DESCALE(tmp0, CONST_BITS);
2087     dataptr[3] = (DCTELEM) DESCALE(tmp1, CONST_BITS);
2088     dataptr[5] = (DCTELEM) DESCALE(tmp2, CONST_BITS);
2089     dataptr[7] = (DCTELEM) DESCALE(tmp3, CONST_BITS);
2090
2091     ctr++;
2092
2093     if (ctr != DCTSIZE) {
2094       if (ctr == 15)
2095         break;                  /* Done. */
2096       dataptr += DCTSIZE;       /* advance pointer to next row */
2097     } else
2098       dataptr = workspace;      /* switch pointer to extended workspace */
2099   }
2100
2101   /* Pass 2: process columns.
2102    * We leave the results scaled up by an overall factor of 8.
2103    * We must also scale the output by (8/15)**2 = 64/225, which we partially
2104    * fold into the constant multipliers and final shifting:
2105    * cK now represents sqrt(2) * cos(K*pi/30) * 256/225.
2106    */
2107
2108   dataptr = data;
2109   wsptr = workspace;
2110   for (ctr = DCTSIZE-1; ctr >= 0; ctr--) {
2111     /* Even part */
2112
2113     tmp0 = dataptr[DCTSIZE*0] + wsptr[DCTSIZE*6];
2114     tmp1 = dataptr[DCTSIZE*1] + wsptr[DCTSIZE*5];
2115     tmp2 = dataptr[DCTSIZE*2] + wsptr[DCTSIZE*4];
2116     tmp3 = dataptr[DCTSIZE*3] + wsptr[DCTSIZE*3];
2117     tmp4 = dataptr[DCTSIZE*4] + wsptr[DCTSIZE*2];
2118     tmp5 = dataptr[DCTSIZE*5] + wsptr[DCTSIZE*1];
2119     tmp6 = dataptr[DCTSIZE*6] + wsptr[DCTSIZE*0];
2120     tmp7 = dataptr[DCTSIZE*7];
2121
2122     tmp10 = dataptr[DCTSIZE*0] - wsptr[DCTSIZE*6];
2123     tmp11 = dataptr[DCTSIZE*1] - wsptr[DCTSIZE*5];
2124     tmp12 = dataptr[DCTSIZE*2] - wsptr[DCTSIZE*4];
2125     tmp13 = dataptr[DCTSIZE*3] - wsptr[DCTSIZE*3];
2126     tmp14 = dataptr[DCTSIZE*4] - wsptr[DCTSIZE*2];
2127     tmp15 = dataptr[DCTSIZE*5] - wsptr[DCTSIZE*1];
2128     tmp16 = dataptr[DCTSIZE*6] - wsptr[DCTSIZE*0];
2129
2130     z1 = tmp0 + tmp4 + tmp5;
2131     z2 = tmp1 + tmp3 + tmp6;
2132     z3 = tmp2 + tmp7;
2133     dataptr[DCTSIZE*0] = (DCTELEM)
2134       DESCALE(MULTIPLY(z1 + z2 + z3, FIX(1.137777778)), /* 256/225 */
2135               CONST_BITS+2);
2136     z3 += z3;
2137     dataptr[DCTSIZE*6] = (DCTELEM)
2138       DESCALE(MULTIPLY(z1 - z3, FIX(1.301757503)) - /* c6 */
2139               MULTIPLY(z2 - z3, FIX(0.497227121)),  /* c12 */
2140               CONST_BITS+2);
2141     tmp2 += ((tmp1 + tmp4) >> 1) - tmp7 - tmp7;
2142     z1 = MULTIPLY(tmp3 - tmp2, FIX(1.742091575)) -  /* c2+c14 */
2143          MULTIPLY(tmp6 - tmp2, FIX(2.546621957));   /* c4+c8 */
2144     z2 = MULTIPLY(tmp5 - tmp2, FIX(0.908479156)) -  /* c8-c14 */
2145          MULTIPLY(tmp0 - tmp2, FIX(0.103948774));   /* c2-c4 */
2146     z3 = MULTIPLY(tmp0 - tmp3, FIX(1.573898926)) +  /* c2 */
2147          MULTIPLY(tmp6 - tmp5, FIX(1.076671805)) +  /* c8 */
2148          MULTIPLY(tmp1 - tmp4, FIX(0.899492312));   /* (c6+c12)/2 */
2149
2150     dataptr[DCTSIZE*2] = (DCTELEM) DESCALE(z1 + z3, CONST_BITS+2);
2151     dataptr[DCTSIZE*4] = (DCTELEM) DESCALE(z2 + z3, CONST_BITS+2);
2152
2153     /* Odd part */
2154
2155     tmp2 = MULTIPLY(tmp10 - tmp12 - tmp13 + tmp15 + tmp16,
2156                     FIX(1.393487498));                         /* c5 */
2157     tmp1 = MULTIPLY(tmp10 - tmp14 - tmp15, FIX(1.530307725)) + /* c3 */
2158            MULTIPLY(tmp11 - tmp13 - tmp16, FIX(0.945782187));  /* c9 */
2159     tmp12 = MULTIPLY(tmp12, FIX(1.393487498));                 /* c5 */
2160     tmp4 = MULTIPLY(tmp10 - tmp16, FIX(1.600246161)) +         /* c1 */
2161            MULTIPLY(tmp11 + tmp14, FIX(1.530307725)) +         /* c3 */
2162            MULTIPLY(tmp13 + tmp15, FIX(0.654463974));          /* c11 */
2163     tmp0 = MULTIPLY(tmp13, FIX(0.541301207)) -                 /* c7-c11 */
2164            MULTIPLY(tmp14, FIX(0.584525538)) +                 /* c3-c9 */
2165            MULTIPLY(tmp16, FIX(1.934788705)) + tmp4 + tmp12;   /* c1+c13 */
2166     tmp3 = MULTIPLY(tmp10, - FIX(0.404480980)) -               /* -(c1-c7) */
2167            MULTIPLY(tmp11, FIX(2.476089912)) -                 /* c3+c9 */
2168            MULTIPLY(tmp15, FIX(0.989006518)) + tmp4 - tmp12;   /* c11+c13 */
2169
2170     dataptr[DCTSIZE*1] = (DCTELEM) DESCALE(tmp0, CONST_BITS+2);
2171     dataptr[DCTSIZE*3] = (DCTELEM) DESCALE(tmp1, CONST_BITS+2);
2172     dataptr[DCTSIZE*5] = (DCTELEM) DESCALE(tmp2, CONST_BITS+2);
2173     dataptr[DCTSIZE*7] = (DCTELEM) DESCALE(tmp3, CONST_BITS+2);
2174
2175     dataptr++;                  /* advance pointer to next column */
2176     wsptr++;                    /* advance pointer to next column */
2177   }
2178 }
2179
2180
2181 /*
2182  * Perform the forward DCT on a 16x16 sample block.
2183  */
2184
2185 GLOBAL(void)
2186 jpeg_fdct_16x16 (DCTELEM * data, JSAMPARRAY sample_data, JDIMENSION start_col)
2187 {
2188   INT32 tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7;
2189   INT32 tmp10, tmp11, tmp12, tmp13, tmp14, tmp15, tmp16, tmp17;
2190   DCTELEM workspace[DCTSIZE2];
2191   DCTELEM *dataptr;
2192   DCTELEM *wsptr;
2193   JSAMPROW elemptr;
2194   int ctr;
2195   SHIFT_TEMPS
2196
2197   /* Pass 1: process rows.
2198    * Note results are scaled up by sqrt(8) compared to a true DCT;
2199    * furthermore, we scale the results by 2**PASS1_BITS.
2200    * cK represents sqrt(2) * cos(K*pi/32).
2201    */
2202
2203   dataptr = data;
2204   ctr = 0;
2205   for (;;) {
2206     elemptr = sample_data[ctr] + start_col;
2207
2208     /* Even part */
2209
2210     tmp0 = GETJSAMPLE(elemptr[0]) + GETJSAMPLE(elemptr[15]);
2211     tmp1 = GETJSAMPLE(elemptr[1]) + GETJSAMPLE(elemptr[14]);
2212     tmp2 = GETJSAMPLE(elemptr[2]) + GETJSAMPLE(elemptr[13]);
2213     tmp3 = GETJSAMPLE(elemptr[3]) + GETJSAMPLE(elemptr[12]);
2214     tmp4 = GETJSAMPLE(elemptr[4]) + GETJSAMPLE(elemptr[11]);
2215     tmp5 = GETJSAMPLE(elemptr[5]) + GETJSAMPLE(elemptr[10]);
2216     tmp6 = GETJSAMPLE(elemptr[6]) + GETJSAMPLE(elemptr[9]);
2217     tmp7 = GETJSAMPLE(elemptr[7]) + GETJSAMPLE(elemptr[8]);
2218
2219     tmp10 = tmp0 + tmp7;
2220     tmp14 = tmp0 - tmp7;
2221     tmp11 = tmp1 + tmp6;
2222     tmp15 = tmp1 - tmp6;
2223     tmp12 = tmp2 + tmp5;
2224     tmp16 = tmp2 - tmp5;
2225     tmp13 = tmp3 + tmp4;
2226     tmp17 = tmp3 - tmp4;
2227
2228     tmp0 = GETJSAMPLE(elemptr[0]) - GETJSAMPLE(elemptr[15]);
2229     tmp1 = GETJSAMPLE(elemptr[1]) - GETJSAMPLE(elemptr[14]);
2230     tmp2 = GETJSAMPLE(elemptr[2]) - GETJSAMPLE(elemptr[13]);
2231     tmp3 = GETJSAMPLE(elemptr[3]) - GETJSAMPLE(elemptr[12]);
2232     tmp4 = GETJSAMPLE(elemptr[4]) - GETJSAMPLE(elemptr[11]);
2233     tmp5 = GETJSAMPLE(elemptr[5]) - GETJSAMPLE(elemptr[10]);
2234     tmp6 = GETJSAMPLE(elemptr[6]) - GETJSAMPLE(elemptr[9]);
2235     tmp7 = GETJSAMPLE(elemptr[7]) - GETJSAMPLE(elemptr[8]);
2236
2237     /* Apply unsigned->signed conversion */
2238     dataptr[0] = (DCTELEM)
2239       ((tmp10 + tmp11 + tmp12 + tmp13 - 16 * CENTERJSAMPLE) << PASS1_BITS);
2240     dataptr[4] = (DCTELEM)
2241       DESCALE(MULTIPLY(tmp10 - tmp13, FIX(1.306562965)) + /* c4[16] = c2[8] */
2242               MULTIPLY(tmp11 - tmp12, FIX_0_541196100),   /* c12[16] = c6[8] */
2243               CONST_BITS-PASS1_BITS);
2244
2245     tmp10 = MULTIPLY(tmp17 - tmp15, FIX(0.275899379)) +   /* c14[16] = c7[8] */
2246             MULTIPLY(tmp14 - tmp16, FIX(1.387039845));    /* c2[16] = c1[8] */
2247
2248     dataptr[2] = (DCTELEM)
2249       DESCALE(tmp10 + MULTIPLY(tmp15, FIX(1.451774982))   /* c6+c14 */
2250               + MULTIPLY(tmp16, FIX(2.172734804)),        /* c2+c10 */
2251               CONST_BITS-PASS1_BITS);
2252     dataptr[6] = (DCTELEM)
2253       DESCALE(tmp10 - MULTIPLY(tmp14, FIX(0.211164243))   /* c2-c6 */
2254               - MULTIPLY(tmp17, FIX(1.061594338)),        /* c10+c14 */
2255               CONST_BITS-PASS1_BITS);
2256
2257     /* Odd part */
2258
2259     tmp11 = MULTIPLY(tmp0 + tmp1, FIX(1.353318001)) +         /* c3 */
2260             MULTIPLY(tmp6 - tmp7, FIX(0.410524528));          /* c13 */
2261     tmp12 = MULTIPLY(tmp0 + tmp2, FIX(1.247225013)) +         /* c5 */
2262             MULTIPLY(tmp5 + tmp7, FIX(0.666655658));          /* c11 */
2263     tmp13 = MULTIPLY(tmp0 + tmp3, FIX(1.093201867)) +         /* c7 */
2264             MULTIPLY(tmp4 - tmp7, FIX(0.897167586));          /* c9 */
2265     tmp14 = MULTIPLY(tmp1 + tmp2, FIX(0.138617169)) +         /* c15 */
2266             MULTIPLY(tmp6 - tmp5, FIX(1.407403738));          /* c1 */
2267     tmp15 = MULTIPLY(tmp1 + tmp3, - FIX(0.666655658)) +       /* -c11 */
2268             MULTIPLY(tmp4 + tmp6, - FIX(1.247225013));        /* -c5 */
2269     tmp16 = MULTIPLY(tmp2 + tmp3, - FIX(1.353318001)) +       /* -c3 */
2270             MULTIPLY(tmp5 - tmp4, FIX(0.410524528));          /* c13 */
2271     tmp10 = tmp11 + tmp12 + tmp13 -
2272             MULTIPLY(tmp0, FIX(2.286341144)) +                /* c7+c5+c3-c1 */
2273             MULTIPLY(tmp7, FIX(0.779653625));                 /* c15+c13-c11+c9 */
2274     tmp11 += tmp14 + tmp15 + MULTIPLY(tmp1, FIX(0.071888074)) /* c9-c3-c15+c11 */
2275              - MULTIPLY(tmp6, FIX(1.663905119));              /* c7+c13+c1-c5 */
2276     tmp12 += tmp14 + tmp16 - MULTIPLY(tmp2, FIX(1.125726048)) /* c7+c5+c15-c3 */
2277              + MULTIPLY(tmp5, FIX(1.227391138));              /* c9-c11+c1-c13 */
2278     tmp13 += tmp15 + tmp16 + MULTIPLY(tmp3, FIX(1.065388962)) /* c15+c3+c11-c7 */
2279              + MULTIPLY(tmp4, FIX(2.167985692));              /* c1+c13+c5-c9 */
2280
2281     dataptr[1] = (DCTELEM) DESCALE(tmp10, CONST_BITS-PASS1_BITS);
2282     dataptr[3] = (DCTELEM) DESCALE(tmp11, CONST_BITS-PASS1_BITS);
2283     dataptr[5] = (DCTELEM) DESCALE(tmp12, CONST_BITS-PASS1_BITS);
2284     dataptr[7] = (DCTELEM) DESCALE(tmp13, CONST_BITS-PASS1_BITS);
2285
2286     ctr++;
2287
2288     if (ctr != DCTSIZE) {
2289       if (ctr == DCTSIZE * 2)
2290         break;                  /* Done. */
2291       dataptr += DCTSIZE;       /* advance pointer to next row */
2292     } else
2293       dataptr = workspace;      /* switch pointer to extended workspace */
2294   }
2295
2296   /* Pass 2: process columns.
2297    * We remove the PASS1_BITS scaling, but leave the results scaled up
2298    * by an overall factor of 8.
2299    * We must also scale the output by (8/16)**2 = 1/2**2.
2300    * cK represents sqrt(2) * cos(K*pi/32).
2301    */
2302
2303   dataptr = data;
2304   wsptr = workspace;
2305   for (ctr = DCTSIZE-1; ctr >= 0; ctr--) {
2306     /* Even part */
2307
2308     tmp0 = dataptr[DCTSIZE*0] + wsptr[DCTSIZE*7];
2309     tmp1 = dataptr[DCTSIZE*1] + wsptr[DCTSIZE*6];
2310     tmp2 = dataptr[DCTSIZE*2] + wsptr[DCTSIZE*5];
2311     tmp3 = dataptr[DCTSIZE*3] + wsptr[DCTSIZE*4];
2312     tmp4 = dataptr[DCTSIZE*4] + wsptr[DCTSIZE*3];
2313     tmp5 = dataptr[DCTSIZE*5] + wsptr[DCTSIZE*2];
2314     tmp6 = dataptr[DCTSIZE*6] + wsptr[DCTSIZE*1];
2315     tmp7 = dataptr[DCTSIZE*7] + wsptr[DCTSIZE*0];
2316
2317     tmp10 = tmp0 + tmp7;
2318     tmp14 = tmp0 - tmp7;
2319     tmp11 = tmp1 + tmp6;
2320     tmp15 = tmp1 - tmp6;
2321     tmp12 = tmp2 + tmp5;
2322     tmp16 = tmp2 - tmp5;
2323     tmp13 = tmp3 + tmp4;
2324     tmp17 = tmp3 - tmp4;
2325
2326     tmp0 = dataptr[DCTSIZE*0] - wsptr[DCTSIZE*7];
2327     tmp1 = dataptr[DCTSIZE*1] - wsptr[DCTSIZE*6];
2328     tmp2 = dataptr[DCTSIZE*2] - wsptr[DCTSIZE*5];
2329     tmp3 = dataptr[DCTSIZE*3] - wsptr[DCTSIZE*4];
2330     tmp4 = dataptr[DCTSIZE*4] - wsptr[DCTSIZE*3];
2331     tmp5 = dataptr[DCTSIZE*5] - wsptr[DCTSIZE*2];
2332     tmp6 = dataptr[DCTSIZE*6] - wsptr[DCTSIZE*1];
2333     tmp7 = dataptr[DCTSIZE*7] - wsptr[DCTSIZE*0];
2334
2335     dataptr[DCTSIZE*0] = (DCTELEM)
2336       DESCALE(tmp10 + tmp11 + tmp12 + tmp13, PASS1_BITS+2);
2337     dataptr[DCTSIZE*4] = (DCTELEM)
2338       DESCALE(MULTIPLY(tmp10 - tmp13, FIX(1.306562965)) + /* c4[16] = c2[8] */
2339               MULTIPLY(tmp11 - tmp12, FIX_0_541196100),   /* c12[16] = c6[8] */
2340               CONST_BITS+PASS1_BITS+2);
2341
2342     tmp10 = MULTIPLY(tmp17 - tmp15, FIX(0.275899379)) +   /* c14[16] = c7[8] */
2343             MULTIPLY(tmp14 - tmp16, FIX(1.387039845));    /* c2[16] = c1[8] */
2344
2345     dataptr[DCTSIZE*2] = (DCTELEM)
2346       DESCALE(tmp10 + MULTIPLY(tmp15, FIX(1.451774982))   /* c6+c14 */
2347               + MULTIPLY(tmp16, FIX(2.172734804)),        /* c2+10 */
2348               CONST_BITS+PASS1_BITS+2);
2349     dataptr[DCTSIZE*6] = (DCTELEM)
2350       DESCALE(tmp10 - MULTIPLY(tmp14, FIX(0.211164243))   /* c2-c6 */
2351               - MULTIPLY(tmp17, FIX(1.061594338)),        /* c10+c14 */
2352               CONST_BITS+PASS1_BITS+2);
2353
2354     /* Odd part */
2355
2356     tmp11 = MULTIPLY(tmp0 + tmp1, FIX(1.353318001)) +         /* c3 */
2357             MULTIPLY(tmp6 - tmp7, FIX(0.410524528));          /* c13 */
2358     tmp12 = MULTIPLY(tmp0 + tmp2, FIX(1.247225013)) +         /* c5 */
2359             MULTIPLY(tmp5 + tmp7, FIX(0.666655658));          /* c11 */
2360     tmp13 = MULTIPLY(tmp0 + tmp3, FIX(1.093201867)) +         /* c7 */
2361             MULTIPLY(tmp4 - tmp7, FIX(0.897167586));          /* c9 */
2362     tmp14 = MULTIPLY(tmp1 + tmp2, FIX(0.138617169)) +         /* c15 */
2363             MULTIPLY(tmp6 - tmp5, FIX(1.407403738));          /* c1 */
2364     tmp15 = MULTIPLY(tmp1 + tmp3, - FIX(0.666655658)) +       /* -c11 */
2365             MULTIPLY(tmp4 + tmp6, - FIX(1.247225013));        /* -c5 */
2366     tmp16 = MULTIPLY(tmp2 + tmp3, - FIX(1.353318001)) +       /* -c3 */
2367             MULTIPLY(tmp5 - tmp4, FIX(0.410524528));          /* c13 */
2368     tmp10 = tmp11 + tmp12 + tmp13 -
2369             MULTIPLY(tmp0, FIX(2.286341144)) +                /* c7+c5+c3-c1 */
2370             MULTIPLY(tmp7, FIX(0.779653625));                 /* c15+c13-c11+c9 */
2371     tmp11 += tmp14 + tmp15 + MULTIPLY(tmp1, FIX(0.071888074)) /* c9-c3-c15+c11 */
2372              - MULTIPLY(tmp6, FIX(1.663905119));              /* c7+c13+c1-c5 */
2373     tmp12 += tmp14 + tmp16 - MULTIPLY(tmp2, FIX(1.125726048)) /* c7+c5+c15-c3 */
2374              + MULTIPLY(tmp5, FIX(1.227391138));              /* c9-c11+c1-c13 */
2375     tmp13 += tmp15 + tmp16 + MULTIPLY(tmp3, FIX(1.065388962)) /* c15+c3+c11-c7 */
2376              + MULTIPLY(tmp4, FIX(2.167985692));              /* c1+c13+c5-c9 */
2377
2378     dataptr[DCTSIZE*1] = (DCTELEM) DESCALE(tmp10, CONST_BITS+PASS1_BITS+2);
2379     dataptr[DCTSIZE*3] = (DCTELEM) DESCALE(tmp11, CONST_BITS+PASS1_BITS+2);
2380     dataptr[DCTSIZE*5] = (DCTELEM) DESCALE(tmp12, CONST_BITS+PASS1_BITS+2);
2381     dataptr[DCTSIZE*7] = (DCTELEM) DESCALE(tmp13, CONST_BITS+PASS1_BITS+2);
2382
2383     dataptr++;                  /* advance pointer to next column */
2384     wsptr++;                    /* advance pointer to next column */
2385   }
2386 }
2387
2388
2389 /*
2390  * Perform the forward DCT on a 16x8 sample block.
2391  *
2392  * 16-point FDCT in pass 1 (rows), 8-point in pass 2 (columns).
2393  */
2394
2395 GLOBAL(void)
2396 jpeg_fdct_16x8 (DCTELEM * data, JSAMPARRAY sample_data, JDIMENSION start_col)
2397 {
2398   INT32 tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7;
2399   INT32 tmp10, tmp11, tmp12, tmp13, tmp14, tmp15, tmp16, tmp17;
2400   INT32 z1;
2401   DCTELEM *dataptr;
2402   JSAMPROW elemptr;
2403   int ctr;
2404   SHIFT_TEMPS
2405
2406   /* Pass 1: process rows.
2407    * Note results are scaled up by sqrt(8) compared to a true DCT;
2408    * furthermore, we scale the results by 2**PASS1_BITS.
2409    * 16-point FDCT kernel, cK represents sqrt(2) * cos(K*pi/32).
2410    */
2411
2412   dataptr = data;
2413   ctr = 0;
2414   for (ctr = 0; ctr < DCTSIZE; ctr++) {
2415     elemptr = sample_data[ctr] + start_col;
2416
2417     /* Even part */
2418
2419     tmp0 = GETJSAMPLE(elemptr[0]) + GETJSAMPLE(elemptr[15]);
2420     tmp1 = GETJSAMPLE(elemptr[1]) + GETJSAMPLE(elemptr[14]);
2421     tmp2 = GETJSAMPLE(elemptr[2]) + GETJSAMPLE(elemptr[13]);
2422     tmp3 = GETJSAMPLE(elemptr[3]) + GETJSAMPLE(elemptr[12]);
2423     tmp4 = GETJSAMPLE(elemptr[4]) + GETJSAMPLE(elemptr[11]);
2424     tmp5 = GETJSAMPLE(elemptr[5]) + GETJSAMPLE(elemptr[10]);
2425     tmp6 = GETJSAMPLE(elemptr[6]) + GETJSAMPLE(elemptr[9]);
2426     tmp7 = GETJSAMPLE(elemptr[7]) + GETJSAMPLE(elemptr[8]);
2427
2428     tmp10 = tmp0 + tmp7;
2429     tmp14 = tmp0 - tmp7;
2430     tmp11 = tmp1 + tmp6;
2431     tmp15 = tmp1 - tmp6;
2432     tmp12 = tmp2 + tmp5;
2433     tmp16 = tmp2 - tmp5;
2434     tmp13 = tmp3 + tmp4;
2435     tmp17 = tmp3 - tmp4;
2436
2437     tmp0 = GETJSAMPLE(elemptr[0]) - GETJSAMPLE(elemptr[15]);
2438     tmp1 = GETJSAMPLE(elemptr[1]) - GETJSAMPLE(elemptr[14]);
2439     tmp2 = GETJSAMPLE(elemptr[2]) - GETJSAMPLE(elemptr[13]);
2440     tmp3 = GETJSAMPLE(elemptr[3]) - GETJSAMPLE(elemptr[12]);
2441     tmp4 = GETJSAMPLE(elemptr[4]) - GETJSAMPLE(elemptr[11]);
2442     tmp5 = GETJSAMPLE(elemptr[5]) - GETJSAMPLE(elemptr[10]);
2443     tmp6 = GETJSAMPLE(elemptr[6]) - GETJSAMPLE(elemptr[9]);
2444     tmp7 = GETJSAMPLE(elemptr[7]) - GETJSAMPLE(elemptr[8]);
2445
2446     /* Apply unsigned->signed conversion */
2447     dataptr[0] = (DCTELEM)
2448       ((tmp10 + tmp11 + tmp12 + tmp13 - 16 * CENTERJSAMPLE) << PASS1_BITS);
2449     dataptr[4] = (DCTELEM)
2450       DESCALE(MULTIPLY(tmp10 - tmp13, FIX(1.306562965)) + /* c4[16] = c2[8] */
2451               MULTIPLY(tmp11 - tmp12, FIX_0_541196100),   /* c12[16] = c6[8] */
2452               CONST_BITS-PASS1_BITS);
2453
2454     tmp10 = MULTIPLY(tmp17 - tmp15, FIX(0.275899379)) +   /* c14[16] = c7[8] */
2455             MULTIPLY(tmp14 - tmp16, FIX(1.387039845));    /* c2[16] = c1[8] */
2456
2457     dataptr[2] = (DCTELEM)
2458       DESCALE(tmp10 + MULTIPLY(tmp15, FIX(1.451774982))   /* c6+c14 */
2459               + MULTIPLY(tmp16, FIX(2.172734804)),        /* c2+c10 */
2460               CONST_BITS-PASS1_BITS);
2461     dataptr[6] = (DCTELEM)
2462       DESCALE(tmp10 - MULTIPLY(tmp14, FIX(0.211164243))   /* c2-c6 */
2463               - MULTIPLY(tmp17, FIX(1.061594338)),        /* c10+c14 */
2464               CONST_BITS-PASS1_BITS);
2465
2466     /* Odd part */
2467
2468     tmp11 = MULTIPLY(tmp0 + tmp1, FIX(1.353318001)) +         /* c3 */
2469             MULTIPLY(tmp6 - tmp7, FIX(0.410524528));          /* c13 */
2470     tmp12 = MULTIPLY(tmp0 + tmp2, FIX(1.247225013)) +         /* c5 */
2471             MULTIPLY(tmp5 + tmp7, FIX(0.666655658));          /* c11 */
2472     tmp13 = MULTIPLY(tmp0 + tmp3, FIX(1.093201867)) +         /* c7 */
2473             MULTIPLY(tmp4 - tmp7, FIX(0.897167586));          /* c9 */
2474     tmp14 = MULTIPLY(tmp1 + tmp2, FIX(0.138617169)) +         /* c15 */
2475             MULTIPLY(tmp6 - tmp5, FIX(1.407403738));          /* c1 */
2476     tmp15 = MULTIPLY(tmp1 + tmp3, - FIX(0.666655658)) +       /* -c11 */
2477             MULTIPLY(tmp4 + tmp6, - FIX(1.247225013));        /* -c5 */
2478     tmp16 = MULTIPLY(tmp2 + tmp3, - FIX(1.353318001)) +       /* -c3 */
2479             MULTIPLY(tmp5 - tmp4, FIX(0.410524528));          /* c13 */
2480     tmp10 = tmp11 + tmp12 + tmp13 -
2481             MULTIPLY(tmp0, FIX(2.286341144)) +                /* c7+c5+c3-c1 */
2482             MULTIPLY(tmp7, FIX(0.779653625));                 /* c15+c13-c11+c9 */
2483     tmp11 += tmp14 + tmp15 + MULTIPLY(tmp1, FIX(0.071888074)) /* c9-c3-c15+c11 */
2484              - MULTIPLY(tmp6, FIX(1.663905119));              /* c7+c13+c1-c5 */
2485     tmp12 += tmp14 + tmp16 - MULTIPLY(tmp2, FIX(1.125726048)) /* c7+c5+c15-c3 */
2486              + MULTIPLY(tmp5, FIX(1.227391138));              /* c9-c11+c1-c13 */
2487     tmp13 += tmp15 + tmp16 + MULTIPLY(tmp3, FIX(1.065388962)) /* c15+c3+c11-c7 */
2488              + MULTIPLY(tmp4, FIX(2.167985692));              /* c1+c13+c5-c9 */
2489
2490     dataptr[1] = (DCTELEM) DESCALE(tmp10, CONST_BITS-PASS1_BITS);
2491     dataptr[3] = (DCTELEM) DESCALE(tmp11, CONST_BITS-PASS1_BITS);
2492     dataptr[5] = (DCTELEM) DESCALE(tmp12, CONST_BITS-PASS1_BITS);
2493     dataptr[7] = (DCTELEM) DESCALE(tmp13, CONST_BITS-PASS1_BITS);
2494
2495     dataptr += DCTSIZE;         /* advance pointer to next row */
2496   }
2497
2498   /* Pass 2: process columns.
2499    * We remove the PASS1_BITS scaling, but leave the results scaled up
2500    * by an overall factor of 8.
2501    * We must also scale the output by 8/16 = 1/2.
2502    * 8-point FDCT kernel, cK represents sqrt(2) * cos(K*pi/16).
2503    */
2504
2505   dataptr = data;
2506   for (ctr = DCTSIZE-1; ctr >= 0; ctr--) {
2507     /* Even part per LL&M figure 1 --- note that published figure is faulty;
2508      * rotator "c1" should be "c6".
2509      */
2510
2511     tmp0 = dataptr[DCTSIZE*0] + dataptr[DCTSIZE*7];
2512     tmp1 = dataptr[DCTSIZE*1] + dataptr[DCTSIZE*6];
2513     tmp2 = dataptr[DCTSIZE*2] + dataptr[DCTSIZE*5];
2514     tmp3 = dataptr[DCTSIZE*3] + dataptr[DCTSIZE*4];
2515
2516     tmp10 = tmp0 + tmp3;
2517     tmp12 = tmp0 - tmp3;
2518     tmp11 = tmp1 + tmp2;
2519     tmp13 = tmp1 - tmp2;
2520
2521     tmp0 = dataptr[DCTSIZE*0] - dataptr[DCTSIZE*7];
2522     tmp1 = dataptr[DCTSIZE*1] - dataptr[DCTSIZE*6];
2523     tmp2 = dataptr[DCTSIZE*2] - dataptr[DCTSIZE*5];
2524     tmp3 = dataptr[DCTSIZE*3] - dataptr[DCTSIZE*4];
2525
2526     dataptr[DCTSIZE*0] = (DCTELEM) DESCALE(tmp10 + tmp11, PASS1_BITS+1);
2527     dataptr[DCTSIZE*4] = (DCTELEM) DESCALE(tmp10 - tmp11, PASS1_BITS+1);
2528
2529     z1 = MULTIPLY(tmp12 + tmp13, FIX_0_541196100);   /* c6 */
2530     dataptr[DCTSIZE*2] = (DCTELEM)
2531       DESCALE(z1 + MULTIPLY(tmp12, FIX_0_765366865), /* c2-c6 */
2532               CONST_BITS+PASS1_BITS+1);
2533     dataptr[DCTSIZE*6] = (DCTELEM)
2534       DESCALE(z1 - MULTIPLY(tmp13, FIX_1_847759065), /* c2+c6 */
2535               CONST_BITS+PASS1_BITS+1);
2536
2537     /* Odd part per figure 8 --- note paper omits factor of sqrt(2).
2538      * i0..i3 in the paper are tmp0..tmp3 here.
2539      */
2540
2541     tmp12 = tmp0 + tmp2;
2542     tmp13 = tmp1 + tmp3;
2543
2544     z1 = MULTIPLY(tmp12 + tmp13, FIX_1_175875602);   /*  c3 */
2545     tmp12 = MULTIPLY(tmp12, - FIX_0_390180644);      /* -c3+c5 */
2546     tmp13 = MULTIPLY(tmp13, - FIX_1_961570560);      /* -c3-c5 */
2547     tmp12 += z1;
2548     tmp13 += z1;
2549
2550     z1 = MULTIPLY(tmp0 + tmp3, - FIX_0_899976223);   /* -c3+c7 */
2551     tmp0 = MULTIPLY(tmp0, FIX_1_501321110);          /*  c1+c3-c5-c7 */
2552     tmp3 = MULTIPLY(tmp3, FIX_0_298631336);          /* -c1+c3+c5-c7 */
2553     tmp0 += z1 + tmp12;
2554     tmp3 += z1 + tmp13;
2555
2556     z1 = MULTIPLY(tmp1 + tmp2, - FIX_2_562915447);   /* -c1-c3 */
2557     tmp1 = MULTIPLY(tmp1, FIX_3_072711026);          /*  c1+c3+c5-c7 */
2558     tmp2 = MULTIPLY(tmp2, FIX_2_053119869);          /*  c1+c3-c5+c7 */
2559     tmp1 += z1 + tmp13;
2560     tmp2 += z1 + tmp12;
2561
2562     dataptr[DCTSIZE*1] = (DCTELEM) DESCALE(tmp0, CONST_BITS+PASS1_BITS+1);
2563     dataptr[DCTSIZE*3] = (DCTELEM) DESCALE(tmp1, CONST_BITS+PASS1_BITS+1);
2564     dataptr[DCTSIZE*5] = (DCTELEM) DESCALE(tmp2, CONST_BITS+PASS1_BITS+1);
2565     dataptr[DCTSIZE*7] = (DCTELEM) DESCALE(tmp3, CONST_BITS+PASS1_BITS+1);
2566
2567     dataptr++;                  /* advance pointer to next column */
2568   }
2569 }
2570
2571
2572 /*
2573  * Perform the forward DCT on a 14x7 sample block.
2574  *
2575  * 14-point FDCT in pass 1 (rows), 7-point in pass 2 (columns).
2576  */
2577
2578 GLOBAL(void)
2579 jpeg_fdct_14x7 (DCTELEM * data, JSAMPARRAY sample_data, JDIMENSION start_col)
2580 {
2581   INT32 tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6;
2582   INT32 tmp10, tmp11, tmp12, tmp13, tmp14, tmp15, tmp16;
2583   INT32 z1, z2, z3;
2584   DCTELEM *dataptr;
2585   JSAMPROW elemptr;
2586   int ctr;
2587   SHIFT_TEMPS
2588
2589   /* Zero bottom row of output coefficient block. */
2590   MEMZERO(&data[DCTSIZE*7], SIZEOF(DCTELEM) * DCTSIZE);
2591
2592   /* Pass 1: process rows.
2593    * Note results are scaled up by sqrt(8) compared to a true DCT;
2594    * furthermore, we scale the results by 2**PASS1_BITS.
2595    * 14-point FDCT kernel, cK represents sqrt(2) * cos(K*pi/28).
2596    */
2597
2598   dataptr = data;
2599   for (ctr = 0; ctr < 7; ctr++) {
2600     elemptr = sample_data[ctr] + start_col;
2601
2602     /* Even part */
2603
2604     tmp0 = GETJSAMPLE(elemptr[0]) + GETJSAMPLE(elemptr[13]);
2605     tmp1 = GETJSAMPLE(elemptr[1]) + GETJSAMPLE(elemptr[12]);
2606     tmp2 = GETJSAMPLE(elemptr[2]) + GETJSAMPLE(elemptr[11]);
2607     tmp13 = GETJSAMPLE(elemptr[3]) + GETJSAMPLE(elemptr[10]);
2608     tmp4 = GETJSAMPLE(elemptr[4]) + GETJSAMPLE(elemptr[9]);
2609     tmp5 = GETJSAMPLE(elemptr[5]) + GETJSAMPLE(elemptr[8]);
2610     tmp6 = GETJSAMPLE(elemptr[6]) + GETJSAMPLE(elemptr[7]);
2611
2612     tmp10 = tmp0 + tmp6;
2613     tmp14 = tmp0 - tmp6;
2614     tmp11 = tmp1 + tmp5;
2615     tmp15 = tmp1 - tmp5;
2616     tmp12 = tmp2 + tmp4;
2617     tmp16 = tmp2 - tmp4;
2618
2619     tmp0 = GETJSAMPLE(elemptr[0]) - GETJSAMPLE(elemptr[13]);
2620     tmp1 = GETJSAMPLE(elemptr[1]) - GETJSAMPLE(elemptr[12]);
2621     tmp2 = GETJSAMPLE(elemptr[2]) - GETJSAMPLE(elemptr[11]);
2622     tmp3 = GETJSAMPLE(elemptr[3]) - GETJSAMPLE(elemptr[10]);
2623     tmp4 = GETJSAMPLE(elemptr[4]) - GETJSAMPLE(elemptr[9]);
2624     tmp5 = GETJSAMPLE(elemptr[5]) - GETJSAMPLE(elemptr[8]);
2625     tmp6 = GETJSAMPLE(elemptr[6]) - GETJSAMPLE(elemptr[7]);
2626
2627     /* Apply unsigned->signed conversion */
2628     dataptr[0] = (DCTELEM)
2629       ((tmp10 + tmp11 + tmp12 + tmp13 - 14 * CENTERJSAMPLE) << PASS1_BITS);
2630     tmp13 += tmp13;
2631     dataptr[4] = (DCTELEM)
2632       DESCALE(MULTIPLY(tmp10 - tmp13, FIX(1.274162392)) + /* c4 */
2633               MULTIPLY(tmp11 - tmp13, FIX(0.314692123)) - /* c12 */
2634               MULTIPLY(tmp12 - tmp13, FIX(0.881747734)),  /* c8 */
2635               CONST_BITS-PASS1_BITS);
2636
2637     tmp10 = MULTIPLY(tmp14 + tmp15, FIX(1.105676686));    /* c6 */
2638
2639     dataptr[2] = (DCTELEM)
2640       DESCALE(tmp10 + MULTIPLY(tmp14, FIX(0.273079590))   /* c2-c6 */
2641               + MULTIPLY(tmp16, FIX(0.613604268)),        /* c10 */
2642               CONST_BITS-PASS1_BITS);
2643     dataptr[6] = (DCTELEM)
2644       DESCALE(tmp10 - MULTIPLY(tmp15, FIX(1.719280954))   /* c6+c10 */
2645               - MULTIPLY(tmp16, FIX(1.378756276)),        /* c2 */
2646               CONST_BITS-PASS1_BITS);
2647
2648     /* Odd part */
2649
2650     tmp10 = tmp1 + tmp2;
2651     tmp11 = tmp5 - tmp4;
2652     dataptr[7] = (DCTELEM) ((tmp0 - tmp10 + tmp3 - tmp11 - tmp6) << PASS1_BITS);
2653     tmp3 <<= CONST_BITS;
2654     tmp10 = MULTIPLY(tmp10, - FIX(0.158341681));          /* -c13 */
2655     tmp11 = MULTIPLY(tmp11, FIX(1.405321284));            /* c1 */
2656     tmp10 += tmp11 - tmp3;
2657     tmp11 = MULTIPLY(tmp0 + tmp2, FIX(1.197448846)) +     /* c5 */
2658             MULTIPLY(tmp4 + tmp6, FIX(0.752406978));      /* c9 */
2659     dataptr[5] = (DCTELEM)
2660       DESCALE(tmp10 + tmp11 - MULTIPLY(tmp2, FIX(2.373959773)) /* c3+c5-c13 */
2661               + MULTIPLY(tmp4, FIX(1.119999435)),         /* c1+c11-c9 */
2662               CONST_BITS-PASS1_BITS);
2663     tmp12 = MULTIPLY(tmp0 + tmp1, FIX(1.334852607)) +     /* c3 */
2664             MULTIPLY(tmp5 - tmp6, FIX(0.467085129));      /* c11 */
2665     dataptr[3] = (DCTELEM)
2666       DESCALE(tmp10 + tmp12 - MULTIPLY(tmp1, FIX(0.424103948)) /* c3-c9-c13 */
2667               - MULTIPLY(tmp5, FIX(3.069855259)),         /* c1+c5+c11 */
2668               CONST_BITS-PASS1_BITS);
2669     dataptr[1] = (DCTELEM)
2670       DESCALE(tmp11 + tmp12 + tmp3 + tmp6 -
2671               MULTIPLY(tmp0 + tmp6, FIX(1.126980169)),    /* c3+c5-c1 */
2672               CONST_BITS-PASS1_BITS);
2673
2674     dataptr += DCTSIZE;         /* advance pointer to next row */
2675   }
2676
2677   /* Pass 2: process columns.
2678    * We remove the PASS1_BITS scaling, but leave the results scaled up
2679    * by an overall factor of 8.
2680    * We must also scale the output by (8/14)*(8/7) = 32/49, which we
2681    * partially fold into the constant multipliers and final shifting:
2682    * 7-point FDCT kernel, cK represents sqrt(2) * cos(K*pi/14) * 64/49.
2683    */
2684
2685   dataptr = data;
2686   for (ctr = DCTSIZE-1; ctr >= 0; ctr--) {
2687     /* Even part */
2688
2689     tmp0 = dataptr[DCTSIZE*0] + dataptr[DCTSIZE*6];
2690     tmp1 = dataptr[DCTSIZE*1] + dataptr[DCTSIZE*5];
2691     tmp2 = dataptr[DCTSIZE*2] + dataptr[DCTSIZE*4];
2692     tmp3 = dataptr[DCTSIZE*3];
2693
2694     tmp10 = dataptr[DCTSIZE*0] - dataptr[DCTSIZE*6];
2695     tmp11 = dataptr[DCTSIZE*1] - dataptr[DCTSIZE*5];
2696     tmp12 = dataptr[DCTSIZE*2] - dataptr[DCTSIZE*4];
2697
2698     z1 = tmp0 + tmp2;
2699     dataptr[DCTSIZE*0] = (DCTELEM)
2700       DESCALE(MULTIPLY(z1 + tmp1 + tmp3, FIX(1.306122449)), /* 64/49 */
2701               CONST_BITS+PASS1_BITS+1);
2702     tmp3 += tmp3;
2703     z1 -= tmp3;
2704     z1 -= tmp3;
2705     z1 = MULTIPLY(z1, FIX(0.461784020));                /* (c2+c6-c4)/2 */
2706     z2 = MULTIPLY(tmp0 - tmp2, FIX(1.202428084));       /* (c2+c4-c6)/2 */
2707     z3 = MULTIPLY(tmp1 - tmp2, FIX(0.411026446));       /* c6 */
2708     dataptr[DCTSIZE*2] = (DCTELEM) DESCALE(z1 + z2 + z3, CONST_BITS+PASS1_BITS+1);
2709     z1 -= z2;
2710     z2 = MULTIPLY(tmp0 - tmp1, FIX(1.151670509));       /* c4 */
2711     dataptr[DCTSIZE*4] = (DCTELEM)
2712       DESCALE(z2 + z3 - MULTIPLY(tmp1 - tmp3, FIX(0.923568041)), /* c2+c6-c4 */
2713               CONST_BITS+PASS1_BITS+1);
2714     dataptr[DCTSIZE*6] = (DCTELEM) DESCALE(z1 + z2, CONST_BITS+PASS1_BITS+1);
2715
2716     /* Odd part */
2717
2718     tmp1 = MULTIPLY(tmp10 + tmp11, FIX(1.221765677));   /* (c3+c1-c5)/2 */
2719     tmp2 = MULTIPLY(tmp10 - tmp11, FIX(0.222383464));   /* (c3+c5-c1)/2 */
2720     tmp0 = tmp1 - tmp2;
2721     tmp1 += tmp2;
2722     tmp2 = MULTIPLY(tmp11 + tmp12, - FIX(1.800824523)); /* -c1 */
2723     tmp1 += tmp2;
2724     tmp3 = MULTIPLY(tmp10 + tmp12, FIX(0.801442310));   /* c5 */
2725     tmp0 += tmp3;
2726     tmp2 += tmp3 + MULTIPLY(tmp12, FIX(2.443531355));   /* c3+c1-c5 */
2727
2728     dataptr[DCTSIZE*1] = (DCTELEM) DESCALE(tmp0, CONST_BITS+PASS1_BITS+1);
2729     dataptr[DCTSIZE*3] = (DCTELEM) DESCALE(tmp1, CONST_BITS+PASS1_BITS+1);
2730     dataptr[DCTSIZE*5] = (DCTELEM) DESCALE(tmp2, CONST_BITS+PASS1_BITS+1);
2731
2732     dataptr++;                  /* advance pointer to next column */
2733   }
2734 }
2735
2736
2737 /*
2738  * Perform the forward DCT on a 12x6 sample block.
2739  *
2740  * 12-point FDCT in pass 1 (rows), 6-point in pass 2 (columns).
2741  */
2742
2743 GLOBAL(void)
2744 jpeg_fdct_12x6 (DCTELEM * data, JSAMPARRAY sample_data, JDIMENSION start_col)
2745 {
2746   INT32 tmp0, tmp1, tmp2, tmp3, tmp4, tmp5;
2747   INT32 tmp10, tmp11, tmp12, tmp13, tmp14, tmp15;
2748   DCTELEM *dataptr;
2749   JSAMPROW elemptr;
2750   int ctr;
2751   SHIFT_TEMPS
2752
2753   /* Zero 2 bottom rows of output coefficient block. */
2754   MEMZERO(&data[DCTSIZE*6], SIZEOF(DCTELEM) * DCTSIZE * 2);
2755
2756   /* Pass 1: process rows.
2757    * Note results are scaled up by sqrt(8) compared to a true DCT;
2758    * furthermore, we scale the results by 2**PASS1_BITS.
2759    * 12-point FDCT kernel, cK represents sqrt(2) * cos(K*pi/24).
2760    */
2761
2762   dataptr = data;
2763   for (ctr = 0; ctr < 6; ctr++) {
2764     elemptr = sample_data[ctr] + start_col;
2765
2766     /* Even part */
2767
2768     tmp0 = GETJSAMPLE(elemptr[0]) + GETJSAMPLE(elemptr[11]);
2769     tmp1 = GETJSAMPLE(elemptr[1]) + GETJSAMPLE(elemptr[10]);
2770     tmp2 = GETJSAMPLE(elemptr[2]) + GETJSAMPLE(elemptr[9]);
2771     tmp3 = GETJSAMPLE(elemptr[3]) + GETJSAMPLE(elemptr[8]);
2772     tmp4 = GETJSAMPLE(elemptr[4]) + GETJSAMPLE(elemptr[7]);
2773     tmp5 = GETJSAMPLE(elemptr[5]) + GETJSAMPLE(elemptr[6]);
2774
2775     tmp10 = tmp0 + tmp5;
2776     tmp13 = tmp0 - tmp5;
2777     tmp11 = tmp1 + tmp4;
2778     tmp14 = tmp1 - tmp4;
2779     tmp12 = tmp2 + tmp3;
2780     tmp15 = tmp2 - tmp3;
2781
2782     tmp0 = GETJSAMPLE(elemptr[0]) - GETJSAMPLE(elemptr[11]);
2783     tmp1 = GETJSAMPLE(elemptr[1]) - GETJSAMPLE(elemptr[10]);
2784     tmp2 = GETJSAMPLE(elemptr[2]) - GETJSAMPLE(elemptr[9]);
2785     tmp3 = GETJSAMPLE(elemptr[3]) - GETJSAMPLE(elemptr[8]);
2786     tmp4 = GETJSAMPLE(elemptr[4]) - GETJSAMPLE(elemptr[7]);
2787     tmp5 = GETJSAMPLE(elemptr[5]) - GETJSAMPLE(elemptr[6]);
2788
2789     /* Apply unsigned->signed conversion */
2790     dataptr[0] = (DCTELEM)
2791       ((tmp10 + tmp11 + tmp12 - 12 * CENTERJSAMPLE) << PASS1_BITS);
2792     dataptr[6] = (DCTELEM) ((tmp13 - tmp14 - tmp15) << PASS1_BITS);
2793     dataptr[4] = (DCTELEM)
2794       DESCALE(MULTIPLY(tmp10 - tmp12, FIX(1.224744871)), /* c4 */
2795               CONST_BITS-PASS1_BITS);
2796     dataptr[2] = (DCTELEM)
2797       DESCALE(tmp14 - tmp15 + MULTIPLY(tmp13 + tmp15, FIX(1.366025404)), /* c2 */
2798               CONST_BITS-PASS1_BITS);
2799
2800     /* Odd part */
2801
2802     tmp10 = MULTIPLY(tmp1 + tmp4, FIX_0_541196100);    /* c9 */
2803     tmp14 = tmp10 + MULTIPLY(tmp1, FIX_0_765366865);   /* c3-c9 */
2804     tmp15 = tmp10 - MULTIPLY(tmp4, FIX_1_847759065);   /* c3+c9 */
2805     tmp12 = MULTIPLY(tmp0 + tmp2, FIX(1.121971054));   /* c5 */
2806     tmp13 = MULTIPLY(tmp0 + tmp3, FIX(0.860918669));   /* c7 */
2807     tmp10 = tmp12 + tmp13 + tmp14 - MULTIPLY(tmp0, FIX(0.580774953)) /* c5+c7-c1 */
2808             + MULTIPLY(tmp5, FIX(0.184591911));        /* c11 */
2809     tmp11 = MULTIPLY(tmp2 + tmp3, - FIX(0.184591911)); /* -c11 */
2810     tmp12 += tmp11 - tmp15 - MULTIPLY(tmp2, FIX(2.339493912)) /* c1+c5-c11 */
2811             + MULTIPLY(tmp5, FIX(0.860918669));        /* c7 */
2812     tmp13 += tmp11 - tmp14 + MULTIPLY(tmp3, FIX(0.725788011)) /* c1+c11-c7 */
2813             - MULTIPLY(tmp5, FIX(1.121971054));        /* c5 */
2814     tmp11 = tmp15 + MULTIPLY(tmp0 - tmp3, FIX(1.306562965)) /* c3 */
2815             - MULTIPLY(tmp2 + tmp5, FIX_0_541196100);  /* c9 */
2816
2817     dataptr[1] = (DCTELEM) DESCALE(tmp10, CONST_BITS-PASS1_BITS);
2818     dataptr[3] = (DCTELEM) DESCALE(tmp11, CONST_BITS-PASS1_BITS);
2819     dataptr[5] = (DCTELEM) DESCALE(tmp12, CONST_BITS-PASS1_BITS);
2820     dataptr[7] = (DCTELEM) DESCALE(tmp13, CONST_BITS-PASS1_BITS);
2821
2822     dataptr += DCTSIZE;         /* advance pointer to next row */
2823   }
2824
2825   /* Pass 2: process columns.
2826    * We remove the PASS1_BITS scaling, but leave the results scaled up
2827    * by an overall factor of 8.
2828    * We must also scale the output by (8/12)*(8/6) = 8/9, which we
2829    * partially fold into the constant multipliers and final shifting:
2830    * 6-point FDCT kernel, cK represents sqrt(2) * cos(K*pi/12) * 16/9.
2831    */
2832
2833   dataptr = data;
2834   for (ctr = DCTSIZE-1; ctr >= 0; ctr--) {
2835     /* Even part */
2836
2837     tmp0 = dataptr[DCTSIZE*0] + dataptr[DCTSIZE*5];
2838     tmp11 = dataptr[DCTSIZE*1] + dataptr[DCTSIZE*4];
2839     tmp2 = dataptr[DCTSIZE*2] + dataptr[DCTSIZE*3];
2840
2841     tmp10 = tmp0 + tmp2;
2842     tmp12 = tmp0 - tmp2;
2843
2844     tmp0 = dataptr[DCTSIZE*0] - dataptr[DCTSIZE*5];
2845     tmp1 = dataptr[DCTSIZE*1] - dataptr[DCTSIZE*4];
2846     tmp2 = dataptr[DCTSIZE*2] - dataptr[DCTSIZE*3];
2847
2848     dataptr[DCTSIZE*0] = (DCTELEM)
2849       DESCALE(MULTIPLY(tmp10 + tmp11, FIX(1.777777778)),         /* 16/9 */
2850               CONST_BITS+PASS1_BITS+1);
2851     dataptr[DCTSIZE*2] = (DCTELEM)
2852       DESCALE(MULTIPLY(tmp12, FIX(2.177324216)),                 /* c2 */
2853               CONST_BITS+PASS1_BITS+1);
2854     dataptr[DCTSIZE*4] = (DCTELEM)
2855       DESCALE(MULTIPLY(tmp10 - tmp11 - tmp11, FIX(1.257078722)), /* c4 */
2856               CONST_BITS+PASS1_BITS+1);
2857
2858     /* Odd part */
2859
2860     tmp10 = MULTIPLY(tmp0 + tmp2, FIX(0.650711829));             /* c5 */
2861
2862     dataptr[DCTSIZE*1] = (DCTELEM)
2863       DESCALE(tmp10 + MULTIPLY(tmp0 + tmp1, FIX(1.777777778)),   /* 16/9 */
2864               CONST_BITS+PASS1_BITS+1);
2865     dataptr[DCTSIZE*3] = (DCTELEM)
2866       DESCALE(MULTIPLY(tmp0 - tmp1 - tmp2, FIX(1.777777778)),    /* 16/9 */
2867               CONST_BITS+PASS1_BITS+1);
2868     dataptr[DCTSIZE*5] = (DCTELEM)
2869       DESCALE(tmp10 + MULTIPLY(tmp2 - tmp1, FIX(1.777777778)),   /* 16/9 */
2870               CONST_BITS+PASS1_BITS+1);
2871
2872     dataptr++;                  /* advance pointer to next column */
2873   }
2874 }
2875
2876
2877 /*
2878  * Perform the forward DCT on a 10x5 sample block.
2879  *
2880  * 10-point FDCT in pass 1 (rows), 5-point in pass 2 (columns).
2881  */
2882
2883 GLOBAL(void)
2884 jpeg_fdct_10x5 (DCTELEM * data, JSAMPARRAY sample_data, JDIMENSION start_col)
2885 {
2886   INT32 tmp0, tmp1, tmp2, tmp3, tmp4;
2887   INT32 tmp10, tmp11, tmp12, tmp13, tmp14;
2888   DCTELEM *dataptr;
2889   JSAMPROW elemptr;
2890   int ctr;
2891   SHIFT_TEMPS
2892
2893   /* Zero 3 bottom rows of output coefficient block. */
2894   MEMZERO(&data[DCTSIZE*5], SIZEOF(DCTELEM) * DCTSIZE * 3);
2895
2896   /* Pass 1: process rows.
2897    * Note results are scaled up by sqrt(8) compared to a true DCT;
2898    * furthermore, we scale the results by 2**PASS1_BITS.
2899    * 10-point FDCT kernel, cK represents sqrt(2) * cos(K*pi/20).
2900    */
2901
2902   dataptr = data;
2903   for (ctr = 0; ctr < 5; ctr++) {
2904     elemptr = sample_data[ctr] + start_col;
2905
2906     /* Even part */
2907
2908     tmp0 = GETJSAMPLE(elemptr[0]) + GETJSAMPLE(elemptr[9]);
2909     tmp1 = GETJSAMPLE(elemptr[1]) + GETJSAMPLE(elemptr[8]);
2910     tmp12 = GETJSAMPLE(elemptr[2]) + GETJSAMPLE(elemptr[7]);
2911     tmp3 = GETJSAMPLE(elemptr[3]) + GETJSAMPLE(elemptr[6]);
2912     tmp4 = GETJSAMPLE(elemptr[4]) + GETJSAMPLE(elemptr[5]);
2913
2914     tmp10 = tmp0 + tmp4;
2915     tmp13 = tmp0 - tmp4;
2916     tmp11 = tmp1 + tmp3;
2917     tmp14 = tmp1 - tmp3;
2918
2919     tmp0 = GETJSAMPLE(elemptr[0]) - GETJSAMPLE(elemptr[9]);
2920     tmp1 = GETJSAMPLE(elemptr[1]) - GETJSAMPLE(elemptr[8]);
2921     tmp2 = GETJSAMPLE(elemptr[2]) - GETJSAMPLE(elemptr[7]);
2922     tmp3 = GETJSAMPLE(elemptr[3]) - GETJSAMPLE(elemptr[6]);
2923     tmp4 = GETJSAMPLE(elemptr[4]) - GETJSAMPLE(elemptr[5]);
2924
2925     /* Apply unsigned->signed conversion */
2926     dataptr[0] = (DCTELEM)
2927       ((tmp10 + tmp11 + tmp12 - 10 * CENTERJSAMPLE) << PASS1_BITS);
2928     tmp12 += tmp12;
2929     dataptr[4] = (DCTELEM)
2930       DESCALE(MULTIPLY(tmp10 - tmp12, FIX(1.144122806)) - /* c4 */
2931               MULTIPLY(tmp11 - tmp12, FIX(0.437016024)),  /* c8 */
2932               CONST_BITS-PASS1_BITS);
2933     tmp10 = MULTIPLY(tmp13 + tmp14, FIX(0.831253876));    /* c6 */
2934     dataptr[2] = (DCTELEM)
2935       DESCALE(tmp10 + MULTIPLY(tmp13, FIX(0.513743148)),  /* c2-c6 */
2936               CONST_BITS-PASS1_BITS);
2937     dataptr[6] = (DCTELEM)
2938       DESCALE(tmp10 - MULTIPLY(tmp14, FIX(2.176250899)),  /* c2+c6 */
2939               CONST_BITS-PASS1_BITS);
2940
2941     /* Odd part */
2942
2943     tmp10 = tmp0 + tmp4;
2944     tmp11 = tmp1 - tmp3;
2945     dataptr[5] = (DCTELEM) ((tmp10 - tmp11 - tmp2) << PASS1_BITS);
2946     tmp2 <<= CONST_BITS;
2947     dataptr[1] = (DCTELEM)
2948       DESCALE(MULTIPLY(tmp0, FIX(1.396802247)) +          /* c1 */
2949               MULTIPLY(tmp1, FIX(1.260073511)) + tmp2 +   /* c3 */
2950               MULTIPLY(tmp3, FIX(0.642039522)) +          /* c7 */
2951               MULTIPLY(tmp4, FIX(0.221231742)),           /* c9 */
2952               CONST_BITS-PASS1_BITS);
2953     tmp12 = MULTIPLY(tmp0 - tmp4, FIX(0.951056516)) -     /* (c3+c7)/2 */
2954             MULTIPLY(tmp1 + tmp3, FIX(0.587785252));      /* (c1-c9)/2 */
2955     tmp13 = MULTIPLY(tmp10 + tmp11, FIX(0.309016994)) +   /* (c3-c7)/2 */
2956             (tmp11 << (CONST_BITS - 1)) - tmp2;
2957     dataptr[3] = (DCTELEM) DESCALE(tmp12 + tmp13, CONST_BITS-PASS1_BITS);
2958     dataptr[7] = (DCTELEM) DESCALE(tmp12 - tmp13, CONST_BITS-PASS1_BITS);
2959
2960     dataptr += DCTSIZE;         /* advance pointer to next row */
2961   }
2962
2963   /* Pass 2: process columns.
2964    * We remove the PASS1_BITS scaling, but leave the results scaled up
2965    * by an overall factor of 8.
2966    * We must also scale the output by (8/10)*(8/5) = 32/25, which we
2967    * fold into the constant multipliers:
2968    * 5-point FDCT kernel, cK represents sqrt(2) * cos(K*pi/10) * 32/25.
2969    */
2970
2971   dataptr = data;
2972   for (ctr = DCTSIZE-1; ctr >= 0; ctr--) {
2973     /* Even part */
2974
2975     tmp0 = dataptr[DCTSIZE*0] + dataptr[DCTSIZE*4];
2976     tmp1 = dataptr[DCTSIZE*1] + dataptr[DCTSIZE*3];
2977     tmp2 = dataptr[DCTSIZE*2];
2978
2979     tmp10 = tmp0 + tmp1;
2980     tmp11 = tmp0 - tmp1;
2981
2982     tmp0 = dataptr[DCTSIZE*0] - dataptr[DCTSIZE*4];
2983     tmp1 = dataptr[DCTSIZE*1] - dataptr[DCTSIZE*3];
2984
2985     dataptr[DCTSIZE*0] = (DCTELEM)
2986       DESCALE(MULTIPLY(tmp10 + tmp2, FIX(1.28)),        /* 32/25 */
2987               CONST_BITS+PASS1_BITS);
2988     tmp11 = MULTIPLY(tmp11, FIX(1.011928851));          /* (c2+c4)/2 */
2989     tmp10 -= tmp2 << 2;
2990     tmp10 = MULTIPLY(tmp10, FIX(0.452548340));          /* (c2-c4)/2 */
2991     dataptr[DCTSIZE*2] = (DCTELEM) DESCALE(tmp11 + tmp10, CONST_BITS+PASS1_BITS);
2992     dataptr[DCTSIZE*4] = (DCTELEM) DESCALE(tmp11 - tmp10, CONST_BITS+PASS1_BITS);
2993
2994     /* Odd part */
2995
2996     tmp10 = MULTIPLY(tmp0 + tmp1, FIX(1.064004961));    /* c3 */
2997
2998     dataptr[DCTSIZE*1] = (DCTELEM)
2999       DESCALE(tmp10 + MULTIPLY(tmp0, FIX(0.657591230)), /* c1-c3 */
3000               CONST_BITS+PASS1_BITS);
3001     dataptr[DCTSIZE*3] = (DCTELEM)
3002       DESCALE(tmp10 - MULTIPLY(tmp1, FIX(2.785601151)), /* c1+c3 */
3003               CONST_BITS+PASS1_BITS);
3004
3005     dataptr++;                  /* advance pointer to next column */
3006   }
3007 }
3008
3009
3010 /*
3011  * Perform the forward DCT on an 8x4 sample block.
3012  *
3013  * 8-point FDCT in pass 1 (rows), 4-point in pass 2 (columns).
3014  */
3015
3016 GLOBAL(void)
3017 jpeg_fdct_8x4 (DCTELEM * data, JSAMPARRAY sample_data, JDIMENSION start_col)
3018 {
3019   INT32 tmp0, tmp1, tmp2, tmp3;
3020   INT32 tmp10, tmp11, tmp12, tmp13;
3021   INT32 z1;
3022   DCTELEM *dataptr;
3023   JSAMPROW elemptr;
3024   int ctr;
3025   SHIFT_TEMPS
3026
3027   /* Zero 4 bottom rows of output coefficient block. */
3028   MEMZERO(&data[DCTSIZE*4], SIZEOF(DCTELEM) * DCTSIZE * 4);
3029
3030   /* Pass 1: process rows.
3031    * Note results are scaled up by sqrt(8) compared to a true DCT;
3032    * furthermore, we scale the results by 2**PASS1_BITS.
3033    * We must also scale the output by 8/4 = 2, which we add here.
3034    * 8-point FDCT kernel, cK represents sqrt(2) * cos(K*pi/16).
3035    */
3036
3037   dataptr = data;
3038   for (ctr = 0; ctr < 4; ctr++) {
3039     elemptr = sample_data[ctr] + start_col;
3040
3041     /* Even part per LL&M figure 1 --- note that published figure is faulty;
3042      * rotator "c1" should be "c6".
3043      */
3044
3045     tmp0 = GETJSAMPLE(elemptr[0]) + GETJSAMPLE(elemptr[7]);
3046     tmp1 = GETJSAMPLE(elemptr[1]) + GETJSAMPLE(elemptr[6]);
3047     tmp2 = GETJSAMPLE(elemptr[2]) + GETJSAMPLE(elemptr[5]);
3048     tmp3 = GETJSAMPLE(elemptr[3]) + GETJSAMPLE(elemptr[4]);
3049
3050     tmp10 = tmp0 + tmp3;
3051     tmp12 = tmp0 - tmp3;
3052     tmp11 = tmp1 + tmp2;
3053     tmp13 = tmp1 - tmp2;
3054
3055     tmp0 = GETJSAMPLE(elemptr[0]) - GETJSAMPLE(elemptr[7]);
3056     tmp1 = GETJSAMPLE(elemptr[1]) - GETJSAMPLE(elemptr[6]);
3057     tmp2 = GETJSAMPLE(elemptr[2]) - GETJSAMPLE(elemptr[5]);
3058     tmp3 = GETJSAMPLE(elemptr[3]) - GETJSAMPLE(elemptr[4]);
3059
3060     /* Apply unsigned->signed conversion */
3061     dataptr[0] = (DCTELEM)
3062       ((tmp10 + tmp11 - 8 * CENTERJSAMPLE) << (PASS1_BITS+1));
3063     dataptr[4] = (DCTELEM) ((tmp10 - tmp11) << (PASS1_BITS+1));
3064
3065     z1 = MULTIPLY(tmp12 + tmp13, FIX_0_541196100);       /* c6 */
3066     /* Add fudge factor here for final descale. */
3067     z1 += ONE << (CONST_BITS-PASS1_BITS-2);
3068
3069     dataptr[2] = (DCTELEM)
3070       RIGHT_SHIFT(z1 + MULTIPLY(tmp12, FIX_0_765366865), /* c2-c6 */
3071                   CONST_BITS-PASS1_BITS-1);
3072     dataptr[6] = (DCTELEM)
3073       RIGHT_SHIFT(z1 - MULTIPLY(tmp13, FIX_1_847759065), /* c2+c6 */
3074                   CONST_BITS-PASS1_BITS-1);
3075
3076     /* Odd part per figure 8 --- note paper omits factor of sqrt(2).
3077      * i0..i3 in the paper are tmp0..tmp3 here.
3078      */
3079
3080     tmp12 = tmp0 + tmp2;
3081     tmp13 = tmp1 + tmp3;
3082
3083     z1 = MULTIPLY(tmp12 + tmp13, FIX_1_175875602);       /*  c3 */
3084     /* Add fudge factor here for final descale. */
3085     z1 += ONE << (CONST_BITS-PASS1_BITS-2);
3086
3087     tmp12 = MULTIPLY(tmp12, - FIX_0_390180644);          /* -c3+c5 */
3088     tmp13 = MULTIPLY(tmp13, - FIX_1_961570560);          /* -c3-c5 */
3089     tmp12 += z1;
3090     tmp13 += z1;
3091
3092     z1 = MULTIPLY(tmp0 + tmp3, - FIX_0_899976223);       /* -c3+c7 */
3093     tmp0 = MULTIPLY(tmp0, FIX_1_501321110);              /*  c1+c3-c5-c7 */
3094     tmp3 = MULTIPLY(tmp3, FIX_0_298631336);              /* -c1+c3+c5-c7 */
3095     tmp0 += z1 + tmp12;
3096     tmp3 += z1 + tmp13;
3097
3098     z1 = MULTIPLY(tmp1 + tmp2, - FIX_2_562915447);       /* -c1-c3 */
3099     tmp1 = MULTIPLY(tmp1, FIX_3_072711026);              /*  c1+c3+c5-c7 */
3100     tmp2 = MULTIPLY(tmp2, FIX_2_053119869);              /*  c1+c3-c5+c7 */
3101     tmp1 += z1 + tmp13;
3102     tmp2 += z1 + tmp12;
3103
3104     dataptr[1] = (DCTELEM) RIGHT_SHIFT(tmp0, CONST_BITS-PASS1_BITS-1);
3105     dataptr[3] = (DCTELEM) RIGHT_SHIFT(tmp1, CONST_BITS-PASS1_BITS-1);
3106     dataptr[5] = (DCTELEM) RIGHT_SHIFT(tmp2, CONST_BITS-PASS1_BITS-1);
3107     dataptr[7] = (DCTELEM) RIGHT_SHIFT(tmp3, CONST_BITS-PASS1_BITS-1);
3108
3109     dataptr += DCTSIZE;         /* advance pointer to next row */
3110   }
3111
3112   /* Pass 2: process columns.
3113    * We remove the PASS1_BITS scaling, but leave the results scaled up
3114    * by an overall factor of 8.
3115    * 4-point FDCT kernel,
3116    * cK represents sqrt(2) * cos(K*pi/16) [refers to 8-point FDCT].
3117    */
3118
3119   dataptr = data;
3120   for (ctr = DCTSIZE-1; ctr >= 0; ctr--) {
3121     /* Even part */
3122
3123     /* Add fudge factor here for final descale. */
3124     tmp0 = dataptr[DCTSIZE*0] + dataptr[DCTSIZE*3] + (ONE << (PASS1_BITS-1));
3125     tmp1 = dataptr[DCTSIZE*1] + dataptr[DCTSIZE*2];
3126
3127     tmp10 = dataptr[DCTSIZE*0] - dataptr[DCTSIZE*3];
3128     tmp11 = dataptr[DCTSIZE*1] - dataptr[DCTSIZE*2];
3129
3130     dataptr[DCTSIZE*0] = (DCTELEM) RIGHT_SHIFT(tmp0 + tmp1, PASS1_BITS);
3131     dataptr[DCTSIZE*2] = (DCTELEM) RIGHT_SHIFT(tmp0 - tmp1, PASS1_BITS);
3132
3133     /* Odd part */
3134
3135     tmp0 = MULTIPLY(tmp10 + tmp11, FIX_0_541196100);       /* c6 */
3136     /* Add fudge factor here for final descale. */
3137     tmp0 += ONE << (CONST_BITS+PASS1_BITS-1);
3138
3139     dataptr[DCTSIZE*1] = (DCTELEM)
3140       RIGHT_SHIFT(tmp0 + MULTIPLY(tmp10, FIX_0_765366865), /* c2-c6 */
3141                   CONST_BITS+PASS1_BITS);
3142     dataptr[DCTSIZE*3] = (DCTELEM)
3143       RIGHT_SHIFT(tmp0 - MULTIPLY(tmp11, FIX_1_847759065), /* c2+c6 */
3144                   CONST_BITS+PASS1_BITS);
3145
3146     dataptr++;                  /* advance pointer to next column */
3147   }
3148 }
3149
3150
3151 /*
3152  * Perform the forward DCT on a 6x3 sample block.
3153  *
3154  * 6-point FDCT in pass 1 (rows), 3-point in pass 2 (columns).
3155  */
3156
3157 GLOBAL(void)
3158 jpeg_fdct_6x3 (DCTELEM * data, JSAMPARRAY sample_data, JDIMENSION start_col)
3159 {
3160   INT32 tmp0, tmp1, tmp2;
3161   INT32 tmp10, tmp11, tmp12;
3162   DCTELEM *dataptr;
3163   JSAMPROW elemptr;
3164   int ctr;
3165   SHIFT_TEMPS
3166
3167   /* Pre-zero output coefficient block. */
3168   MEMZERO(data, SIZEOF(DCTELEM) * DCTSIZE2);
3169
3170   /* Pass 1: process rows.
3171    * Note results are scaled up by sqrt(8) compared to a true DCT;
3172    * furthermore, we scale the results by 2**PASS1_BITS.
3173    * We scale the results further by 2 as part of output adaption
3174    * scaling for different DCT size.
3175    * 6-point FDCT kernel, cK represents sqrt(2) * cos(K*pi/12).
3176    */
3177
3178   dataptr = data;
3179   for (ctr = 0; ctr < 3; ctr++) {
3180     elemptr = sample_data[ctr] + start_col;
3181
3182     /* Even part */
3183
3184     tmp0 = GETJSAMPLE(elemptr[0]) + GETJSAMPLE(elemptr[5]);
3185     tmp11 = GETJSAMPLE(elemptr[1]) + GETJSAMPLE(elemptr[4]);
3186     tmp2 = GETJSAMPLE(elemptr[2]) + GETJSAMPLE(elemptr[3]);
3187
3188     tmp10 = tmp0 + tmp2;
3189     tmp12 = tmp0 - tmp2;
3190
3191     tmp0 = GETJSAMPLE(elemptr[0]) - GETJSAMPLE(elemptr[5]);
3192     tmp1 = GETJSAMPLE(elemptr[1]) - GETJSAMPLE(elemptr[4]);
3193     tmp2 = GETJSAMPLE(elemptr[2]) - GETJSAMPLE(elemptr[3]);
3194
3195     /* Apply unsigned->signed conversion */
3196     dataptr[0] = (DCTELEM)
3197       ((tmp10 + tmp11 - 6 * CENTERJSAMPLE) << (PASS1_BITS+1));
3198     dataptr[2] = (DCTELEM)
3199       DESCALE(MULTIPLY(tmp12, FIX(1.224744871)),                 /* c2 */
3200               CONST_BITS-PASS1_BITS-1);
3201     dataptr[4] = (DCTELEM)
3202       DESCALE(MULTIPLY(tmp10 - tmp11 - tmp11, FIX(0.707106781)), /* c4 */
3203               CONST_BITS-PASS1_BITS-1);
3204
3205     /* Odd part */
3206
3207     tmp10 = DESCALE(MULTIPLY(tmp0 + tmp2, FIX(0.366025404)),     /* c5 */
3208                     CONST_BITS-PASS1_BITS-1);
3209
3210     dataptr[1] = (DCTELEM) (tmp10 + ((tmp0 + tmp1) << (PASS1_BITS+1)));
3211     dataptr[3] = (DCTELEM) ((tmp0 - tmp1 - tmp2) << (PASS1_BITS+1));
3212     dataptr[5] = (DCTELEM) (tmp10 + ((tmp2 - tmp1) << (PASS1_BITS+1)));
3213
3214     dataptr += DCTSIZE;         /* advance pointer to next row */
3215   }
3216
3217   /* Pass 2: process columns.
3218    * We remove the PASS1_BITS scaling, but leave the results scaled up
3219    * by an overall factor of 8.
3220    * We must also scale the output by (8/6)*(8/3) = 32/9, which we partially
3221    * fold into the constant multipliers (other part was done in pass 1):
3222    * 3-point FDCT kernel, cK represents sqrt(2) * cos(K*pi/6) * 16/9.
3223    */
3224
3225   dataptr = data;
3226   for (ctr = 0; ctr < 6; ctr++) {
3227     /* Even part */
3228
3229     tmp0 = dataptr[DCTSIZE*0] + dataptr[DCTSIZE*2];
3230     tmp1 = dataptr[DCTSIZE*1];
3231
3232     tmp2 = dataptr[DCTSIZE*0] - dataptr[DCTSIZE*2];
3233
3234     dataptr[DCTSIZE*0] = (DCTELEM)
3235       DESCALE(MULTIPLY(tmp0 + tmp1, FIX(1.777777778)),        /* 16/9 */
3236               CONST_BITS+PASS1_BITS);
3237     dataptr[DCTSIZE*2] = (DCTELEM)
3238       DESCALE(MULTIPLY(tmp0 - tmp1 - tmp1, FIX(1.257078722)), /* c2 */
3239               CONST_BITS+PASS1_BITS);
3240
3241     /* Odd part */
3242
3243     dataptr[DCTSIZE*1] = (DCTELEM)
3244       DESCALE(MULTIPLY(tmp2, FIX(2.177324216)),               /* c1 */
3245               CONST_BITS+PASS1_BITS);
3246
3247     dataptr++;                  /* advance pointer to next column */
3248   }
3249 }
3250
3251
3252 /*
3253  * Perform the forward DCT on a 4x2 sample block.
3254  *
3255  * 4-point FDCT in pass 1 (rows), 2-point in pass 2 (columns).
3256  */
3257
3258 GLOBAL(void)
3259 jpeg_fdct_4x2 (DCTELEM * data, JSAMPARRAY sample_data, JDIMENSION start_col)
3260 {
3261   INT32 tmp0, tmp1;
3262   INT32 tmp10, tmp11;
3263   DCTELEM *dataptr;
3264   JSAMPROW elemptr;
3265   int ctr;
3266   SHIFT_TEMPS
3267
3268   /* Pre-zero output coefficient block. */
3269   MEMZERO(data, SIZEOF(DCTELEM) * DCTSIZE2);
3270
3271   /* Pass 1: process rows.
3272    * Note results are scaled up by sqrt(8) compared to a true DCT;
3273    * furthermore, we scale the results by 2**PASS1_BITS.
3274    * We must also scale the output by (8/4)*(8/2) = 2**3, which we add here.
3275    * 4-point FDCT kernel,
3276    * cK represents sqrt(2) * cos(K*pi/16) [refers to 8-point FDCT].
3277    */
3278
3279   dataptr = data;
3280   for (ctr = 0; ctr < 2; ctr++) {
3281     elemptr = sample_data[ctr] + start_col;
3282
3283     /* Even part */
3284
3285     tmp0 = GETJSAMPLE(elemptr[0]) + GETJSAMPLE(elemptr[3]);
3286     tmp1 = GETJSAMPLE(elemptr[1]) + GETJSAMPLE(elemptr[2]);
3287
3288     tmp10 = GETJSAMPLE(elemptr[0]) - GETJSAMPLE(elemptr[3]);
3289     tmp11 = GETJSAMPLE(elemptr[1]) - GETJSAMPLE(elemptr[2]);
3290
3291     /* Apply unsigned->signed conversion */
3292     dataptr[0] = (DCTELEM)
3293       ((tmp0 + tmp1 - 4 * CENTERJSAMPLE) << (PASS1_BITS+3));
3294     dataptr[2] = (DCTELEM) ((tmp0 - tmp1) << (PASS1_BITS+3));
3295
3296     /* Odd part */
3297
3298     tmp0 = MULTIPLY(tmp10 + tmp11, FIX_0_541196100);       /* c6 */
3299     /* Add fudge factor here for final descale. */
3300     tmp0 += ONE << (CONST_BITS-PASS1_BITS-4);
3301
3302     dataptr[1] = (DCTELEM)
3303       RIGHT_SHIFT(tmp0 + MULTIPLY(tmp10, FIX_0_765366865), /* c2-c6 */
3304                   CONST_BITS-PASS1_BITS-3);
3305     dataptr[3] = (DCTELEM)
3306       RIGHT_SHIFT(tmp0 - MULTIPLY(tmp11, FIX_1_847759065), /* c2+c6 */
3307                   CONST_BITS-PASS1_BITS-3);
3308
3309     dataptr += DCTSIZE;         /* advance pointer to next row */
3310   }
3311
3312   /* Pass 2: process columns.
3313    * We remove the PASS1_BITS scaling, but leave the results scaled up
3314    * by an overall factor of 8.
3315    */
3316
3317   dataptr = data;
3318   for (ctr = 0; ctr < 4; ctr++) {
3319     /* Even part */
3320
3321     /* Add fudge factor here for final descale. */
3322     tmp0 = dataptr[DCTSIZE*0] + (ONE << (PASS1_BITS-1));
3323     tmp1 = dataptr[DCTSIZE*1];
3324
3325     dataptr[DCTSIZE*0] = (DCTELEM) RIGHT_SHIFT(tmp0 + tmp1, PASS1_BITS);
3326
3327     /* Odd part */
3328
3329     dataptr[DCTSIZE*1] = (DCTELEM) RIGHT_SHIFT(tmp0 - tmp1, PASS1_BITS);
3330
3331     dataptr++;                  /* advance pointer to next column */
3332   }
3333 }
3334
3335
3336 /*
3337  * Perform the forward DCT on a 2x1 sample block.
3338  *
3339  * 2-point FDCT in pass 1 (rows), 1-point in pass 2 (columns).
3340  */
3341
3342 GLOBAL(void)
3343 jpeg_fdct_2x1 (DCTELEM * data, JSAMPARRAY sample_data, JDIMENSION start_col)
3344 {
3345   INT32 tmp0, tmp1;
3346   JSAMPROW elemptr;
3347
3348   /* Pre-zero output coefficient block. */
3349   MEMZERO(data, SIZEOF(DCTELEM) * DCTSIZE2);
3350
3351   elemptr = sample_data[0] + start_col;
3352
3353   tmp0 = GETJSAMPLE(elemptr[0]);
3354   tmp1 = GETJSAMPLE(elemptr[1]);
3355
3356   /* We leave the results scaled up by an overall factor of 8.
3357    * We must also scale the output by (8/2)*(8/1) = 2**5.
3358    */
3359
3360   /* Even part */
3361
3362   /* Apply unsigned->signed conversion */
3363   data[0] = (DCTELEM) ((tmp0 + tmp1 - 2 * CENTERJSAMPLE) << 5);
3364
3365   /* Odd part */
3366
3367   data[1] = (DCTELEM) ((tmp0 - tmp1) << 5);
3368 }
3369
3370
3371 /*
3372  * Perform the forward DCT on an 8x16 sample block.
3373  *
3374  * 8-point FDCT in pass 1 (rows), 16-point in pass 2 (columns).
3375  */
3376
3377 GLOBAL(void)
3378 jpeg_fdct_8x16 (DCTELEM * data, JSAMPARRAY sample_data, JDIMENSION start_col)
3379 {
3380   INT32 tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7;
3381   INT32 tmp10, tmp11, tmp12, tmp13, tmp14, tmp15, tmp16, tmp17;
3382   INT32 z1;
3383   DCTELEM workspace[DCTSIZE2];
3384   DCTELEM *dataptr;
3385   DCTELEM *wsptr;
3386   JSAMPROW elemptr;
3387   int ctr;
3388   SHIFT_TEMPS
3389
3390   /* Pass 1: process rows.
3391    * Note results are scaled up by sqrt(8) compared to a true DCT;
3392    * furthermore, we scale the results by 2**PASS1_BITS.
3393    * 8-point FDCT kernel, cK represents sqrt(2) * cos(K*pi/16).
3394    */
3395
3396   dataptr = data;
3397   ctr = 0;
3398   for (;;) {
3399     elemptr = sample_data[ctr] + start_col;
3400
3401     /* Even part per LL&M figure 1 --- note that published figure is faulty;
3402      * rotator "c1" should be "c6".
3403      */
3404
3405     tmp0 = GETJSAMPLE(elemptr[0]) + GETJSAMPLE(elemptr[7]);
3406     tmp1 = GETJSAMPLE(elemptr[1]) + GETJSAMPLE(elemptr[6]);
3407     tmp2 = GETJSAMPLE(elemptr[2]) + GETJSAMPLE(elemptr[5]);
3408     tmp3 = GETJSAMPLE(elemptr[3]) + GETJSAMPLE(elemptr[4]);
3409
3410     tmp10 = tmp0 + tmp3;
3411     tmp12 = tmp0 - tmp3;
3412     tmp11 = tmp1 + tmp2;
3413     tmp13 = tmp1 - tmp2;
3414
3415     tmp0 = GETJSAMPLE(elemptr[0]) - GETJSAMPLE(elemptr[7]);
3416     tmp1 = GETJSAMPLE(elemptr[1]) - GETJSAMPLE(elemptr[6]);
3417     tmp2 = GETJSAMPLE(elemptr[2]) - GETJSAMPLE(elemptr[5]);
3418     tmp3 = GETJSAMPLE(elemptr[3]) - GETJSAMPLE(elemptr[4]);
3419
3420     /* Apply unsigned->signed conversion */
3421     dataptr[0] = (DCTELEM) ((tmp10 + tmp11 - 8 * CENTERJSAMPLE) << PASS1_BITS);
3422     dataptr[4] = (DCTELEM) ((tmp10 - tmp11) << PASS1_BITS);
3423
3424     z1 = MULTIPLY(tmp12 + tmp13, FIX_0_541196100);   /* c6 */
3425     dataptr[2] = (DCTELEM)
3426       DESCALE(z1 + MULTIPLY(tmp12, FIX_0_765366865), /* c2-c6 */
3427               CONST_BITS-PASS1_BITS);
3428     dataptr[6] = (DCTELEM)
3429       DESCALE(z1 - MULTIPLY(tmp13, FIX_1_847759065), /* c2+c6 */
3430               CONST_BITS-PASS1_BITS);
3431
3432     /* Odd part per figure 8 --- note paper omits factor of sqrt(2).
3433      * i0..i3 in the paper are tmp0..tmp3 here.
3434      */
3435
3436     tmp12 = tmp0 + tmp2;
3437     tmp13 = tmp1 + tmp3;
3438
3439     z1 = MULTIPLY(tmp12 + tmp13, FIX_1_175875602);   /*  c3 */
3440     tmp12 = MULTIPLY(tmp12, - FIX_0_390180644);      /* -c3+c5 */
3441     tmp13 = MULTIPLY(tmp13, - FIX_1_961570560);      /* -c3-c5 */
3442     tmp12 += z1;
3443     tmp13 += z1;
3444
3445     z1 = MULTIPLY(tmp0 + tmp3, - FIX_0_899976223);   /* -c3+c7 */
3446     tmp0 = MULTIPLY(tmp0, FIX_1_501321110);          /*  c1+c3-c5-c7 */
3447     tmp3 = MULTIPLY(tmp3, FIX_0_298631336);          /* -c1+c3+c5-c7 */
3448     tmp0 += z1 + tmp12;
3449     tmp3 += z1 + tmp13;
3450
3451     z1 = MULTIPLY(tmp1 + tmp2, - FIX_2_562915447);   /* -c1-c3 */
3452     tmp1 = MULTIPLY(tmp1, FIX_3_072711026);          /*  c1+c3+c5-c7 */
3453     tmp2 = MULTIPLY(tmp2, FIX_2_053119869);          /*  c1+c3-c5+c7 */
3454     tmp1 += z1 + tmp13;
3455     tmp2 += z1 + tmp12;
3456
3457     dataptr[1] = (DCTELEM) DESCALE(tmp0, CONST_BITS-PASS1_BITS);
3458     dataptr[3] = (DCTELEM) DESCALE(tmp1, CONST_BITS-PASS1_BITS);
3459     dataptr[5] = (DCTELEM) DESCALE(tmp2, CONST_BITS-PASS1_BITS);
3460     dataptr[7] = (DCTELEM) DESCALE(tmp3, CONST_BITS-PASS1_BITS);
3461
3462     ctr++;
3463
3464     if (ctr != DCTSIZE) {
3465       if (ctr == DCTSIZE * 2)
3466         break;                  /* Done. */
3467       dataptr += DCTSIZE;       /* advance pointer to next row */
3468     } else
3469       dataptr = workspace;      /* switch pointer to extended workspace */
3470   }
3471
3472   /* Pass 2: process columns.
3473    * We remove the PASS1_BITS scaling, but leave the results scaled up
3474    * by an overall factor of 8.
3475    * We must also scale the output by 8/16 = 1/2.
3476    * 16-point FDCT kernel, cK represents sqrt(2) * cos(K*pi/32).
3477    */
3478
3479   dataptr = data;
3480   wsptr = workspace;
3481   for (ctr = DCTSIZE-1; ctr >= 0; ctr--) {
3482     /* Even part */
3483
3484     tmp0 = dataptr[DCTSIZE*0] + wsptr[DCTSIZE*7];
3485     tmp1 = dataptr[DCTSIZE*1] + wsptr[DCTSIZE*6];
3486     tmp2 = dataptr[DCTSIZE*2] + wsptr[DCTSIZE*5];
3487     tmp3 = dataptr[DCTSIZE*3] + wsptr[DCTSIZE*4];
3488     tmp4 = dataptr[DCTSIZE*4] + wsptr[DCTSIZE*3];
3489     tmp5 = dataptr[DCTSIZE*5] + wsptr[DCTSIZE*2];
3490     tmp6 = dataptr[DCTSIZE*6] + wsptr[DCTSIZE*1];
3491     tmp7 = dataptr[DCTSIZE*7] + wsptr[DCTSIZE*0];
3492
3493     tmp10 = tmp0 + tmp7;
3494     tmp14 = tmp0 - tmp7;
3495     tmp11 = tmp1 + tmp6;
3496     tmp15 = tmp1 - tmp6;
3497     tmp12 = tmp2 + tmp5;
3498     tmp16 = tmp2 - tmp5;
3499     tmp13 = tmp3 + tmp4;
3500     tmp17 = tmp3 - tmp4;
3501
3502     tmp0 = dataptr[DCTSIZE*0] - wsptr[DCTSIZE*7];
3503     tmp1 = dataptr[DCTSIZE*1] - wsptr[DCTSIZE*6];
3504     tmp2 = dataptr[DCTSIZE*2] - wsptr[DCTSIZE*5];
3505     tmp3 = dataptr[DCTSIZE*3] - wsptr[DCTSIZE*4];
3506     tmp4 = dataptr[DCTSIZE*4] - wsptr[DCTSIZE*3];
3507     tmp5 = dataptr[DCTSIZE*5] - wsptr[DCTSIZE*2];
3508     tmp6 = dataptr[DCTSIZE*6] - wsptr[DCTSIZE*1];
3509     tmp7 = dataptr[DCTSIZE*7] - wsptr[DCTSIZE*0];
3510
3511     dataptr[DCTSIZE*0] = (DCTELEM)
3512       DESCALE(tmp10 + tmp11 + tmp12 + tmp13, PASS1_BITS+1);
3513     dataptr[DCTSIZE*4] = (DCTELEM)
3514       DESCALE(MULTIPLY(tmp10 - tmp13, FIX(1.306562965)) + /* c4[16] = c2[8] */
3515               MULTIPLY(tmp11 - tmp12, FIX_0_541196100),   /* c12[16] = c6[8] */
3516               CONST_BITS+PASS1_BITS+1);
3517
3518     tmp10 = MULTIPLY(tmp17 - tmp15, FIX(0.275899379)) +   /* c14[16] = c7[8] */
3519             MULTIPLY(tmp14 - tmp16, FIX(1.387039845));    /* c2[16] = c1[8] */
3520
3521     dataptr[DCTSIZE*2] = (DCTELEM)
3522       DESCALE(tmp10 + MULTIPLY(tmp15, FIX(1.451774982))   /* c6+c14 */
3523               + MULTIPLY(tmp16, FIX(2.172734804)),        /* c2+c10 */
3524               CONST_BITS+PASS1_BITS+1);
3525     dataptr[DCTSIZE*6] = (DCTELEM)
3526       DESCALE(tmp10 - MULTIPLY(tmp14, FIX(0.211164243))   /* c2-c6 */
3527               - MULTIPLY(tmp17, FIX(1.061594338)),        /* c10+c14 */
3528               CONST_BITS+PASS1_BITS+1);
3529
3530     /* Odd part */
3531
3532     tmp11 = MULTIPLY(tmp0 + tmp1, FIX(1.353318001)) +         /* c3 */
3533             MULTIPLY(tmp6 - tmp7, FIX(0.410524528));          /* c13 */
3534     tmp12 = MULTIPLY(tmp0 + tmp2, FIX(1.247225013)) +         /* c5 */
3535             MULTIPLY(tmp5 + tmp7, FIX(0.666655658));          /* c11 */
3536     tmp13 = MULTIPLY(tmp0 + tmp3, FIX(1.093201867)) +         /* c7 */
3537             MULTIPLY(tmp4 - tmp7, FIX(0.897167586));          /* c9 */
3538     tmp14 = MULTIPLY(tmp1 + tmp2, FIX(0.138617169)) +         /* c15 */
3539             MULTIPLY(tmp6 - tmp5, FIX(1.407403738));          /* c1 */
3540     tmp15 = MULTIPLY(tmp1 + tmp3, - FIX(0.666655658)) +       /* -c11 */
3541             MULTIPLY(tmp4 + tmp6, - FIX(1.247225013));        /* -c5 */
3542     tmp16 = MULTIPLY(tmp2 + tmp3, - FIX(1.353318001)) +       /* -c3 */
3543             MULTIPLY(tmp5 - tmp4, FIX(0.410524528));          /* c13 */
3544     tmp10 = tmp11 + tmp12 + tmp13 -
3545             MULTIPLY(tmp0, FIX(2.286341144)) +                /* c7+c5+c3-c1 */
3546             MULTIPLY(tmp7, FIX(0.779653625));                 /* c15+c13-c11+c9 */
3547     tmp11 += tmp14 + tmp15 + MULTIPLY(tmp1, FIX(0.071888074)) /* c9-c3-c15+c11 */
3548              - MULTIPLY(tmp6, FIX(1.663905119));              /* c7+c13+c1-c5 */
3549     tmp12 += tmp14 + tmp16 - MULTIPLY(tmp2, FIX(1.125726048)) /* c7+c5+c15-c3 */
3550              + MULTIPLY(tmp5, FIX(1.227391138));              /* c9-c11+c1-c13 */
3551     tmp13 += tmp15 + tmp16 + MULTIPLY(tmp3, FIX(1.065388962)) /* c15+c3+c11-c7 */
3552              + MULTIPLY(tmp4, FIX(2.167985692));              /* c1+c13+c5-c9 */
3553
3554     dataptr[DCTSIZE*1] = (DCTELEM) DESCALE(tmp10, CONST_BITS+PASS1_BITS+1);
3555     dataptr[DCTSIZE*3] = (DCTELEM) DESCALE(tmp11, CONST_BITS+PASS1_BITS+1);
3556     dataptr[DCTSIZE*5] = (DCTELEM) DESCALE(tmp12, CONST_BITS+PASS1_BITS+1);
3557     dataptr[DCTSIZE*7] = (DCTELEM) DESCALE(tmp13, CONST_BITS+PASS1_BITS+1);
3558
3559     dataptr++;                  /* advance pointer to next column */
3560     wsptr++;                    /* advance pointer to next column */
3561   }
3562 }
3563
3564
3565 /*
3566  * Perform the forward DCT on a 7x14 sample block.
3567  *
3568  * 7-point FDCT in pass 1 (rows), 14-point in pass 2 (columns).
3569  */
3570
3571 GLOBAL(void)
3572 jpeg_fdct_7x14 (DCTELEM * data, JSAMPARRAY sample_data, JDIMENSION start_col)
3573 {
3574   INT32 tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6;
3575   INT32 tmp10, tmp11, tmp12, tmp13, tmp14, tmp15, tmp16;
3576   INT32 z1, z2, z3;
3577   DCTELEM workspace[8*6];
3578   DCTELEM *dataptr;
3579   DCTELEM *wsptr;
3580   JSAMPROW elemptr;
3581   int ctr;
3582   SHIFT_TEMPS
3583
3584   /* Pre-zero output coefficient block. */
3585   MEMZERO(data, SIZEOF(DCTELEM) * DCTSIZE2);
3586
3587   /* Pass 1: process rows.
3588    * Note results are scaled up by sqrt(8) compared to a true DCT;
3589    * furthermore, we scale the results by 2**PASS1_BITS.
3590    * 7-point FDCT kernel, cK represents sqrt(2) * cos(K*pi/14).
3591    */
3592
3593   dataptr = data;
3594   ctr = 0;
3595   for (;;) {
3596     elemptr = sample_data[ctr] + start_col;
3597
3598     /* Even part */
3599
3600     tmp0 = GETJSAMPLE(elemptr[0]) + GETJSAMPLE(elemptr[6]);
3601     tmp1 = GETJSAMPLE(elemptr[1]) + GETJSAMPLE(elemptr[5]);
3602     tmp2 = GETJSAMPLE(elemptr[2]) + GETJSAMPLE(elemptr[4]);
3603     tmp3 = GETJSAMPLE(elemptr[3]);
3604
3605     tmp10 = GETJSAMPLE(elemptr[0]) - GETJSAMPLE(elemptr[6]);
3606     tmp11 = GETJSAMPLE(elemptr[1]) - GETJSAMPLE(elemptr[5]);
3607     tmp12 = GETJSAMPLE(elemptr[2]) - GETJSAMPLE(elemptr[4]);
3608
3609     z1 = tmp0 + tmp2;
3610     /* Apply unsigned->signed conversion */
3611     dataptr[0] = (DCTELEM)
3612       ((z1 + tmp1 + tmp3 - 7 * CENTERJSAMPLE) << PASS1_BITS);
3613     tmp3 += tmp3;
3614     z1 -= tmp3;
3615     z1 -= tmp3;
3616     z1 = MULTIPLY(z1, FIX(0.353553391));                /* (c2+c6-c4)/2 */
3617     z2 = MULTIPLY(tmp0 - tmp2, FIX(0.920609002));       /* (c2+c4-c6)/2 */
3618     z3 = MULTIPLY(tmp1 - tmp2, FIX(0.314692123));       /* c6 */
3619     dataptr[2] = (DCTELEM) DESCALE(z1 + z2 + z3, CONST_BITS-PASS1_BITS);
3620     z1 -= z2;
3621     z2 = MULTIPLY(tmp0 - tmp1, FIX(0.881747734));       /* c4 */
3622     dataptr[4] = (DCTELEM)
3623       DESCALE(z2 + z3 - MULTIPLY(tmp1 - tmp3, FIX(0.707106781)), /* c2+c6-c4 */
3624               CONST_BITS-PASS1_BITS);
3625     dataptr[6] = (DCTELEM) DESCALE(z1 + z2, CONST_BITS-PASS1_BITS);
3626
3627     /* Odd part */
3628
3629     tmp1 = MULTIPLY(tmp10 + tmp11, FIX(0.935414347));   /* (c3+c1-c5)/2 */
3630     tmp2 = MULTIPLY(tmp10 - tmp11, FIX(0.170262339));   /* (c3+c5-c1)/2 */
3631     tmp0 = tmp1 - tmp2;
3632     tmp1 += tmp2;
3633     tmp2 = MULTIPLY(tmp11 + tmp12, - FIX(1.378756276)); /* -c1 */
3634     tmp1 += tmp2;
3635     tmp3 = MULTIPLY(tmp10 + tmp12, FIX(0.613604268));   /* c5 */
3636     tmp0 += tmp3;
3637     tmp2 += tmp3 + MULTIPLY(tmp12, FIX(1.870828693));   /* c3+c1-c5 */
3638
3639     dataptr[1] = (DCTELEM) DESCALE(tmp0, CONST_BITS-PASS1_BITS);
3640     dataptr[3] = (DCTELEM) DESCALE(tmp1, CONST_BITS-PASS1_BITS);
3641     dataptr[5] = (DCTELEM) DESCALE(tmp2, CONST_BITS-PASS1_BITS);
3642
3643     ctr++;
3644
3645     if (ctr != DCTSIZE) {
3646       if (ctr == 14)
3647         break;                  /* Done. */
3648       dataptr += DCTSIZE;       /* advance pointer to next row */
3649     } else
3650       dataptr = workspace;      /* switch pointer to extended workspace */
3651   }
3652
3653   /* Pass 2: process columns.
3654    * We remove the PASS1_BITS scaling, but leave the results scaled up
3655    * by an overall factor of 8.
3656    * We must also scale the output by (8/7)*(8/14) = 32/49, which we
3657    * fold into the constant multipliers:
3658    * 14-point FDCT kernel, cK represents sqrt(2) * cos(K*pi/28) * 32/49.
3659    */
3660
3661   dataptr = data;
3662   wsptr = workspace;
3663   for (ctr = 0; ctr < 7; ctr++) {
3664     /* Even part */
3665
3666     tmp0 = dataptr[DCTSIZE*0] + wsptr[DCTSIZE*5];
3667     tmp1 = dataptr[DCTSIZE*1] + wsptr[DCTSIZE*4];
3668     tmp2 = dataptr[DCTSIZE*2] + wsptr[DCTSIZE*3];
3669     tmp13 = dataptr[DCTSIZE*3] + wsptr[DCTSIZE*2];
3670     tmp4 = dataptr[DCTSIZE*4] + wsptr[DCTSIZE*1];
3671     tmp5 = dataptr[DCTSIZE*5] + wsptr[DCTSIZE*0];
3672     tmp6 = dataptr[DCTSIZE*6] + dataptr[DCTSIZE*7];
3673
3674     tmp10 = tmp0 + tmp6;
3675     tmp14 = tmp0 - tmp6;
3676     tmp11 = tmp1 + tmp5;
3677     tmp15 = tmp1 - tmp5;
3678     tmp12 = tmp2 + tmp4;
3679     tmp16 = tmp2 - tmp4;
3680
3681     tmp0 = dataptr[DCTSIZE*0] - wsptr[DCTSIZE*5];
3682     tmp1 = dataptr[DCTSIZE*1] - wsptr[DCTSIZE*4];
3683     tmp2 = dataptr[DCTSIZE*2] - wsptr[DCTSIZE*3];
3684     tmp3 = dataptr[DCTSIZE*3] - wsptr[DCTSIZE*2];
3685     tmp4 = dataptr[DCTSIZE*4] - wsptr[DCTSIZE*1];
3686     tmp5 = dataptr[DCTSIZE*5] - wsptr[DCTSIZE*0];
3687     tmp6 = dataptr[DCTSIZE*6] - dataptr[DCTSIZE*7];
3688
3689     dataptr[DCTSIZE*0] = (DCTELEM)
3690       DESCALE(MULTIPLY(tmp10 + tmp11 + tmp12 + tmp13,
3691                        FIX(0.653061224)),                 /* 32/49 */
3692               CONST_BITS+PASS1_BITS);
3693     tmp13 += tmp13;
3694     dataptr[DCTSIZE*4] = (DCTELEM)
3695       DESCALE(MULTIPLY(tmp10 - tmp13, FIX(0.832106052)) + /* c4 */
3696               MULTIPLY(tmp11 - tmp13, FIX(0.205513223)) - /* c12 */
3697               MULTIPLY(tmp12 - tmp13, FIX(0.575835255)),  /* c8 */
3698               CONST_BITS+PASS1_BITS);
3699
3700     tmp10 = MULTIPLY(tmp14 + tmp15, FIX(0.722074570));    /* c6 */
3701
3702     dataptr[DCTSIZE*2] = (DCTELEM)
3703       DESCALE(tmp10 + MULTIPLY(tmp14, FIX(0.178337691))   /* c2-c6 */
3704               + MULTIPLY(tmp16, FIX(0.400721155)),        /* c10 */
3705               CONST_BITS+PASS1_BITS);
3706     dataptr[DCTSIZE*6] = (DCTELEM)
3707       DESCALE(tmp10 - MULTIPLY(tmp15, FIX(1.122795725))   /* c6+c10 */
3708               - MULTIPLY(tmp16, FIX(0.900412262)),        /* c2 */
3709               CONST_BITS+PASS1_BITS);
3710
3711     /* Odd part */
3712
3713     tmp10 = tmp1 + tmp2;
3714     tmp11 = tmp5 - tmp4;
3715     dataptr[DCTSIZE*7] = (DCTELEM)
3716       DESCALE(MULTIPLY(tmp0 - tmp10 + tmp3 - tmp11 - tmp6,
3717                        FIX(0.653061224)),                 /* 32/49 */
3718               CONST_BITS+PASS1_BITS);
3719     tmp3  = MULTIPLY(tmp3 , FIX(0.653061224));            /* 32/49 */
3720     tmp10 = MULTIPLY(tmp10, - FIX(0.103406812));          /* -c13 */
3721     tmp11 = MULTIPLY(tmp11, FIX(0.917760839));            /* c1 */
3722     tmp10 += tmp11 - tmp3;
3723     tmp11 = MULTIPLY(tmp0 + tmp2, FIX(0.782007410)) +     /* c5 */
3724             MULTIPLY(tmp4 + tmp6, FIX(0.491367823));      /* c9 */
3725     dataptr[DCTSIZE*5] = (DCTELEM)
3726       DESCALE(tmp10 + tmp11 - MULTIPLY(tmp2, FIX(1.550341076)) /* c3+c5-c13 */
3727               + MULTIPLY(tmp4, FIX(0.731428202)),         /* c1+c11-c9 */
3728               CONST_BITS+PASS1_BITS);
3729     tmp12 = MULTIPLY(tmp0 + tmp1, FIX(0.871740478)) +     /* c3 */
3730             MULTIPLY(tmp5 - tmp6, FIX(0.305035186));      /* c11 */
3731     dataptr[DCTSIZE*3] = (DCTELEM)
3732       DESCALE(tmp10 + tmp12 - MULTIPLY(tmp1, FIX(0.276965844)) /* c3-c9-c13 */
3733               - MULTIPLY(tmp5, FIX(2.004803435)),         /* c1+c5+c11 */
3734               CONST_BITS+PASS1_BITS);
3735     dataptr[DCTSIZE*1] = (DCTELEM)
3736       DESCALE(tmp11 + tmp12 + tmp3
3737               - MULTIPLY(tmp0, FIX(0.735987049))          /* c3+c5-c1 */
3738               - MULTIPLY(tmp6, FIX(0.082925825)),         /* c9-c11-c13 */
3739               CONST_BITS+PASS1_BITS);
3740
3741     dataptr++;                  /* advance pointer to next column */
3742     wsptr++;                    /* advance pointer to next column */
3743   }
3744 }
3745
3746
3747 /*
3748  * Perform the forward DCT on a 6x12 sample block.
3749  *
3750  * 6-point FDCT in pass 1 (rows), 12-point in pass 2 (columns).
3751  */
3752
3753 GLOBAL(void)
3754 jpeg_fdct_6x12 (DCTELEM * data, JSAMPARRAY sample_data, JDIMENSION start_col)
3755 {
3756   INT32 tmp0, tmp1, tmp2, tmp3, tmp4, tmp5;
3757   INT32 tmp10, tmp11, tmp12, tmp13, tmp14, tmp15;
3758   DCTELEM workspace[8*4];
3759   DCTELEM *dataptr;
3760   DCTELEM *wsptr;
3761   JSAMPROW elemptr;
3762   int ctr;
3763   SHIFT_TEMPS
3764
3765   /* Pre-zero output coefficient block. */
3766   MEMZERO(data, SIZEOF(DCTELEM) * DCTSIZE2);
3767
3768   /* Pass 1: process rows.
3769    * Note results are scaled up by sqrt(8) compared to a true DCT;
3770    * furthermore, we scale the results by 2**PASS1_BITS.
3771    * 6-point FDCT kernel, cK represents sqrt(2) * cos(K*pi/12).
3772    */
3773
3774   dataptr = data;
3775   ctr = 0;
3776   for (;;) {
3777     elemptr = sample_data[ctr] + start_col;
3778
3779     /* Even part */
3780
3781     tmp0 = GETJSAMPLE(elemptr[0]) + GETJSAMPLE(elemptr[5]);
3782     tmp11 = GETJSAMPLE(elemptr[1]) + GETJSAMPLE(elemptr[4]);
3783     tmp2 = GETJSAMPLE(elemptr[2]) + GETJSAMPLE(elemptr[3]);
3784
3785     tmp10 = tmp0 + tmp2;
3786     tmp12 = tmp0 - tmp2;
3787
3788     tmp0 = GETJSAMPLE(elemptr[0]) - GETJSAMPLE(elemptr[5]);
3789     tmp1 = GETJSAMPLE(elemptr[1]) - GETJSAMPLE(elemptr[4]);
3790     tmp2 = GETJSAMPLE(elemptr[2]) - GETJSAMPLE(elemptr[3]);
3791
3792     /* Apply unsigned->signed conversion */
3793     dataptr[0] = (DCTELEM)
3794       ((tmp10 + tmp11 - 6 * CENTERJSAMPLE) << PASS1_BITS);
3795     dataptr[2] = (DCTELEM)
3796       DESCALE(MULTIPLY(tmp12, FIX(1.224744871)),                 /* c2 */
3797               CONST_BITS-PASS1_BITS);
3798     dataptr[4] = (DCTELEM)
3799       DESCALE(MULTIPLY(tmp10 - tmp11 - tmp11, FIX(0.707106781)), /* c4 */
3800               CONST_BITS-PASS1_BITS);
3801
3802     /* Odd part */
3803
3804     tmp10 = DESCALE(MULTIPLY(tmp0 + tmp2, FIX(0.366025404)),     /* c5 */
3805                     CONST_BITS-PASS1_BITS);
3806
3807     dataptr[1] = (DCTELEM) (tmp10 + ((tmp0 + tmp1) << PASS1_BITS));
3808     dataptr[3] = (DCTELEM) ((tmp0 - tmp1 - tmp2) << PASS1_BITS);
3809     dataptr[5] = (DCTELEM) (tmp10 + ((tmp2 - tmp1) << PASS1_BITS));
3810
3811     ctr++;
3812
3813     if (ctr != DCTSIZE) {
3814       if (ctr == 12)
3815         break;                  /* Done. */
3816       dataptr += DCTSIZE;       /* advance pointer to next row */
3817     } else
3818       dataptr = workspace;      /* switch pointer to extended workspace */
3819   }
3820
3821   /* Pass 2: process columns.
3822    * We remove the PASS1_BITS scaling, but leave the results scaled up
3823    * by an overall factor of 8.
3824    * We must also scale the output by (8/6)*(8/12) = 8/9, which we
3825    * fold into the constant multipliers:
3826    * 12-point FDCT kernel, cK represents sqrt(2) * cos(K*pi/24) * 8/9.
3827    */
3828
3829   dataptr = data;
3830   wsptr = workspace;
3831   for (ctr = 0; ctr < 6; ctr++) {
3832     /* Even part */
3833
3834     tmp0 = dataptr[DCTSIZE*0] + wsptr[DCTSIZE*3];
3835     tmp1 = dataptr[DCTSIZE*1] + wsptr[DCTSIZE*2];
3836     tmp2 = dataptr[DCTSIZE*2] + wsptr[DCTSIZE*1];
3837     tmp3 = dataptr[DCTSIZE*3] + wsptr[DCTSIZE*0];
3838     tmp4 = dataptr[DCTSIZE*4] + dataptr[DCTSIZE*7];
3839     tmp5 = dataptr[DCTSIZE*5] + dataptr[DCTSIZE*6];
3840
3841     tmp10 = tmp0 + tmp5;
3842     tmp13 = tmp0 - tmp5;
3843     tmp11 = tmp1 + tmp4;
3844     tmp14 = tmp1 - tmp4;
3845     tmp12 = tmp2 + tmp3;
3846     tmp15 = tmp2 - tmp3;
3847
3848     tmp0 = dataptr[DCTSIZE*0] - wsptr[DCTSIZE*3];
3849     tmp1 = dataptr[DCTSIZE*1] - wsptr[DCTSIZE*2];
3850     tmp2 = dataptr[DCTSIZE*2] - wsptr[DCTSIZE*1];
3851     tmp3 = dataptr[DCTSIZE*3] - wsptr[DCTSIZE*0];
3852     tmp4 = dataptr[DCTSIZE*4] - dataptr[DCTSIZE*7];
3853     tmp5 = dataptr[DCTSIZE*5] - dataptr[DCTSIZE*6];
3854
3855     dataptr[DCTSIZE*0] = (DCTELEM)
3856       DESCALE(MULTIPLY(tmp10 + tmp11 + tmp12, FIX(0.888888889)), /* 8/9 */
3857               CONST_BITS+PASS1_BITS);
3858     dataptr[DCTSIZE*6] = (DCTELEM)
3859       DESCALE(MULTIPLY(tmp13 - tmp14 - tmp15, FIX(0.888888889)), /* 8/9 */
3860               CONST_BITS+PASS1_BITS);
3861     dataptr[DCTSIZE*4] = (DCTELEM)
3862       DESCALE(MULTIPLY(tmp10 - tmp12, FIX(1.088662108)),         /* c4 */
3863               CONST_BITS+PASS1_BITS);
3864     dataptr[DCTSIZE*2] = (DCTELEM)
3865       DESCALE(MULTIPLY(tmp14 - tmp15, FIX(0.888888889)) +        /* 8/9 */
3866               MULTIPLY(tmp13 + tmp15, FIX(1.214244803)),         /* c2 */
3867               CONST_BITS+PASS1_BITS);
3868
3869     /* Odd part */
3870
3871     tmp10 = MULTIPLY(tmp1 + tmp4, FIX(0.481063200));   /* c9 */
3872     tmp14 = tmp10 + MULTIPLY(tmp1, FIX(0.680326102));  /* c3-c9 */
3873     tmp15 = tmp10 - MULTIPLY(tmp4, FIX(1.642452502));  /* c3+c9 */
3874     tmp12 = MULTIPLY(tmp0 + tmp2, FIX(0.997307603));   /* c5 */
3875     tmp13 = MULTIPLY(tmp0 + tmp3, FIX(0.765261039));   /* c7 */
3876     tmp10 = tmp12 + tmp13 + tmp14 - MULTIPLY(tmp0, FIX(0.516244403)) /* c5+c7-c1 */
3877             + MULTIPLY(tmp5, FIX(0.164081699));        /* c11 */
3878     tmp11 = MULTIPLY(tmp2 + tmp3, - FIX(0.164081699)); /* -c11 */
3879     tmp12 += tmp11 - tmp15 - MULTIPLY(tmp2, FIX(2.079550144)) /* c1+c5-c11 */
3880             + MULTIPLY(tmp5, FIX(0.765261039));        /* c7 */
3881     tmp13 += tmp11 - tmp14 + MULTIPLY(tmp3, FIX(0.645144899)) /* c1+c11-c7 */
3882             - MULTIPLY(tmp5, FIX(0.997307603));        /* c5 */
3883     tmp11 = tmp15 + MULTIPLY(tmp0 - tmp3, FIX(1.161389302)) /* c3 */
3884             - MULTIPLY(tmp2 + tmp5, FIX(0.481063200)); /* c9 */
3885
3886     dataptr[DCTSIZE*1] = (DCTELEM) DESCALE(tmp10, CONST_BITS+PASS1_BITS);
3887     dataptr[DCTSIZE*3] = (DCTELEM) DESCALE(tmp11, CONST_BITS+PASS1_BITS);
3888     dataptr[DCTSIZE*5] = (DCTELEM) DESCALE(tmp12, CONST_BITS+PASS1_BITS);
3889     dataptr[DCTSIZE*7] = (DCTELEM) DESCALE(tmp13, CONST_BITS+PASS1_BITS);
3890
3891     dataptr++;                  /* advance pointer to next column */
3892     wsptr++;                    /* advance pointer to next column */
3893   }
3894 }
3895
3896
3897 /*
3898  * Perform the forward DCT on a 5x10 sample block.
3899  *
3900  * 5-point FDCT in pass 1 (rows), 10-point in pass 2 (columns).
3901  */
3902
3903 GLOBAL(void)
3904 jpeg_fdct_5x10 (DCTELEM * data, JSAMPARRAY sample_data, JDIMENSION start_col)
3905 {
3906   INT32 tmp0, tmp1, tmp2, tmp3, tmp4;
3907   INT32 tmp10, tmp11, tmp12, tmp13, tmp14;
3908   DCTELEM workspace[8*2];
3909   DCTELEM *dataptr;
3910   DCTELEM *wsptr;
3911   JSAMPROW elemptr;
3912   int ctr;
3913   SHIFT_TEMPS
3914
3915   /* Pre-zero output coefficient block. */
3916   MEMZERO(data, SIZEOF(DCTELEM) * DCTSIZE2);
3917
3918   /* Pass 1: process rows.
3919    * Note results are scaled up by sqrt(8) compared to a true DCT;
3920    * furthermore, we scale the results by 2**PASS1_BITS.
3921    * 5-point FDCT kernel, cK represents sqrt(2) * cos(K*pi/10).
3922    */
3923
3924   dataptr = data;
3925   ctr = 0;
3926   for (;;) {
3927     elemptr = sample_data[ctr] + start_col;
3928
3929     /* Even part */
3930
3931     tmp0 = GETJSAMPLE(elemptr[0]) + GETJSAMPLE(elemptr[4]);
3932     tmp1 = GETJSAMPLE(elemptr[1]) + GETJSAMPLE(elemptr[3]);
3933     tmp2 = GETJSAMPLE(elemptr[2]);
3934
3935     tmp10 = tmp0 + tmp1;
3936     tmp11 = tmp0 - tmp1;
3937
3938     tmp0 = GETJSAMPLE(elemptr[0]) - GETJSAMPLE(elemptr[4]);
3939     tmp1 = GETJSAMPLE(elemptr[1]) - GETJSAMPLE(elemptr[3]);
3940
3941     /* Apply unsigned->signed conversion */
3942     dataptr[0] = (DCTELEM)
3943       ((tmp10 + tmp2 - 5 * CENTERJSAMPLE) << PASS1_BITS);
3944     tmp11 = MULTIPLY(tmp11, FIX(0.790569415));          /* (c2+c4)/2 */
3945     tmp10 -= tmp2 << 2;
3946     tmp10 = MULTIPLY(tmp10, FIX(0.353553391));          /* (c2-c4)/2 */
3947     dataptr[2] = (DCTELEM) DESCALE(tmp11 + tmp10, CONST_BITS-PASS1_BITS);
3948     dataptr[4] = (DCTELEM) DESCALE(tmp11 - tmp10, CONST_BITS-PASS1_BITS);
3949
3950     /* Odd part */
3951
3952     tmp10 = MULTIPLY(tmp0 + tmp1, FIX(0.831253876));    /* c3 */
3953
3954     dataptr[1] = (DCTELEM)
3955       DESCALE(tmp10 + MULTIPLY(tmp0, FIX(0.513743148)), /* c1-c3 */
3956               CONST_BITS-PASS1_BITS);
3957     dataptr[3] = (DCTELEM)
3958       DESCALE(tmp10 - MULTIPLY(tmp1, FIX(2.176250899)), /* c1+c3 */
3959               CONST_BITS-PASS1_BITS);
3960
3961     ctr++;
3962
3963     if (ctr != DCTSIZE) {
3964       if (ctr == 10)
3965         break;                  /* Done. */
3966       dataptr += DCTSIZE;       /* advance pointer to next row */
3967     } else
3968       dataptr = workspace;      /* switch pointer to extended workspace */
3969   }
3970
3971   /* Pass 2: process columns.
3972    * We remove the PASS1_BITS scaling, but leave the results scaled up
3973    * by an overall factor of 8.
3974    * We must also scale the output by (8/5)*(8/10) = 32/25, which we
3975    * fold into the constant multipliers:
3976    * 10-point FDCT kernel, cK represents sqrt(2) * cos(K*pi/20) * 32/25.
3977    */
3978
3979   dataptr = data;
3980   wsptr = workspace;
3981   for (ctr = 0; ctr < 5; ctr++) {
3982     /* Even part */
3983
3984     tmp0 = dataptr[DCTSIZE*0] + wsptr[DCTSIZE*1];
3985     tmp1 = dataptr[DCTSIZE*1] + wsptr[DCTSIZE*0];
3986     tmp12 = dataptr[DCTSIZE*2] + dataptr[DCTSIZE*7];
3987     tmp3 = dataptr[DCTSIZE*3] + dataptr[DCTSIZE*6];
3988     tmp4 = dataptr[DCTSIZE*4] + dataptr[DCTSIZE*5];
3989
3990     tmp10 = tmp0 + tmp4;
3991     tmp13 = tmp0 - tmp4;
3992     tmp11 = tmp1 + tmp3;
3993     tmp14 = tmp1 - tmp3;
3994
3995     tmp0 = dataptr[DCTSIZE*0] - wsptr[DCTSIZE*1];
3996     tmp1 = dataptr[DCTSIZE*1] - wsptr[DCTSIZE*0];
3997     tmp2 = dataptr[DCTSIZE*2] - dataptr[DCTSIZE*7];
3998     tmp3 = dataptr[DCTSIZE*3] - dataptr[DCTSIZE*6];
3999     tmp4 = dataptr[DCTSIZE*4] - dataptr[DCTSIZE*5];
4000
4001     dataptr[DCTSIZE*0] = (DCTELEM)
4002       DESCALE(MULTIPLY(tmp10 + tmp11 + tmp12, FIX(1.28)), /* 32/25 */
4003               CONST_BITS+PASS1_BITS);
4004     tmp12 += tmp12;
4005     dataptr[DCTSIZE*4] = (DCTELEM)
4006       DESCALE(MULTIPLY(tmp10 - tmp12, FIX(1.464477191)) - /* c4 */
4007               MULTIPLY(tmp11 - tmp12, FIX(0.559380511)),  /* c8 */
4008               CONST_BITS+PASS1_BITS);
4009     tmp10 = MULTIPLY(tmp13 + tmp14, FIX(1.064004961));    /* c6 */
4010     dataptr[DCTSIZE*2] = (DCTELEM)
4011       DESCALE(tmp10 + MULTIPLY(tmp13, FIX(0.657591230)),  /* c2-c6 */
4012               CONST_BITS+PASS1_BITS);
4013     dataptr[DCTSIZE*6] = (DCTELEM)
4014       DESCALE(tmp10 - MULTIPLY(tmp14, FIX(2.785601151)),  /* c2+c6 */
4015               CONST_BITS+PASS1_BITS);
4016
4017     /* Odd part */
4018
4019     tmp10 = tmp0 + tmp4;
4020     tmp11 = tmp1 - tmp3;
4021     dataptr[DCTSIZE*5] = (DCTELEM)
4022       DESCALE(MULTIPLY(tmp10 - tmp11 - tmp2, FIX(1.28)),  /* 32/25 */
4023               CONST_BITS+PASS1_BITS);
4024     tmp2 = MULTIPLY(tmp2, FIX(1.28));                     /* 32/25 */
4025     dataptr[DCTSIZE*1] = (DCTELEM)
4026       DESCALE(MULTIPLY(tmp0, FIX(1.787906876)) +          /* c1 */
4027               MULTIPLY(tmp1, FIX(1.612894094)) + tmp2 +   /* c3 */
4028               MULTIPLY(tmp3, FIX(0.821810588)) +          /* c7 */
4029               MULTIPLY(tmp4, FIX(0.283176630)),           /* c9 */
4030               CONST_BITS+PASS1_BITS);
4031     tmp12 = MULTIPLY(tmp0 - tmp4, FIX(1.217352341)) -     /* (c3+c7)/2 */
4032             MULTIPLY(tmp1 + tmp3, FIX(0.752365123));      /* (c1-c9)/2 */
4033     tmp13 = MULTIPLY(tmp10 + tmp11, FIX(0.395541753)) +   /* (c3-c7)/2 */
4034             MULTIPLY(tmp11, FIX(0.64)) - tmp2;            /* 16/25 */
4035     dataptr[DCTSIZE*3] = (DCTELEM) DESCALE(tmp12 + tmp13, CONST_BITS+PASS1_BITS);
4036     dataptr[DCTSIZE*7] = (DCTELEM) DESCALE(tmp12 - tmp13, CONST_BITS+PASS1_BITS);
4037
4038     dataptr++;                  /* advance pointer to next column */
4039     wsptr++;                    /* advance pointer to next column */
4040   }
4041 }
4042
4043
4044 /*
4045  * Perform the forward DCT on a 4x8 sample block.
4046  *
4047  * 4-point FDCT in pass 1 (rows), 8-point in pass 2 (columns).
4048  */
4049
4050 GLOBAL(void)
4051 jpeg_fdct_4x8 (DCTELEM * data, JSAMPARRAY sample_data, JDIMENSION start_col)
4052 {
4053   INT32 tmp0, tmp1, tmp2, tmp3;
4054   INT32 tmp10, tmp11, tmp12, tmp13;
4055   INT32 z1;
4056   DCTELEM *dataptr;
4057   JSAMPROW elemptr;
4058   int ctr;
4059   SHIFT_TEMPS
4060
4061   /* Pre-zero output coefficient block. */
4062   MEMZERO(data, SIZEOF(DCTELEM) * DCTSIZE2);
4063
4064   /* Pass 1: process rows.
4065    * Note results are scaled up by sqrt(8) compared to a true DCT;
4066    * furthermore, we scale the results by 2**PASS1_BITS.
4067    * We must also scale the output by 8/4 = 2, which we add here.
4068    * 4-point FDCT kernel,
4069    * cK represents sqrt(2) * cos(K*pi/16) [refers to 8-point FDCT].
4070    */
4071
4072   dataptr = data;
4073   for (ctr = 0; ctr < DCTSIZE; ctr++) {
4074     elemptr = sample_data[ctr] + start_col;
4075
4076     /* Even part */
4077
4078     tmp0 = GETJSAMPLE(elemptr[0]) + GETJSAMPLE(elemptr[3]);
4079     tmp1 = GETJSAMPLE(elemptr[1]) + GETJSAMPLE(elemptr[2]);
4080
4081     tmp10 = GETJSAMPLE(elemptr[0]) - GETJSAMPLE(elemptr[3]);
4082     tmp11 = GETJSAMPLE(elemptr[1]) - GETJSAMPLE(elemptr[2]);
4083
4084     /* Apply unsigned->signed conversion */
4085     dataptr[0] = (DCTELEM)
4086       ((tmp0 + tmp1 - 4 * CENTERJSAMPLE) << (PASS1_BITS+1));
4087     dataptr[2] = (DCTELEM) ((tmp0 - tmp1) << (PASS1_BITS+1));
4088
4089     /* Odd part */
4090
4091     tmp0 = MULTIPLY(tmp10 + tmp11, FIX_0_541196100);       /* c6 */
4092     /* Add fudge factor here for final descale. */
4093     tmp0 += ONE << (CONST_BITS-PASS1_BITS-2);
4094
4095     dataptr[1] = (DCTELEM)
4096       RIGHT_SHIFT(tmp0 + MULTIPLY(tmp10, FIX_0_765366865), /* c2-c6 */
4097                   CONST_BITS-PASS1_BITS-1);
4098     dataptr[3] = (DCTELEM)
4099       RIGHT_SHIFT(tmp0 - MULTIPLY(tmp11, FIX_1_847759065), /* c2+c6 */
4100                   CONST_BITS-PASS1_BITS-1);
4101
4102     dataptr += DCTSIZE;         /* advance pointer to next row */
4103   }
4104
4105   /* Pass 2: process columns.
4106    * We remove the PASS1_BITS scaling, but leave the results scaled up
4107    * by an overall factor of 8.
4108    * 8-point FDCT kernel, cK represents sqrt(2) * cos(K*pi/16).
4109    */
4110
4111   dataptr = data;
4112   for (ctr = 0; ctr < 4; ctr++) {
4113     /* Even part per LL&M figure 1 --- note that published figure is faulty;
4114      * rotator "c1" should be "c6".
4115      */
4116
4117     tmp0 = dataptr[DCTSIZE*0] + dataptr[DCTSIZE*7];
4118     tmp1 = dataptr[DCTSIZE*1] + dataptr[DCTSIZE*6];
4119     tmp2 = dataptr[DCTSIZE*2] + dataptr[DCTSIZE*5];
4120     tmp3 = dataptr[DCTSIZE*3] + dataptr[DCTSIZE*4];
4121
4122     /* Add fudge factor here for final descale. */
4123     tmp10 = tmp0 + tmp3 + (ONE << (PASS1_BITS-1));
4124     tmp12 = tmp0 - tmp3;
4125     tmp11 = tmp1 + tmp2;
4126     tmp13 = tmp1 - tmp2;
4127
4128     tmp0 = dataptr[DCTSIZE*0] - dataptr[DCTSIZE*7];
4129     tmp1 = dataptr[DCTSIZE*1] - dataptr[DCTSIZE*6];
4130     tmp2 = dataptr[DCTSIZE*2] - dataptr[DCTSIZE*5];
4131     tmp3 = dataptr[DCTSIZE*3] - dataptr[DCTSIZE*4];
4132
4133     dataptr[DCTSIZE*0] = (DCTELEM) RIGHT_SHIFT(tmp10 + tmp11, PASS1_BITS);
4134     dataptr[DCTSIZE*4] = (DCTELEM) RIGHT_SHIFT(tmp10 - tmp11, PASS1_BITS);
4135
4136     z1 = MULTIPLY(tmp12 + tmp13, FIX_0_541196100);       /* c6 */
4137     /* Add fudge factor here for final descale. */
4138     z1 += ONE << (CONST_BITS+PASS1_BITS-1);
4139
4140     dataptr[DCTSIZE*2] = (DCTELEM)
4141       RIGHT_SHIFT(z1 + MULTIPLY(tmp12, FIX_0_765366865), /* c2-c6 */
4142                   CONST_BITS+PASS1_BITS);
4143     dataptr[DCTSIZE*6] = (DCTELEM)
4144       RIGHT_SHIFT(z1 - MULTIPLY(tmp13, FIX_1_847759065), /* c2+c6 */
4145                   CONST_BITS+PASS1_BITS);
4146
4147     /* Odd part per figure 8 --- note paper omits factor of sqrt(2).
4148      * i0..i3 in the paper are tmp0..tmp3 here.
4149      */
4150
4151     tmp12 = tmp0 + tmp2;
4152     tmp13 = tmp1 + tmp3;
4153
4154     z1 = MULTIPLY(tmp12 + tmp13, FIX_1_175875602);       /*  c3 */
4155     /* Add fudge factor here for final descale. */
4156     z1 += ONE << (CONST_BITS+PASS1_BITS-1);
4157
4158     tmp12 = MULTIPLY(tmp12, - FIX_0_390180644);          /* -c3+c5 */
4159     tmp13 = MULTIPLY(tmp13, - FIX_1_961570560);          /* -c3-c5 */
4160     tmp12 += z1;
4161     tmp13 += z1;
4162
4163     z1 = MULTIPLY(tmp0 + tmp3, - FIX_0_899976223);       /* -c3+c7 */
4164     tmp0 = MULTIPLY(tmp0, FIX_1_501321110);              /*  c1+c3-c5-c7 */
4165     tmp3 = MULTIPLY(tmp3, FIX_0_298631336);              /* -c1+c3+c5-c7 */
4166     tmp0 += z1 + tmp12;
4167     tmp3 += z1 + tmp13;
4168
4169     z1 = MULTIPLY(tmp1 + tmp2, - FIX_2_562915447);       /* -c1-c3 */
4170     tmp1 = MULTIPLY(tmp1, FIX_3_072711026);              /*  c1+c3+c5-c7 */
4171     tmp2 = MULTIPLY(tmp2, FIX_2_053119869);              /*  c1+c3-c5+c7 */
4172     tmp1 += z1 + tmp13;
4173     tmp2 += z1 + tmp12;
4174
4175     dataptr[DCTSIZE*1] = (DCTELEM) RIGHT_SHIFT(tmp0, CONST_BITS+PASS1_BITS);
4176     dataptr[DCTSIZE*3] = (DCTELEM) RIGHT_SHIFT(tmp1, CONST_BITS+PASS1_BITS);
4177     dataptr[DCTSIZE*5] = (DCTELEM) RIGHT_SHIFT(tmp2, CONST_BITS+PASS1_BITS);
4178     dataptr[DCTSIZE*7] = (DCTELEM) RIGHT_SHIFT(tmp3, CONST_BITS+PASS1_BITS);
4179
4180     dataptr++;                  /* advance pointer to next column */
4181   }
4182 }
4183
4184
4185 /*
4186  * Perform the forward DCT on a 3x6 sample block.
4187  *
4188  * 3-point FDCT in pass 1 (rows), 6-point in pass 2 (columns).
4189  */
4190
4191 GLOBAL(void)
4192 jpeg_fdct_3x6 (DCTELEM * data, JSAMPARRAY sample_data, JDIMENSION start_col)
4193 {
4194   INT32 tmp0, tmp1, tmp2;
4195   INT32 tmp10, tmp11, tmp12;
4196   DCTELEM *dataptr;
4197   JSAMPROW elemptr;
4198   int ctr;
4199   SHIFT_TEMPS
4200
4201   /* Pre-zero output coefficient block. */
4202   MEMZERO(data, SIZEOF(DCTELEM) * DCTSIZE2);
4203
4204   /* Pass 1: process rows.
4205    * Note results are scaled up by sqrt(8) compared to a true DCT;
4206    * furthermore, we scale the results by 2**PASS1_BITS.
4207    * We scale the results further by 2 as part of output adaption
4208    * scaling for different DCT size.
4209    * 3-point FDCT kernel, cK represents sqrt(2) * cos(K*pi/6).
4210    */
4211
4212   dataptr = data;
4213   for (ctr = 0; ctr < 6; ctr++) {
4214     elemptr = sample_data[ctr] + start_col;
4215
4216     /* Even part */
4217
4218     tmp0 = GETJSAMPLE(elemptr[0]) + GETJSAMPLE(elemptr[2]);
4219     tmp1 = GETJSAMPLE(elemptr[1]);
4220
4221     tmp2 = GETJSAMPLE(elemptr[0]) - GETJSAMPLE(elemptr[2]);
4222
4223     /* Apply unsigned->signed conversion */
4224     dataptr[0] = (DCTELEM)
4225       ((tmp0 + tmp1 - 3 * CENTERJSAMPLE) << (PASS1_BITS+1));
4226     dataptr[2] = (DCTELEM)
4227       DESCALE(MULTIPLY(tmp0 - tmp1 - tmp1, FIX(0.707106781)), /* c2 */
4228               CONST_BITS-PASS1_BITS-1);
4229
4230     /* Odd part */
4231
4232     dataptr[1] = (DCTELEM)
4233       DESCALE(MULTIPLY(tmp2, FIX(1.224744871)),               /* c1 */
4234               CONST_BITS-PASS1_BITS-1);
4235
4236     dataptr += DCTSIZE;         /* advance pointer to next row */
4237   }
4238
4239   /* Pass 2: process columns.
4240    * We remove the PASS1_BITS scaling, but leave the results scaled up
4241    * by an overall factor of 8.
4242    * We must also scale the output by (8/6)*(8/3) = 32/9, which we partially
4243    * fold into the constant multipliers (other part was done in pass 1):
4244    * 6-point FDCT kernel, cK represents sqrt(2) * cos(K*pi/12) * 16/9.
4245    */
4246
4247   dataptr = data;
4248   for (ctr = 0; ctr < 3; ctr++) {
4249     /* Even part */
4250
4251     tmp0 = dataptr[DCTSIZE*0] + dataptr[DCTSIZE*5];
4252     tmp11 = dataptr[DCTSIZE*1] + dataptr[DCTSIZE*4];
4253     tmp2 = dataptr[DCTSIZE*2] + dataptr[DCTSIZE*3];
4254
4255     tmp10 = tmp0 + tmp2;
4256     tmp12 = tmp0 - tmp2;
4257
4258     tmp0 = dataptr[DCTSIZE*0] - dataptr[DCTSIZE*5];
4259     tmp1 = dataptr[DCTSIZE*1] - dataptr[DCTSIZE*4];
4260     tmp2 = dataptr[DCTSIZE*2] - dataptr[DCTSIZE*3];
4261
4262     dataptr[DCTSIZE*0] = (DCTELEM)
4263       DESCALE(MULTIPLY(tmp10 + tmp11, FIX(1.777777778)),         /* 16/9 */
4264               CONST_BITS+PASS1_BITS);
4265     dataptr[DCTSIZE*2] = (DCTELEM)
4266       DESCALE(MULTIPLY(tmp12, FIX(2.177324216)),                 /* c2 */
4267               CONST_BITS+PASS1_BITS);
4268     dataptr[DCTSIZE*4] = (DCTELEM)
4269       DESCALE(MULTIPLY(tmp10 - tmp11 - tmp11, FIX(1.257078722)), /* c4 */
4270               CONST_BITS+PASS1_BITS);
4271
4272     /* Odd part */
4273
4274     tmp10 = MULTIPLY(tmp0 + tmp2, FIX(0.650711829));             /* c5 */
4275
4276     dataptr[DCTSIZE*1] = (DCTELEM)
4277       DESCALE(tmp10 + MULTIPLY(tmp0 + tmp1, FIX(1.777777778)),   /* 16/9 */
4278               CONST_BITS+PASS1_BITS);
4279     dataptr[DCTSIZE*3] = (DCTELEM)
4280       DESCALE(MULTIPLY(tmp0 - tmp1 - tmp2, FIX(1.777777778)),    /* 16/9 */
4281               CONST_BITS+PASS1_BITS);
4282     dataptr[DCTSIZE*5] = (DCTELEM)
4283       DESCALE(tmp10 + MULTIPLY(tmp2 - tmp1, FIX(1.777777778)),   /* 16/9 */
4284               CONST_BITS+PASS1_BITS);
4285
4286     dataptr++;                  /* advance pointer to next column */
4287   }
4288 }
4289
4290
4291 /*
4292  * Perform the forward DCT on a 2x4 sample block.
4293  *
4294  * 2-point FDCT in pass 1 (rows), 4-point in pass 2 (columns).
4295  */
4296
4297 GLOBAL(void)
4298 jpeg_fdct_2x4 (DCTELEM * data, JSAMPARRAY sample_data, JDIMENSION start_col)
4299 {
4300   INT32 tmp0, tmp1;
4301   INT32 tmp10, tmp11;
4302   DCTELEM *dataptr;
4303   JSAMPROW elemptr;
4304   int ctr;
4305   SHIFT_TEMPS
4306
4307   /* Pre-zero output coefficient block. */
4308   MEMZERO(data, SIZEOF(DCTELEM) * DCTSIZE2);
4309
4310   /* Pass 1: process rows.
4311    * Note results are scaled up by sqrt(8) compared to a true DCT.
4312    * We must also scale the output by (8/2)*(8/4) = 2**3, which we add here.
4313    */
4314
4315   dataptr = data;
4316   for (ctr = 0; ctr < 4; ctr++) {
4317     elemptr = sample_data[ctr] + start_col;
4318
4319     /* Even part */
4320
4321     tmp0 = GETJSAMPLE(elemptr[0]);
4322     tmp1 = GETJSAMPLE(elemptr[1]);
4323
4324     /* Apply unsigned->signed conversion */
4325     dataptr[0] = (DCTELEM) ((tmp0 + tmp1 - 2 * CENTERJSAMPLE) << 3);
4326
4327     /* Odd part */
4328
4329     dataptr[1] = (DCTELEM) ((tmp0 - tmp1) << 3);
4330
4331     dataptr += DCTSIZE;         /* advance pointer to next row */
4332   }
4333
4334   /* Pass 2: process columns.
4335    * We leave the results scaled up by an overall factor of 8.
4336    * 4-point FDCT kernel,
4337    * cK represents sqrt(2) * cos(K*pi/16) [refers to 8-point FDCT].
4338    */
4339
4340   dataptr = data;
4341   for (ctr = 0; ctr < 2; ctr++) {
4342     /* Even part */
4343
4344     tmp0 = dataptr[DCTSIZE*0] + dataptr[DCTSIZE*3];
4345     tmp1 = dataptr[DCTSIZE*1] + dataptr[DCTSIZE*2];
4346
4347     tmp10 = dataptr[DCTSIZE*0] - dataptr[DCTSIZE*3];
4348     tmp11 = dataptr[DCTSIZE*1] - dataptr[DCTSIZE*2];
4349
4350     dataptr[DCTSIZE*0] = (DCTELEM) (tmp0 + tmp1);
4351     dataptr[DCTSIZE*2] = (DCTELEM) (tmp0 - tmp1);
4352
4353     /* Odd part */
4354
4355     tmp0 = MULTIPLY(tmp10 + tmp11, FIX_0_541196100);       /* c6 */
4356     /* Add fudge factor here for final descale. */
4357     tmp0 += ONE << (CONST_BITS-1);
4358
4359     dataptr[DCTSIZE*1] = (DCTELEM)
4360       RIGHT_SHIFT(tmp0 + MULTIPLY(tmp10, FIX_0_765366865), /* c2-c6 */
4361                   CONST_BITS);
4362     dataptr[DCTSIZE*3] = (DCTELEM)
4363       RIGHT_SHIFT(tmp0 - MULTIPLY(tmp11, FIX_1_847759065), /* c2+c6 */
4364                   CONST_BITS);
4365
4366     dataptr++;                  /* advance pointer to next column */
4367   }
4368 }
4369
4370
4371 /*
4372  * Perform the forward DCT on a 1x2 sample block.
4373  *
4374  * 1-point FDCT in pass 1 (rows), 2-point in pass 2 (columns).
4375  */
4376
4377 GLOBAL(void)
4378 jpeg_fdct_1x2 (DCTELEM * data, JSAMPARRAY sample_data, JDIMENSION start_col)
4379 {
4380   INT32 tmp0, tmp1;
4381
4382   /* Pre-zero output coefficient block. */
4383   MEMZERO(data, SIZEOF(DCTELEM) * DCTSIZE2);
4384
4385   /* Pass 1: empty. */
4386
4387   /* Pass 2: process columns.
4388    * We leave the results scaled up by an overall factor of 8.
4389    * We must also scale the output by (8/1)*(8/2) = 2**5.
4390    */
4391
4392   /* Even part */
4393
4394   tmp0 = GETJSAMPLE(sample_data[0][start_col]);
4395   tmp1 = GETJSAMPLE(sample_data[1][start_col]);
4396
4397   /* Apply unsigned->signed conversion */
4398   data[DCTSIZE*0] = (DCTELEM) ((tmp0 + tmp1 - 2 * CENTERJSAMPLE) << 5);
4399
4400   /* Odd part */
4401
4402   data[DCTSIZE*1] = (DCTELEM) ((tmp0 - tmp1) << 5);
4403 }
4404
4405 #endif /* DCT_SCALING_SUPPORTED */
4406 #endif /* DCT_ISLOW_SUPPORTED */