1 
/* 
2 
* jfdctfst.c 
3 
* 
4 
* Copyright (C) 19941996, Thomas G. Lane. 
5 
* Modified 20032009 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 fast, not so accurate integer implementation of the 
10 
* forward DCT (Discrete Cosine Transform). 
11 
* 
12 
* A 2D DCT can be done by 1D DCT on each row followed by 1D 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 Arai, Agui, and Nakajima's algorithm for 
17 
* scaled DCT. Their original paper (Trans. IEICE E71(11):1095) is in 
18 
* Japanese, but the algorithm is described in the Pennebaker & Mitchell 
19 
* JPEG textbook (see REFERENCES section in file README). The following code 
20 
* is based directly on figure 48 in P&M. 
21 
* While an 8point DCT cannot be done in less than 11 multiplies, it is 
22 
* possible to arrange the computation so that many of the multiplies are 
23 
* simple scalings of the final outputs. These multiplies can then be 
24 
* folded into the multiplications or divisions by the JPEG quantization 
25 
* table entries. The AA&N method leaves only 5 multiplies and 29 adds 
26 
* to be done in the DCT itself. 
27 
* The primary disadvantage of this method is that with fixedpoint math, 
28 
* accuracy is lost due to imprecise representation of the scaled 
29 
* quantization values. The smaller the quantization table entry, the less 
30 
* precise the scaled value, so this implementation does worse with high 
31 
* qualitysetting files than with lowquality ones. 
32 
*/ 
33 

34 
#define JPEG_INTERNALS 
35 
#include "jinclude.h" 
36 
#include "jpeglib.h" 
37 
#include "jdct.h" /* Private declarations for DCT subsystem */ 
38 

39 
#ifdef DCT_IFAST_SUPPORTED 
40 

41 

42 
/* 
43 
* This module is specialized to the case DCTSIZE = 8. 
44 
*/ 
45 

46 
#if DCTSIZE != 8 
47 
Sorry, this code only copes with 8x8 DCTs. /* deliberate syntax err */ 
48 
#endif 
49 

50 

51 
/* Scaling decisions are generally the same as in the LL&M algorithm; 
52 
* see jfdctint.c for more details. However, we choose to descale 
53 
* (right shift) multiplication products as soon as they are formed, 
54 
* rather than carrying additional fractional bits into subsequent additions. 
55 
* This compromises accuracy slightly, but it lets us save a few shifts. 
56 
* More importantly, 16bit arithmetic is then adequate (for 8bit samples) 
57 
* everywhere except in the multiplications proper; this saves a good deal 
58 
* of work on 16bitint machines. 
59 
* 
60 
* Again to save a few shifts, the intermediate results between pass 1 and 
61 
* pass 2 are not upscaled, but are represented only to integral precision. 
62 
* 
63 
* A final compromise is to represent the multiplicative constants to only 
64 
* 8 fractional bits, rather than 13. This saves some shifting work on some 
65 
* machines, and may also reduce the cost of multiplication (since there 
66 
* are fewer onebits in the constants). 
67 
*/ 
68 

69 
#define CONST_BITS 8 
70 

71 

72 
/* Some C compilers fail to reduce "FIX(constant)" at compile time, thus 
73 
* causing a lot of useless floatingpoint operations at run time. 
74 
* To get around this we use the following precalculated constants. 
75 
* If you change CONST_BITS you may want to add appropriate values. 
76 
* (With a reasonable C compiler, you can just rely on the FIX() macro...) 
77 
*/ 
78 

79 
#if CONST_BITS == 8 
80 
#define FIX_0_382683433 ((INT32) 98) /* FIX(0.382683433) */ 
81 
#define FIX_0_541196100 ((INT32) 139) /* FIX(0.541196100) */ 
82 
#define FIX_0_707106781 ((INT32) 181) /* FIX(0.707106781) */ 
83 
#define FIX_1_306562965 ((INT32) 334) /* FIX(1.306562965) */ 
84 
#else 
85 
#define FIX_0_382683433 FIX(0.382683433) 
86 
#define FIX_0_541196100 FIX(0.541196100) 
87 
#define FIX_0_707106781 FIX(0.707106781) 
88 
#define FIX_1_306562965 FIX(1.306562965) 
89 
#endif 
90 

