* UINTX -> uintx_t INTX -> intx_t
[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 (DCTELEM * 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(sizeof(DCTELEM)==2){
86             if ( ((((uint64_t *)row)[0] & ~ROW0_MASK) | 
87                   ((uint64_t *)row)[1]) == 0) {
88                 temp = (row[0] << 3) & 0xffff;
89                 temp += temp << 16;
90                 temp += temp << 32;
91                 ((uint64_t *)row)[0] = temp;
92                 ((uint64_t *)row)[1] = temp;
93                 return;
94             }
95         }else{
96             if (!(row[1]|row[2]|row[3]|row[4]|row[5]|row[6]|row[7])) {
97                 row[0]=row[1]=row[2]=row[3]=row[4]=row[5]=row[6]=row[7]= row[0] << 3;
98                 return;
99             }
100         }
101 #else
102         if(sizeof(DCTELEM)==2){
103             if (!(((uint32_t*)row)[1] |
104                   ((uint32_t*)row)[2] |
105                   ((uint32_t*)row)[3] | 
106                   row[1])) {
107                 temp = (row[0] << 3) & 0xffff;
108                 temp += temp << 16;
109                 ((uint32_t*)row)[0]=((uint32_t*)row)[1] =
110                 ((uint32_t*)row)[2]=((uint32_t*)row)[3] = temp;
111                 return;
112             }
113         }else{
114             if (!(row[1]|row[2]|row[3]|row[4]|row[5]|row[6]|row[7])) {
115                 row[0]=row[1]=row[2]=row[3]=row[4]=row[5]=row[6]=row[7]= row[0] << 3;
116                 return;
117             }
118         }
119 #endif
120
121         a0 = (W4 * row[0]) + (1 << (ROW_SHIFT - 1));
122         a1 = a0;
123         a2 = a0;
124         a3 = a0;
125
126         /* no need to optimize : gcc does it */
127         a0 += W2 * row[2];
128         a1 += W6 * row[2];
129         a2 -= W6 * row[2];
130         a3 -= W2 * row[2];
131
132         MUL16(b0, W1, row[1]);
133         MAC16(b0, W3, row[3]);
134         MUL16(b1, W3, row[1]);
135         MAC16(b1, -W7, row[3]);
136         MUL16(b2, W5, row[1]);
137         MAC16(b2, -W1, row[3]);
138         MUL16(b3, W7, row[1]);
139         MAC16(b3, -W5, row[3]);
140
141 #ifdef FAST_64BIT
142         temp = ((uint64_t*)row)[1];
143 #else
144         temp = ((uint32_t*)row)[2] | ((uint32_t*)row)[3];
145 #endif
146         if (temp != 0) {
147             a0 += W4*row[4] + W6*row[6];
148             a1 += - W4*row[4] - W2*row[6];
149             a2 += - W4*row[4] + W2*row[6];
150             a3 += W4*row[4] - W6*row[6];
151
152             MAC16(b0, W5, row[5]);
153             MAC16(b0, W7, row[7]);
154             
155             MAC16(b1, -W1, row[5]);
156             MAC16(b1, -W5, row[7]);
157             
158             MAC16(b2, W7, row[5]);
159             MAC16(b2, W3, row[7]);
160             
161             MAC16(b3, W3, row[5]);
162             MAC16(b3, -W1, row[7]);
163         }
164
165         row[0] = (a0 + b0) >> ROW_SHIFT;
166         row[7] = (a0 - b0) >> ROW_SHIFT;
167         row[1] = (a1 + b1) >> ROW_SHIFT;
168         row[6] = (a1 - b1) >> ROW_SHIFT;
169         row[2] = (a2 + b2) >> ROW_SHIFT;
170         row[5] = (a2 - b2) >> ROW_SHIFT;
171         row[3] = (a3 + b3) >> ROW_SHIFT;
172         row[4] = (a3 - b3) >> ROW_SHIFT;
173 }
174
175 static inline void idctSparseColPut (uint8_t *dest, int line_size, 
176                                      DCTELEM * col)
177 {
178         int a0, a1, a2, a3, b0, b1, b2, b3;
179         uint8_t *cm = cropTbl + MAX_NEG_CROP;
180
181         /* XXX: I did that only to give same values as previous code */
182         a0 = W4 * (col[8*0] + ((1<<(COL_SHIFT-1))/W4));
183         a1 = a0;
184         a2 = a0;
185         a3 = a0;
186
187         a0 +=  + W2*col[8*2];
188         a1 +=  + W6*col[8*2];
189         a2 +=  - W6*col[8*2];
190         a3 +=  - W2*col[8*2];
191
192         MUL16(b0, W1, col[8*1]);
193         MUL16(b1, W3, col[8*1]);
194         MUL16(b2, W5, col[8*1]);
195         MUL16(b3, W7, col[8*1]);
196
197         MAC16(b0, + W3, col[8*3]);
198         MAC16(b1, - W7, col[8*3]);
199         MAC16(b2, - W1, col[8*3]);
200         MAC16(b3, - W5, col[8*3]);
201
202         if(col[8*4]){
203             a0 += + W4*col[8*4];
204             a1 += - W4*col[8*4];
205             a2 += - W4*col[8*4];
206             a3 += + W4*col[8*4];
207         }
208
209         if (col[8*5]) {
210             MAC16(b0, + W5, col[8*5]);
211             MAC16(b1, - W1, col[8*5]);
212             MAC16(b2, + W7, col[8*5]);
213             MAC16(b3, + W3, col[8*5]);
214         }
215
216         if(col[8*6]){
217             a0 += + W6*col[8*6];
218             a1 += - W2*col[8*6];
219             a2 += + W2*col[8*6];
220             a3 += - W6*col[8*6];
221         }
222
223         if (col[8*7]) {
224             MAC16(b0, + W7, col[8*7]);
225             MAC16(b1, - W5, col[8*7]);
226             MAC16(b2, + W3, col[8*7]);
227             MAC16(b3, - W1, col[8*7]);
228         }
229
230         dest[0] = cm[(a0 + b0) >> COL_SHIFT];
231         dest += line_size;
232         dest[0] = cm[(a1 + b1) >> COL_SHIFT];
233         dest += line_size;
234         dest[0] = cm[(a2 + b2) >> COL_SHIFT];
235         dest += line_size;
236         dest[0] = cm[(a3 + b3) >> COL_SHIFT];
237         dest += line_size;
238         dest[0] = cm[(a3 - b3) >> COL_SHIFT];
239         dest += line_size;
240         dest[0] = cm[(a2 - b2) >> COL_SHIFT];
241         dest += line_size;
242         dest[0] = cm[(a1 - b1) >> COL_SHIFT];
243         dest += line_size;
244         dest[0] = cm[(a0 - b0) >> COL_SHIFT];
245 }
246
247 static inline void idctSparseColAdd (uint8_t *dest, int line_size, 
248                                      DCTELEM * col)
249 {
250         int a0, a1, a2, a3, b0, b1, b2, b3;
251         uint8_t *cm = cropTbl + MAX_NEG_CROP;
252
253         /* XXX: I did that only to give same values as previous code */
254         a0 = W4 * (col[8*0] + ((1<<(COL_SHIFT-1))/W4));
255         a1 = a0;
256         a2 = a0;
257         a3 = a0;
258
259         a0 +=  + W2*col[8*2];
260         a1 +=  + W6*col[8*2];
261         a2 +=  - W6*col[8*2];
262         a3 +=  - W2*col[8*2];
263
264         MUL16(b0, W1, col[8*1]);
265         MUL16(b1, W3, col[8*1]);
266         MUL16(b2, W5, col[8*1]);
267         MUL16(b3, W7, col[8*1]);
268
269         MAC16(b0, + W3, col[8*3]);
270         MAC16(b1, - W7, col[8*3]);
271         MAC16(b2, - W1, col[8*3]);
272         MAC16(b3, - W5, col[8*3]);
273
274         if(col[8*4]){
275             a0 += + W4*col[8*4];
276             a1 += - W4*col[8*4];
277             a2 += - W4*col[8*4];
278             a3 += + W4*col[8*4];
279         }
280
281         if (col[8*5]) {
282             MAC16(b0, + W5, col[8*5]);
283             MAC16(b1, - W1, col[8*5]);
284             MAC16(b2, + W7, col[8*5]);
285             MAC16(b3, + W3, col[8*5]);
286         }
287
288         if(col[8*6]){
289             a0 += + W6*col[8*6];
290             a1 += - W2*col[8*6];
291             a2 += + W2*col[8*6];
292             a3 += - W6*col[8*6];
293         }
294
295         if (col[8*7]) {
296             MAC16(b0, + W7, col[8*7]);
297             MAC16(b1, - W5, col[8*7]);
298             MAC16(b2, + W3, col[8*7]);
299             MAC16(b3, - W1, col[8*7]);
300         }
301
302         dest[0] = cm[dest[0] + ((a0 + b0) >> COL_SHIFT)];
303         dest += line_size;
304         dest[0] = cm[dest[0] + ((a1 + b1) >> COL_SHIFT)];
305         dest += line_size;
306         dest[0] = cm[dest[0] + ((a2 + b2) >> COL_SHIFT)];
307         dest += line_size;
308         dest[0] = cm[dest[0] + ((a3 + b3) >> COL_SHIFT)];
309         dest += line_size;
310         dest[0] = cm[dest[0] + ((a3 - b3) >> COL_SHIFT)];
311         dest += line_size;
312         dest[0] = cm[dest[0] + ((a2 - b2) >> COL_SHIFT)];
313         dest += line_size;
314         dest[0] = cm[dest[0] + ((a1 - b1) >> COL_SHIFT)];
315         dest += line_size;
316         dest[0] = cm[dest[0] + ((a0 - b0) >> COL_SHIFT)];
317 }
318
319 static inline void idctSparseCol (DCTELEM * col)
320 {
321         int a0, a1, a2, a3, b0, b1, b2, b3;
322
323         /* XXX: I did that only to give same values as previous code */
324         a0 = W4 * (col[8*0] + ((1<<(COL_SHIFT-1))/W4));
325         a1 = a0;
326         a2 = a0;
327         a3 = a0;
328
329         a0 +=  + W2*col[8*2];
330         a1 +=  + W6*col[8*2];
331         a2 +=  - W6*col[8*2];
332         a3 +=  - W2*col[8*2];
333
334         MUL16(b0, W1, col[8*1]);
335         MUL16(b1, W3, col[8*1]);
336         MUL16(b2, W5, col[8*1]);
337         MUL16(b3, W7, col[8*1]);
338
339         MAC16(b0, + W3, col[8*3]);
340         MAC16(b1, - W7, col[8*3]);
341         MAC16(b2, - W1, col[8*3]);
342         MAC16(b3, - W5, col[8*3]);
343
344         if(col[8*4]){
345             a0 += + W4*col[8*4];
346             a1 += - W4*col[8*4];
347             a2 += - W4*col[8*4];
348             a3 += + W4*col[8*4];
349         }
350
351         if (col[8*5]) {
352             MAC16(b0, + W5, col[8*5]);
353             MAC16(b1, - W1, col[8*5]);
354             MAC16(b2, + W7, col[8*5]);
355             MAC16(b3, + W3, col[8*5]);
356         }
357
358         if(col[8*6]){
359             a0 += + W6*col[8*6];
360             a1 += - W2*col[8*6];
361             a2 += + W2*col[8*6];
362             a3 += - W6*col[8*6];
363         }
364
365         if (col[8*7]) {
366             MAC16(b0, + W7, col[8*7]);
367             MAC16(b1, - W5, col[8*7]);
368             MAC16(b2, + W3, col[8*7]);
369             MAC16(b3, - W1, col[8*7]);
370         }
371
372         col[0 ] = ((a0 + b0) >> COL_SHIFT);
373         col[8 ] = ((a1 + b1) >> COL_SHIFT);
374         col[16] = ((a2 + b2) >> COL_SHIFT);
375         col[24] = ((a3 + b3) >> COL_SHIFT);
376         col[32] = ((a3 - b3) >> COL_SHIFT);
377         col[40] = ((a2 - b2) >> COL_SHIFT);
378         col[48] = ((a1 - b1) >> COL_SHIFT);
379         col[56] = ((a0 - b0) >> COL_SHIFT);
380 }
381
382 void simple_idct_put(uint8_t *dest, int line_size, DCTELEM *block)
383 {
384     int i;
385     for(i=0; i<8; i++)
386         idctRowCondDC(block + i*8);
387     
388     for(i=0; i<8; i++)
389         idctSparseColPut(dest + i, line_size, block + i);
390 }
391
392 void simple_idct_add(uint8_t *dest, int line_size, DCTELEM *block)
393 {
394     int i;
395     for(i=0; i<8; i++)
396         idctRowCondDC(block + i*8);
397     
398     for(i=0; i<8; i++)
399         idctSparseColAdd(dest + i, line_size, block + i);
400 }
401
402 void simple_idct(DCTELEM *block)
403 {
404     int i;
405     for(i=0; i<8; i++)
406         idctRowCondDC(block + i*8);
407     
408     for(i=0; i<8; i++)
409         idctSparseCol(block + i);
410 }
411
412 /* 2x4x8 idct */
413
414 #define CN_SHIFT 12
415 #define C_FIX(x) ((int)((x) * (1 << CN_SHIFT) + 0.5))
416 #define C1 C_FIX(0.6532814824)
417 #define C2 C_FIX(0.2705980501)
418
419 /* row idct is multiple by 16 * sqrt(2.0), col idct4 is normalized,
420    and the butterfly must be multiplied by 0.5 * sqrt(2.0) */
421 #define C_SHIFT (4+1+12)
422
423 static inline void idct4col(uint8_t *dest, int line_size, const DCTELEM *col)
424 {
425     int c0, c1, c2, c3, a0, a1, a2, a3;
426     const uint8_t *cm = cropTbl + MAX_NEG_CROP;
427
428     a0 = col[8*0];
429     a1 = col[8*2];
430     a2 = col[8*4];
431     a3 = col[8*6];
432     c0 = ((a0 + a2) << (CN_SHIFT - 1)) + (1 << (C_SHIFT - 1));
433     c2 = ((a0 - a2) << (CN_SHIFT - 1)) + (1 << (C_SHIFT - 1));
434     c1 = a1 * C1 + a3 * C2;
435     c3 = a1 * C2 - a3 * C1;
436     dest[0] = cm[(c0 + c1) >> C_SHIFT];
437     dest += line_size;
438     dest[0] = cm[(c2 + c3) >> C_SHIFT];
439     dest += line_size;
440     dest[0] = cm[(c2 - c3) >> C_SHIFT];
441     dest += line_size;
442     dest[0] = cm[(c0 - c1) >> C_SHIFT];
443 }
444
445 #define BF(k) \
446 {\
447     int a0, a1;\
448     a0 = ptr[k];\
449     a1 = ptr[8 + k];\
450     ptr[k] = a0 + a1;\
451     ptr[8 + k] = a0 - a1;\
452 }
453
454 /* only used by DV codec. The input must be interlaced. 128 is added
455    to the pixels before clamping to avoid systematic error
456    (1024*sqrt(2)) offset would be needed otherwise. */
457 /* XXX: I think a 1.0/sqrt(2) normalization should be needed to
458    compensate the extra butterfly stage - I don't have the full DV
459    specification */
460 void simple_idct248_put(uint8_t *dest, int line_size, DCTELEM *block)
461 {
462     int i;
463     DCTELEM *ptr;
464     
465     /* butterfly */
466     ptr = block;
467     for(i=0;i<4;i++) {
468         BF(0);
469         BF(1);
470         BF(2);
471         BF(3);
472         BF(4);
473         BF(5);
474         BF(6);
475         BF(7);
476         ptr += 2 * 8;
477     }
478
479     /* IDCT8 on each line */
480     for(i=0; i<8; i++) {
481         idctRowCondDC(block + i*8);
482     }
483
484     /* IDCT4 and store */
485     for(i=0;i<8;i++) {
486         idct4col(dest + i, 2 * line_size, block + i);
487         idct4col(dest + line_size + i, 2 * line_size, block + 8 + i);
488     }
489 }
490
491 /* 8x4 & 4x8 WMV2 IDCT */
492 #undef CN_SHIFT
493 #undef C_SHIFT
494 #undef C_FIX
495 #undef C1
496 #undef C2
497 #define CN_SHIFT 12
498 #define C_FIX(x) ((int)((x) * 1.414213562 * (1 << CN_SHIFT) + 0.5))
499 #define C1 C_FIX(0.6532814824)
500 #define C2 C_FIX(0.2705980501)
501 #define C3 C_FIX(0.5)
502 #define C_SHIFT (4+1+12)
503 static inline void idct4col_add(uint8_t *dest, int line_size, const DCTELEM *col)
504 {
505     int c0, c1, c2, c3, a0, a1, a2, a3;
506     const uint8_t *cm = cropTbl + MAX_NEG_CROP;
507
508     a0 = col[8*0];
509     a1 = col[8*1];
510     a2 = col[8*2];
511     a3 = col[8*3];
512     c0 = (a0 + a2)*C3 + (1 << (C_SHIFT - 1));
513     c2 = (a0 - a2)*C3 + (1 << (C_SHIFT - 1));
514     c1 = a1 * C1 + a3 * C2;
515     c3 = a1 * C2 - a3 * C1;
516     dest[0] = cm[dest[0] + ((c0 + c1) >> C_SHIFT)];
517     dest += line_size;
518     dest[0] = cm[dest[0] + ((c2 + c3) >> C_SHIFT)];
519     dest += line_size;
520     dest[0] = cm[dest[0] + ((c2 - c3) >> C_SHIFT)];
521     dest += line_size;
522     dest[0] = cm[dest[0] + ((c0 - c1) >> C_SHIFT)];
523 }
524
525 #define RN_SHIFT 15
526 #define R_FIX(x) ((int)((x) * 1.414213562 * (1 << RN_SHIFT) + 0.5))
527 #define R1 R_FIX(0.6532814824)
528 #define R2 R_FIX(0.2705980501)
529 #define R3 R_FIX(0.5)
530 #define R_SHIFT 11
531 static inline void idct4row(DCTELEM *row)
532 {
533     int c0, c1, c2, c3, a0, a1, a2, a3;
534     //const uint8_t *cm = cropTbl + MAX_NEG_CROP;
535
536     a0 = row[0];
537     a1 = row[1];
538     a2 = row[2];
539     a3 = row[3];
540     c0 = (a0 + a2)*R3 + (1 << (R_SHIFT - 1));
541     c2 = (a0 - a2)*R3 + (1 << (R_SHIFT - 1));
542     c1 = a1 * R1 + a3 * R2;
543     c3 = a1 * R2 - a3 * R1;
544     row[0]= (c0 + c1) >> R_SHIFT;
545     row[1]= (c2 + c3) >> R_SHIFT;
546     row[2]= (c2 - c3) >> R_SHIFT;
547     row[3]= (c0 - c1) >> R_SHIFT;
548 }
549
550 void simple_idct84_add(uint8_t *dest, int line_size, DCTELEM *block)
551 {
552     int i;
553
554     /* IDCT8 on each line */
555     for(i=0; i<4; i++) {
556         idctRowCondDC(block + i*8);
557     }
558
559     /* IDCT4 and store */
560     for(i=0;i<8;i++) {
561         idct4col_add(dest + i, line_size, block + i);
562     }
563 }
564
565 void simple_idct48_add(uint8_t *dest, int line_size, DCTELEM *block)
566 {
567     int i;
568
569     /* IDCT4 on each line */
570     for(i=0; i<8; i++) {
571         idct4row(block + i*8);
572     }
573
574     /* IDCT8 and store */
575     for(i=0; i<4; i++){
576         idctSparseColAdd(dest + i, line_size, block + i);
577     }
578 }
579