Move Alpha optimized IDCT to own file. Based on a patch by Måns
[ffmpeg.git] / libavcodec / simple_idct.c
1 /*
2  * Simple IDCT
3  *
4  * Copyright (c) 2001 Michael Niedermayer <michaelni@gmx.at>
5  *
6  * This library is free software; you can redistribute it and/or
7  * modify it under the terms of the GNU Lesser General Public
8  * License as published by the Free Software Foundation; either
9  * version 2 of the License, or (at your option) any later version.
10  *
11  * This library is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
14  * Lesser General Public License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public
17  * License along with this library; if not, write to the Free Software
18  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
19  */
20 /*
21   based upon some outcommented c code from mpeg2dec (idct_mmx.c
22   written by Aaron Holtzman <aholtzma@ess.engr.uvic.ca>) 
23  */
24 #include "avcodec.h"
25 #include "dsputil.h"
26 #include "simple_idct.h"
27
28 #if 0
29 #define W1 2841 /* 2048*sqrt (2)*cos (1*pi/16) */
30 #define W2 2676 /* 2048*sqrt (2)*cos (2*pi/16) */
31 #define W3 2408 /* 2048*sqrt (2)*cos (3*pi/16) */
32 #define W4 2048 /* 2048*sqrt (2)*cos (4*pi/16) */
33 #define W5 1609 /* 2048*sqrt (2)*cos (5*pi/16) */
34 #define W6 1108 /* 2048*sqrt (2)*cos (6*pi/16) */
35 #define W7 565  /* 2048*sqrt (2)*cos (7*pi/16) */
36 #define ROW_SHIFT 8
37 #define COL_SHIFT 17
38 #else
39 #define W1  22725  //cos(i*M_PI/16)*sqrt(2)*(1<<14) + 0.5
40 #define W2  21407  //cos(i*M_PI/16)*sqrt(2)*(1<<14) + 0.5
41 #define W3  19266  //cos(i*M_PI/16)*sqrt(2)*(1<<14) + 0.5
42 #define W4  16383  //cos(i*M_PI/16)*sqrt(2)*(1<<14) + 0.5
43 #define W5  12873  //cos(i*M_PI/16)*sqrt(2)*(1<<14) + 0.5
44 #define W6  8867   //cos(i*M_PI/16)*sqrt(2)*(1<<14) + 0.5
45 #define W7  4520   //cos(i*M_PI/16)*sqrt(2)*(1<<14) + 0.5
46 #define ROW_SHIFT 11
47 #define COL_SHIFT 20 // 6
48 #endif
49
50 #if defined(ARCH_POWERPC_405)
51
52 /* signed 16x16 -> 32 multiply add accumulate */
53 #define MAC16(rt, ra, rb) \
54     asm ("maclhw %0, %2, %3" : "=r" (rt) : "0" (rt), "r" (ra), "r" (rb));
55
56 /* signed 16x16 -> 32 multiply */
57 #define MUL16(rt, ra, rb) \
58     asm ("mullhw %0, %1, %2" : "=r" (rt) : "r" (ra), "r" (rb));
59
60 #else
61
62 /* signed 16x16 -> 32 multiply add accumulate */
63 #define MAC16(rt, ra, rb) rt += (ra) * (rb)
64
65 /* signed 16x16 -> 32 multiply */
66 #define MUL16(rt, ra, rb) rt = (ra) * (rb)
67
68 #endif
69
70 static inline void idctRowCondDC (int16_t * row)
71 {
72         int a0, a1, a2, a3, b0, b1, b2, b3;
73 #ifdef FAST_64BIT
74         uint64_t temp;
75 #else
76         uint32_t temp;
77 #endif
78
79 #ifdef FAST_64BIT
80 #ifdef WORDS_BIGENDIAN
81 #define ROW0_MASK 0xffff000000000000LL
82 #else
83 #define ROW0_MASK 0xffffLL
84 #endif
85         if ( ((((uint64_t *)row)[0] & ~ROW0_MASK) | 
86               ((uint64_t *)row)[1]) == 0) {
87             temp = (row[0] << 3) & 0xffff;
88             temp += temp << 16;
89             temp += temp << 32;
90             ((uint64_t *)row)[0] = temp;
91             ((uint64_t *)row)[1] = temp;
92             return;
93         }
94 #else
95         if (!(((uint32_t*)row)[1] |
96               ((uint32_t*)row)[2] |
97               ((uint32_t*)row)[3] | 
98               row[1])) {
99             temp = (row[0] << 3) & 0xffff;
100             temp += temp << 16;
101             ((uint32_t*)row)[0]=((uint32_t*)row)[1] =
102                 ((uint32_t*)row)[2]=((uint32_t*)row)[3] = temp;
103                 return;
104         }
105 #endif
106
107         a0 = (W4 * row[0]) + (1 << (ROW_SHIFT - 1));
108         a1 = a0;
109         a2 = a0;
110         a3 = a0;
111
112         /* no need to optimize : gcc does it */
113         a0 += W2 * row[2];
114         a1 += W6 * row[2];
115         a2 -= W6 * row[2];
116         a3 -= W2 * row[2];
117
118         MUL16(b0, W1, row[1]);
119         MAC16(b0, W3, row[3]);
120         MUL16(b1, W3, row[1]);
121         MAC16(b1, -W7, row[3]);
122         MUL16(b2, W5, row[1]);
123         MAC16(b2, -W1, row[3]);
124         MUL16(b3, W7, row[1]);
125         MAC16(b3, -W5, row[3]);
126
127 #ifdef FAST_64BIT
128         temp = ((uint64_t*)row)[1];
129 #else
130         temp = ((uint32_t*)row)[2] | ((uint32_t*)row)[3];
131 #endif
132         if (temp != 0) {
133             a0 += W4*row[4] + W6*row[6];
134             a1 += - W4*row[4] - W2*row[6];
135             a2 += - W4*row[4] + W2*row[6];
136             a3 += W4*row[4] - W6*row[6];
137
138             MAC16(b0, W5, row[5]);
139             MAC16(b0, W7, row[7]);
140             
141             MAC16(b1, -W1, row[5]);
142             MAC16(b1, -W5, row[7]);
143             
144             MAC16(b2, W7, row[5]);
145             MAC16(b2, W3, row[7]);
146             
147             MAC16(b3, W3, row[5]);
148             MAC16(b3, -W1, row[7]);
149         }
150
151         row[0] = (a0 + b0) >> ROW_SHIFT;
152         row[7] = (a0 - b0) >> ROW_SHIFT;
153         row[1] = (a1 + b1) >> ROW_SHIFT;
154         row[6] = (a1 - b1) >> ROW_SHIFT;
155         row[2] = (a2 + b2) >> ROW_SHIFT;
156         row[5] = (a2 - b2) >> ROW_SHIFT;
157         row[3] = (a3 + b3) >> ROW_SHIFT;
158         row[4] = (a3 - b3) >> ROW_SHIFT;
159 }
160
161 static inline void idctSparseColPut (UINT8 *dest, int line_size, 
162                                      int16_t * col)
163 {
164         int a0, a1, a2, a3, b0, b1, b2, b3;
165         UINT8 *cm = cropTbl + MAX_NEG_CROP;
166
167         /* XXX: I did that only to give same values as previous code */
168         a0 = W4 * (col[8*0] + ((1<<(COL_SHIFT-1))/W4));
169         a1 = a0;
170         a2 = a0;
171         a3 = a0;
172
173         a0 +=  + W2*col[8*2];
174         a1 +=  + W6*col[8*2];
175         a2 +=  - W6*col[8*2];
176         a3 +=  - W2*col[8*2];
177
178         MUL16(b0, W1, col[8*1]);
179         MUL16(b1, W3, col[8*1]);
180         MUL16(b2, W5, col[8*1]);
181         MUL16(b3, W7, col[8*1]);
182
183         MAC16(b0, + W3, col[8*3]);
184         MAC16(b1, - W7, col[8*3]);
185         MAC16(b2, - W1, col[8*3]);
186         MAC16(b3, - W5, col[8*3]);
187
188         if(col[8*4]){
189             a0 += + W4*col[8*4];
190             a1 += - W4*col[8*4];
191             a2 += - W4*col[8*4];
192             a3 += + W4*col[8*4];
193         }
194
195         if (col[8*5]) {
196             MAC16(b0, + W5, col[8*5]);
197             MAC16(b1, - W1, col[8*5]);
198             MAC16(b2, + W7, col[8*5]);
199             MAC16(b3, + W3, col[8*5]);
200         }
201
202         if(col[8*6]){
203             a0 += + W6*col[8*6];
204             a1 += - W2*col[8*6];
205             a2 += + W2*col[8*6];
206             a3 += - W6*col[8*6];
207         }
208
209         if (col[8*7]) {
210             MAC16(b0, + W7, col[8*7]);
211             MAC16(b1, - W5, col[8*7]);
212             MAC16(b2, + W3, col[8*7]);
213             MAC16(b3, - W1, col[8*7]);
214         }
215
216         dest[0] = cm[(a0 + b0) >> COL_SHIFT];
217         dest += line_size;
218         dest[0] = cm[(a1 + b1) >> COL_SHIFT];
219         dest += line_size;
220         dest[0] = cm[(a2 + b2) >> COL_SHIFT];
221         dest += line_size;
222         dest[0] = cm[(a3 + b3) >> COL_SHIFT];
223         dest += line_size;
224         dest[0] = cm[(a3 - b3) >> COL_SHIFT];
225         dest += line_size;
226         dest[0] = cm[(a2 - b2) >> COL_SHIFT];
227         dest += line_size;
228         dest[0] = cm[(a1 - b1) >> COL_SHIFT];
229         dest += line_size;
230         dest[0] = cm[(a0 - b0) >> COL_SHIFT];
231 }
232
233 static inline void idctSparseColAdd (UINT8 *dest, int line_size, 
234                                      int16_t * col)
235 {
236         int a0, a1, a2, a3, b0, b1, b2, b3;
237         UINT8 *cm = cropTbl + MAX_NEG_CROP;
238
239         /* XXX: I did that only to give same values as previous code */
240         a0 = W4 * (col[8*0] + ((1<<(COL_SHIFT-1))/W4));
241         a1 = a0;
242         a2 = a0;
243         a3 = a0;
244
245         a0 +=  + W2*col[8*2];
246         a1 +=  + W6*col[8*2];
247         a2 +=  - W6*col[8*2];
248         a3 +=  - W2*col[8*2];
249
250         MUL16(b0, W1, col[8*1]);
251         MUL16(b1, W3, col[8*1]);
252         MUL16(b2, W5, col[8*1]);
253         MUL16(b3, W7, col[8*1]);
254
255         MAC16(b0, + W3, col[8*3]);
256         MAC16(b1, - W7, col[8*3]);
257         MAC16(b2, - W1, col[8*3]);
258         MAC16(b3, - W5, col[8*3]);
259
260         if(col[8*4]){
261             a0 += + W4*col[8*4];
262             a1 += - W4*col[8*4];
263             a2 += - W4*col[8*4];
264             a3 += + W4*col[8*4];
265         }
266
267         if (col[8*5]) {
268             MAC16(b0, + W5, col[8*5]);
269             MAC16(b1, - W1, col[8*5]);
270             MAC16(b2, + W7, col[8*5]);
271             MAC16(b3, + W3, col[8*5]);
272         }
273
274         if(col[8*6]){
275             a0 += + W6*col[8*6];
276             a1 += - W2*col[8*6];
277             a2 += + W2*col[8*6];
278             a3 += - W6*col[8*6];
279         }
280
281         if (col[8*7]) {
282             MAC16(b0, + W7, col[8*7]);
283             MAC16(b1, - W5, col[8*7]);
284             MAC16(b2, + W3, col[8*7]);
285             MAC16(b3, - W1, col[8*7]);
286         }
287
288         dest[0] = cm[dest[0] + ((a0 + b0) >> COL_SHIFT)];
289         dest += line_size;
290         dest[0] = cm[dest[0] + ((a1 + b1) >> COL_SHIFT)];
291         dest += line_size;
292         dest[0] = cm[dest[0] + ((a2 + b2) >> COL_SHIFT)];
293         dest += line_size;
294         dest[0] = cm[dest[0] + ((a3 + b3) >> COL_SHIFT)];
295         dest += line_size;
296         dest[0] = cm[dest[0] + ((a3 - b3) >> COL_SHIFT)];
297         dest += line_size;
298         dest[0] = cm[dest[0] + ((a2 - b2) >> COL_SHIFT)];
299         dest += line_size;
300         dest[0] = cm[dest[0] + ((a1 - b1) >> COL_SHIFT)];
301         dest += line_size;
302         dest[0] = cm[dest[0] + ((a0 - b0) >> COL_SHIFT)];
303 }
304
305 static inline void idctSparseCol (int16_t * col)
306 {
307         int a0, a1, a2, a3, b0, b1, b2, b3;
308
309         /* XXX: I did that only to give same values as previous code */
310         a0 = W4 * (col[8*0] + ((1<<(COL_SHIFT-1))/W4));
311         a1 = a0;
312         a2 = a0;
313         a3 = a0;
314
315         a0 +=  + W2*col[8*2];
316         a1 +=  + W6*col[8*2];
317         a2 +=  - W6*col[8*2];
318         a3 +=  - W2*col[8*2];
319
320         MUL16(b0, W1, col[8*1]);
321         MUL16(b1, W3, col[8*1]);
322         MUL16(b2, W5, col[8*1]);
323         MUL16(b3, W7, col[8*1]);
324
325         MAC16(b0, + W3, col[8*3]);
326         MAC16(b1, - W7, col[8*3]);
327         MAC16(b2, - W1, col[8*3]);
328         MAC16(b3, - W5, col[8*3]);
329
330         if(col[8*4]){
331             a0 += + W4*col[8*4];
332             a1 += - W4*col[8*4];
333             a2 += - W4*col[8*4];
334             a3 += + W4*col[8*4];
335         }
336
337         if (col[8*5]) {
338             MAC16(b0, + W5, col[8*5]);
339             MAC16(b1, - W1, col[8*5]);
340             MAC16(b2, + W7, col[8*5]);
341             MAC16(b3, + W3, col[8*5]);
342         }
343
344         if(col[8*6]){
345             a0 += + W6*col[8*6];
346             a1 += - W2*col[8*6];
347             a2 += + W2*col[8*6];
348             a3 += - W6*col[8*6];
349         }
350
351         if (col[8*7]) {
352             MAC16(b0, + W7, col[8*7]);
353             MAC16(b1, - W5, col[8*7]);
354             MAC16(b2, + W3, col[8*7]);
355             MAC16(b3, - W1, col[8*7]);
356         }
357
358         col[0 ] = ((a0 + b0) >> COL_SHIFT);
359         col[8 ] = ((a1 + b1) >> COL_SHIFT);
360         col[16] = ((a2 + b2) >> COL_SHIFT);
361         col[24] = ((a3 + b3) >> COL_SHIFT);
362         col[32] = ((a3 - b3) >> COL_SHIFT);
363         col[40] = ((a2 - b2) >> COL_SHIFT);
364         col[48] = ((a1 - b1) >> COL_SHIFT);
365         col[56] = ((a0 - b0) >> COL_SHIFT);
366 }
367
368 void simple_idct_put(UINT8 *dest, int line_size, INT16 *block)
369 {
370     int i;
371     for(i=0; i<8; i++)
372         idctRowCondDC(block + i*8);
373     
374     for(i=0; i<8; i++)
375         idctSparseColPut(dest + i, line_size, block + i);
376 }
377
378 void simple_idct_add(UINT8 *dest, int line_size, INT16 *block)
379 {
380     int i;
381     for(i=0; i<8; i++)
382         idctRowCondDC(block + i*8);
383     
384     for(i=0; i<8; i++)
385         idctSparseColAdd(dest + i, line_size, block + i);
386 }
387
388 void simple_idct(INT16 *block)
389 {
390     int i;
391     for(i=0; i<8; i++)
392         idctRowCondDC(block + i*8);
393     
394     for(i=0; i<8; i++)
395         idctSparseCol(block + i);
396 }
397
398 /* 2x4x8 idct */
399
400 #define CN_SHIFT 12
401 #define C_FIX(x) ((int)((x) * (1 << CN_SHIFT) + 0.5))
402 #define C1 C_FIX(0.6532814824)
403 #define C2 C_FIX(0.2705980501)
404
405 /* row idct is multiple by 16 * sqrt(2.0), col idct4 is normalized,
406    and the butterfly must be multiplied by 0.5 * sqrt(2.0) */
407 #define C_SHIFT (4+1+12)
408
409 static inline void idct4col(UINT8 *dest, int line_size, const INT16 *col)
410 {
411     int c0, c1, c2, c3, a0, a1, a2, a3;
412     const UINT8 *cm = cropTbl + MAX_NEG_CROP;
413
414     a0 = col[8*0];
415     a1 = col[8*2];
416     a2 = col[8*4];
417     a3 = col[8*6];
418     c0 = ((a0 + a2) << (CN_SHIFT - 1)) + (1 << (C_SHIFT - 1));
419     c2 = ((a0 - a2) << (CN_SHIFT - 1)) + (1 << (C_SHIFT - 1));
420     c1 = a1 * C1 + a3 * C2;
421     c3 = a1 * C2 - a3 * C1;
422     dest[0] = cm[(c0 + c1) >> C_SHIFT];
423     dest += line_size;
424     dest[0] = cm[(c2 + c3) >> C_SHIFT];
425     dest += line_size;
426     dest[0] = cm[(c2 - c3) >> C_SHIFT];
427     dest += line_size;
428     dest[0] = cm[(c0 - c1) >> C_SHIFT];
429 }
430
431 #define BF(k) \
432 {\
433     int a0, a1;\
434     a0 = ptr[k];\
435     a1 = ptr[8 + k];\
436     ptr[k] = a0 + a1;\
437     ptr[8 + k] = a0 - a1;\
438 }
439
440 /* only used by DV codec. The input must be interlaced. 128 is added
441    to the pixels before clamping to avoid systematic error
442    (1024*sqrt(2)) offset would be needed otherwise. */
443 /* XXX: I think a 1.0/sqrt(2) normalization should be needed to
444    compensate the extra butterfly stage - I don't have the full DV
445    specification */
446 void simple_idct248_put(UINT8 *dest, int line_size, INT16 *block)
447 {
448     int i;
449     INT16 *ptr;
450     
451     /* butterfly */
452     ptr = block;
453     for(i=0;i<4;i++) {
454         BF(0);
455         BF(1);
456         BF(2);
457         BF(3);
458         BF(4);
459         BF(5);
460         BF(6);
461         BF(7);
462         ptr += 2 * 8;
463     }
464
465     /* IDCT8 on each line */
466     for(i=0; i<8; i++) {
467         idctRowCondDC(block + i*8);
468     }
469
470     /* IDCT4 and store */
471     for(i=0;i<8;i++) {
472         idct4col(dest + i, 2 * line_size, block + i);
473         idct4col(dest + line_size + i, 2 * line_size, block + 8 + i);
474     }
475 }