91 

92 
/* We can gain a little more speed, with a further compromise in accuracy, 
93 
* by omitting the addition in a descaling shift. This yields an incorrectly 
94 
* rounded result half the time... 
95 
*/ 
96 

97 
#ifndef USE_ACCURATE_ROUNDING 
98 
#undef DESCALE 
99 
#define DESCALE(x,n) RIGHT_SHIFT(x, n) 
100 
#endif 
101 

102 

103 
/* Multiply a DCTELEM variable by an INT32 constant, and immediately 
104 
* descale to yield a DCTELEM result. 
105 
*/ 
106 

107 
#define MULTIPLY(var,const) ((DCTELEM) DESCALE((var) * (const), CONST_BITS)) 
108 

109 

110 
/* 
111 
* Perform the forward DCT on one block of samples. 
112 
*/ 
113 

114 
GLOBAL(void) 
115 
jpeg_fdct_ifast (DCTELEM * data, JSAMPARRAY sample_data, JDIMENSION start_col) 
116 
{ 
117 
DCTELEM tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7; 
118 
DCTELEM tmp10, tmp11, tmp12, tmp13; 
119 
DCTELEM z1, z2, z3, z4, z5, z11, z13; 
120 
DCTELEM *dataptr; 
121 
JSAMPROW elemptr; 
122 
int ctr; 
123 
SHIFT_TEMPS 
124 

125 
/* Pass 1: process rows. */ 
126 

127 
dataptr = data; 
128 
for (ctr = 0; ctr < DCTSIZE; ctr++) { 
129 
elemptr = sample_data[ctr] + start_col; 
130 

131 
/* Load data into workspace */ 
132 
tmp0 = GETJSAMPLE(elemptr[0]) + GETJSAMPLE(elemptr[7]); 
133 
tmp7 = GETJSAMPLE(elemptr[0])  GETJSAMPLE(elemptr[7]); 
134 
tmp1 = GETJSAMPLE(elemptr[1]) + GETJSAMPLE(elemptr[6]); 
135 
tmp6 = GETJSAMPLE(elemptr[1])  GETJSAMPLE(elemptr[6]); 
136 
tmp2 = GETJSAMPLE(elemptr[2]) + GETJSAMPLE(elemptr[5]); 
137 
tmp5 = GETJSAMPLE(elemptr[2])  GETJSAMPLE(elemptr[5]); 
138 
tmp3 = GETJSAMPLE(elemptr[3]) + GETJSAMPLE(elemptr[4]); 
139 
tmp4 = GETJSAMPLE(elemptr[3])  GETJSAMPLE(elemptr[4]); 
140 

141 
/* Even part */ 
142 

143 
tmp10 = tmp0 + tmp3; /* phase 2 */ 
144 
tmp13 = tmp0  tmp3; 
145 
tmp11 = tmp1 + tmp2; 
146 
tmp12 = tmp1  tmp2; 
147 

148 
/* Apply unsigned>signed conversion */ 
149 
dataptr[0] = tmp10 + tmp11  8 * CENTERJSAMPLE; /* phase 3 */ 
150 
dataptr[4] = tmp10  tmp11; 
151 

152 
z1 = MULTIPLY(tmp12 + tmp13, FIX_0_707106781); /* c4 */ 
153 
dataptr[2] = tmp13 + z1; /* phase 5 */ 
154 
dataptr[6] = tmp13  z1; 
155 

156 
/* Odd part */ 
157 

158 
tmp10 = tmp4 + tmp5; /* phase 2 */ 
159 
tmp11 = tmp5 + tmp6; 
160 
tmp12 = tmp6 + tmp7; 
161 

162 
/* The rotator is modified from fig 48 to avoid extra negations. */ 
163 
z5 = MULTIPLY(tmp10  tmp12, FIX_0_382683433); /* c6 */ 
164 
z2 = MULTIPLY(tmp10, FIX_0_541196100) + z5; /* c2c6 */ 
165 
z4 = MULTIPLY(tmp12, FIX_1_306562965) + z5; /* c2+c6 */ 
166 
z3 = MULTIPLY(tmp11, FIX_0_707106781); /* c4 */ 
167 

168 
z11 = tmp7 + z3; /* phase 5 */ 
169 
z13 = tmp7  z3; 
170 

171 
dataptr[5] = z13 + z2; /* phase 6 */ 
172 
dataptr[3] = z13  z2; 
173 
dataptr[1] = z11 + z4; 
174 
dataptr[7] = z11  z4; 
175 

176 
dataptr += DCTSIZE; /* advance pointer to next row */ 
177 
} 
178 

179 
/* Pass 2: process columns. */ 
180 

181 
dataptr = data; 
182 
for (ctr = DCTSIZE1; ctr >= 0; ctr) { 
183 
tmp0 = dataptr[DCTSIZE*0] + dataptr[DCTSIZE*7]; 
184 
tmp7 = dataptr[DCTSIZE*0]  dataptr[DCTSIZE*7]; 
185 
tmp1 = dataptr[DCTSIZE*1] + dataptr[DCTSIZE*6]; 
186 
tmp6 = dataptr[DCTSIZE*1]  dataptr[DCTSIZE*6]; 
187 
tmp2 = dataptr[DCTSIZE*2] + dataptr[DCTSIZE*5]; 
188 
tmp5 = dataptr[DCTSIZE*2]  dataptr[DCTSIZE*5]; 
189 
tmp3 = dataptr[DCTSIZE*3] + dataptr[DCTSIZE*4]; 
190 
tmp4 = dataptr[DCTSIZE*3]  dataptr[DCTSIZE*4]; 
191 

192 
/* Even part */ 
193 

194 
tmp10 = tmp0 + tmp3; /* phase 2 */ 
195 
tmp13 = tmp0  tmp3; 
196 
tmp11 = tmp1 + tmp2; 
197 
tmp12 = tmp1  tmp2; 
198 

199 
dataptr[DCTSIZE*0] = tmp10 + tmp11; /* phase 3 */ 
200 
dataptr[DCTSIZE*4] = tmp10  tmp11; 
201 

202 
z1 = MULTIPLY(tmp12 + tmp13, FIX_0_707106781); /* c4 */ 
203 
dataptr[DCTSIZE*2] = tmp13 + z1; /* phase 5 */ 
204 
dataptr[DCTSIZE*6] = tmp13  z1; 
205 

206 
/* Odd part */ 
207 

208 
tmp10 = tmp4 + tmp5; /* phase 2 */ 
209 
tmp11 = tmp5 + tmp6; 
210 
tmp12 = tmp6 + tmp7; 
211 

212 
/* The rotator is modified from fig 48 to avoid extra negations. */ 
213 
z5 = MULTIPLY(tmp10  tmp12, FIX_0_382683433); /* c6 */ 
214 
z2 = MULTIPLY(tmp10, FIX_0_541196100) + z5; /* c2c6 */ 
215 
z4 = MULTIPLY(tmp12, FIX_1_306562965) + z5; /* c2+c6 */ 
216 
z3 = MULTIPLY(tmp11, FIX_0_707106781); /* c4 */ 
217 

218 
z11 = tmp7 + z3; /* phase 5 */ 
219 
z13 = tmp7  z3; 
220 

221 
dataptr[DCTSIZE*5] = z13 + z2; /* phase 6 */ 
222 
dataptr[DCTSIZE*3] = z13  z2; 
223 
dataptr[DCTSIZE*1] = z11 + z4; 
224 
dataptr[DCTSIZE*7] = z11  z4; 
225 

226 
dataptr++; /* advance pointer to next column */ 
227 
} 
228 
} 
229 

230 
#endif /* DCT_IFAST_SUPPORTED */ 