platform-independent restructuring and code simplification
[ffmpeg.git] / libavcodec / dct-test.c
1 /*
2  * (c) 2001 Fabrice Bellard
3  *     2007 Marc Hoffman <marc.hoffman@analog.com>
4  *
5  * This file is part of FFmpeg.
6  *
7  * FFmpeg is free software; you can redistribute it and/or
8  * modify it under the terms of the GNU Lesser General Public
9  * License as published by the Free Software Foundation; either
10  * version 2.1 of the License, or (at your option) any later version.
11  *
12  * FFmpeg is distributed in the hope that it will be useful,
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
15  * Lesser General Public License for more details.
16  *
17  * You should have received a copy of the GNU Lesser General Public
18  * License along with FFmpeg; if not, write to the Free Software
19  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
20  */
21
22 /**
23  * @file dct-test.c
24  * DCT test. (c) 2001 Fabrice Bellard.
25  * Started from sample code by Juan J. Sierralta P.
26  */
27
28 #include <stdlib.h>
29 #include <stdio.h>
30 #include <string.h>
31 #include <sys/time.h>
32 #include <unistd.h>
33
34 #include "dsputil.h"
35
36 #include "simple_idct.h"
37 #include "faandct.h"
38
39 #ifndef MAX
40 #define MAX(a, b)  (((a) > (b)) ? (a) : (b))
41 #endif
42
43 #undef printf
44
45 void *fast_memcpy(void *a, const void *b, size_t c){return memcpy(a,b,c);};
46
47 /* reference fdct/idct */
48 extern void fdct(DCTELEM *block);
49 extern void idct(DCTELEM *block);
50 extern void ff_idct_xvid_mmx(DCTELEM *block);
51 extern void ff_idct_xvid_mmx2(DCTELEM *block);
52 extern void init_fdct();
53
54 extern void j_rev_dct(DCTELEM *data);
55 extern void ff_mmx_idct(DCTELEM *data);
56 extern void ff_mmxext_idct(DCTELEM *data);
57
58 extern void odivx_idct_c (short *block);
59
60 // BFIN
61 extern void ff_bfin_idct (DCTELEM *block)  ;
62 extern void ff_bfin_fdct (DCTELEM *block) ;
63
64 // ALTIVEC
65 extern void fdct_altivec (DCTELEM *block);
66 //extern void idct_altivec (DCTELEM *block);?? no routine
67
68
69 struct algo {
70   char *name;
71   enum { FDCT, IDCT } is_idct;
72   void (* func) (DCTELEM *block);
73   void (* ref)  (DCTELEM *block);
74   enum formattag { NO_PERM,MMX_PERM, MMX_SIMPLE_PERM, SCALE_PERM } format;
75 };
76
77 #ifndef FAAN_POSTSCALE
78 #define FAAN_SCALE SCALE_PERM
79 #else
80 #define FAAN_SCALE NO_PERM
81 #endif
82
83 #define DCT_ERROR(name,is_idct,func,ref,form) {name,is_idct,func,ref,form}
84
85
86 struct algo algos[] = {
87   DCT_ERROR( "REF-DBL",        0, fdct,               fdct, NO_PERM),
88   DCT_ERROR("IJG-AAN-INT",     0, fdct_ifast,         fdct, SCALE_PERM),
89   DCT_ERROR("IJG-LLM-INT",     0, ff_jpeg_fdct_islow, fdct, NO_PERM),
90   DCT_ERROR("REF-DBL",         1, idct,               idct, NO_PERM),
91   DCT_ERROR("INT",             1, j_rev_dct,          idct, MMX_PERM),
92   DCT_ERROR("SIMPLE-C",        1, simple_idct,        idct, NO_PERM),
93
94 #ifdef ARCH_X86
95   DCT_ERROR("MMX",             0, ff_fdct_mmx,        fdct, NO_PERM),
96   DCT_ERROR("MMX2",            0, ff_fdct_mmx2,       fdct, NO_PERM),
97   DCT_ERROR("FAAN",            0, ff_faandct,         fdct, FAAN_SCALE),
98
99   DCT_ERROR("LIBMPEG2-MMX",    1, ff_mmx_idct,        idct, MMX_PERM),
100   DCT_ERROR("LIBMPEG2-MMXEXT", 1, ff_mmxext_idct,     idct, MMX_PERM),
101   DCT_ERROR("SIMPLE-MMX",      1, ff_simple_idct_mmx, idct, MMX_SIMPLE_PERM),
102   DCT_ERROR("XVID-MMX",        1, ff_idct_xvid_mmx,   idct, NO_PERM),
103   DCT_ERROR("XVID-MMX2",       1, ff_idct_xvid_mmx2,  idct, NO_PERM),
104 #endif
105
106 #ifdef HAVE_ALTIVEC
107   DCT_ERROR("altivecfdct",     0, fdct_altivec,       fdct, NO_PERM),
108 #endif
109
110 #ifdef ARCH_BFIN
111   DCT_ERROR("BFINfdct",        0, ff_bfin_fdct,       fdct, NO_PERM),
112   DCT_ERROR("BFINidct",        1, ff_bfin_idct,       idct, NO_PERM),
113 #endif
114
115   { 0 }
116 };
117
118 #define AANSCALE_BITS 12
119 static const unsigned short aanscales[64] = {
120     /* precomputed values scaled up by 14 bits */
121     16384, 22725, 21407, 19266, 16384, 12873,  8867,  4520,
122     22725, 31521, 29692, 26722, 22725, 17855, 12299,  6270,
123     21407, 29692, 27969, 25172, 21407, 16819, 11585,  5906,
124     19266, 26722, 25172, 22654, 19266, 15137, 10426,  5315,
125     16384, 22725, 21407, 19266, 16384, 12873,  8867,  4520,
126     12873, 17855, 16819, 15137, 12873, 10114,  6967,  3552,
127     8867, 12299, 11585, 10426,  8867,  6967,  4799,  2446,
128     4520,  6270,  5906,  5315,  4520,  3552,  2446,  1247
129 };
130
131 uint8_t cropTbl[256 + 2 * MAX_NEG_CROP];
132
133 int64_t gettime(void)
134 {
135     struct timeval tv;
136     gettimeofday(&tv,NULL);
137     return (int64_t)tv.tv_sec * 1000000 + tv.tv_usec;
138 }
139
140 #define NB_ITS 20000
141 #define NB_ITS_SPEED 50000
142
143 static short idct_mmx_perm[64];
144
145 static short idct_simple_mmx_perm[64]={
146         0x00, 0x08, 0x04, 0x09, 0x01, 0x0C, 0x05, 0x0D,
147         0x10, 0x18, 0x14, 0x19, 0x11, 0x1C, 0x15, 0x1D,
148         0x20, 0x28, 0x24, 0x29, 0x21, 0x2C, 0x25, 0x2D,
149         0x12, 0x1A, 0x16, 0x1B, 0x13, 0x1E, 0x17, 0x1F,
150         0x02, 0x0A, 0x06, 0x0B, 0x03, 0x0E, 0x07, 0x0F,
151         0x30, 0x38, 0x34, 0x39, 0x31, 0x3C, 0x35, 0x3D,
152         0x22, 0x2A, 0x26, 0x2B, 0x23, 0x2E, 0x27, 0x2F,
153         0x32, 0x3A, 0x36, 0x3B, 0x33, 0x3E, 0x37, 0x3F,
154 };
155
156 void idct_mmx_init(void)
157 {
158     int i;
159
160     /* the mmx/mmxext idct uses a reordered input, so we patch scan tables */
161     for (i = 0; i < 64; i++) {
162         idct_mmx_perm[i] = (i & 0x38) | ((i & 6) >> 1) | ((i & 1) << 2);
163 //        idct_simple_mmx_perm[i] = simple_block_permute_op(i);
164     }
165 }
166
167 static DCTELEM block[64] __attribute__ ((aligned (8)));
168 static DCTELEM block1[64] __attribute__ ((aligned (8)));
169 static DCTELEM block_org[64] __attribute__ ((aligned (8)));
170
171 void dct_error(const char *name, int is_idct,
172                void (*fdct_func)(DCTELEM *block),
173                void (*fdct_ref)(DCTELEM *block), int form, int test)
174 {
175     int it, i, scale;
176     int err_inf, v;
177     int64_t err2, ti, ti1, it1;
178     int64_t sysErr[64], sysErrMax=0;
179     int maxout=0;
180     int blockSumErrMax=0, blockSumErr;
181
182     srandom(0);
183
184     err_inf = 0;
185     err2 = 0;
186     for(i=0; i<64; i++) sysErr[i]=0;
187     for(it=0;it<NB_ITS;it++) {
188         for(i=0;i<64;i++)
189             block1[i] = 0;
190         switch(test){
191         case 0:
192             for(i=0;i<64;i++)
193                 block1[i] = (random() % 512) -256;
194             if (is_idct){
195                 fdct(block1);
196
197                 for(i=0;i<64;i++)
198                     block1[i]>>=3;
199             }
200         break;
201         case 1:{
202             int num= (random()%10)+1;
203             for(i=0;i<num;i++)
204                 block1[random()%64] = (random() % 512) -256;
205         }break;
206         case 2:
207             block1[0]= (random()%4096)-2048;
208             block1[63]= (block1[0]&1)^1;
209         break;
210         }
211
212 #if 0 // simulate mismatch control
213 { int sum=0;
214         for(i=0;i<64;i++)
215            sum+=block1[i];
216
217         if((sum&1)==0) block1[63]^=1;
218 }
219 #endif
220
221         for(i=0; i<64; i++)
222             block_org[i]= block1[i];
223
224         if (form == MMX_PERM) {
225             for(i=0;i<64;i++)
226                 block[idct_mmx_perm[i]] = block1[i];
227             } else if (form == MMX_SIMPLE_PERM) {
228             for(i=0;i<64;i++)
229                 block[idct_simple_mmx_perm[i]] = block1[i];
230
231         } else {
232             for(i=0; i<64; i++)
233                 block[i]= block1[i];
234         }
235 #if 0 // simulate mismatch control for tested IDCT but not the ref
236 { int sum=0;
237         for(i=0;i<64;i++)
238            sum+=block[i];
239
240         if((sum&1)==0) block[63]^=1;
241 }
242 #endif
243
244         fdct_func(block);
245         emms_c(); /* for ff_mmx_idct */
246
247         if (form == SCALE_PERM) {
248             for(i=0; i<64; i++) {
249                 scale = 8*(1 << (AANSCALE_BITS + 11)) / aanscales[i];
250                 block[i] = (block[i] * scale /*+ (1<<(AANSCALE_BITS-1))*/) >> AANSCALE_BITS;
251             }
252         }
253
254         fdct_ref(block1);
255
256         blockSumErr=0;
257         for(i=0;i<64;i++) {
258             v = abs(block[i] - block1[i]);
259             if (v > err_inf)
260                 err_inf = v;
261             err2 += v * v;
262             sysErr[i] += block[i] - block1[i];
263             blockSumErr += v;
264             if( abs(block[i])>maxout) maxout=abs(block[i]);
265         }
266         if(blockSumErrMax < blockSumErr) blockSumErrMax= blockSumErr;
267 #if 0 // print different matrix pairs
268         if(blockSumErr){
269             printf("\n");
270             for(i=0; i<64; i++){
271                 if((i&7)==0) printf("\n");
272                 printf("%4d ", block_org[i]);
273             }
274             for(i=0; i<64; i++){
275                 if((i&7)==0) printf("\n");
276                 printf("%4d ", block[i] - block1[i]);
277             }
278         }
279 #endif
280     }
281     for(i=0; i<64; i++) sysErrMax= MAX(sysErrMax, FFABS(sysErr[i]));
282
283 #if 1 // dump systematic errors
284     for(i=0; i<64; i++){
285         if(i%8==0) printf("\n");
286         printf("%5d ", (int)sysErr[i]);
287     }
288     printf("\n");
289 #endif
290
291     printf("%s %s: err_inf=%d err2=%0.8f syserr=%0.8f maxout=%d blockSumErr=%d\n",
292            is_idct ? "IDCT" : "DCT",
293            name, err_inf, (double)err2 / NB_ITS / 64.0, (double)sysErrMax / NB_ITS, maxout, blockSumErrMax);
294 #if 1 //Speed test
295     /* speed test */
296     for(i=0;i<64;i++)
297         block1[i] = 0;
298     switch(test){
299     case 0:
300         for(i=0;i<64;i++)
301             block1[i] = (random() % 512) -256;
302         if (is_idct){
303             fdct(block1);
304
305             for(i=0;i<64;i++)
306                 block1[i]>>=3;
307         }
308     break;
309     case 1:{
310     case 2:
311         block1[0] = (random() % 512) -256;
312         block1[1] = (random() % 512) -256;
313         block1[2] = (random() % 512) -256;
314         block1[3] = (random() % 512) -256;
315     }break;
316     }
317
318     if (form == MMX_PERM) {
319         for(i=0;i<64;i++)
320             block[idct_mmx_perm[i]] = block1[i];
321     } else if(form == MMX_SIMPLE_PERM) {
322         for(i=0;i<64;i++)
323             block[idct_simple_mmx_perm[i]] = block1[i];
324     } else {
325         for(i=0; i<64; i++)
326             block[i]= block1[i];
327     }
328
329     ti = gettime();
330     it1 = 0;
331     do {
332         for(it=0;it<NB_ITS_SPEED;it++) {
333             for(i=0; i<64; i++)
334                 block[i]= block1[i];
335 //            memcpy(block, block1, sizeof(DCTELEM) * 64);
336 // dont memcpy especially not fastmemcpy because it does movntq !!!
337             fdct_func(block);
338         }
339         it1 += NB_ITS_SPEED;
340         ti1 = gettime() - ti;
341     } while (ti1 < 1000000);
342     emms_c();
343
344     printf("%s %s: %0.1f kdct/s\n",
345            is_idct ? "IDCT" : "DCT",
346            name, (double)it1 * 1000.0 / (double)ti1);
347 #endif
348 }
349
350 static uint8_t img_dest[64] __attribute__ ((aligned (8)));
351 static uint8_t img_dest1[64] __attribute__ ((aligned (8)));
352
353 void idct248_ref(uint8_t *dest, int linesize, int16_t *block)
354 {
355     static int init;
356     static double c8[8][8];
357     static double c4[4][4];
358     double block1[64], block2[64], block3[64];
359     double s, sum, v;
360     int i, j, k;
361
362     if (!init) {
363         init = 1;
364
365         for(i=0;i<8;i++) {
366             sum = 0;
367             for(j=0;j<8;j++) {
368                 s = (i==0) ? sqrt(1.0/8.0) : sqrt(1.0/4.0);
369                 c8[i][j] = s * cos(M_PI * i * (j + 0.5) / 8.0);
370                 sum += c8[i][j] * c8[i][j];
371             }
372         }
373
374         for(i=0;i<4;i++) {
375             sum = 0;
376             for(j=0;j<4;j++) {
377                 s = (i==0) ? sqrt(1.0/4.0) : sqrt(1.0/2.0);
378                 c4[i][j] = s * cos(M_PI * i * (j + 0.5) / 4.0);
379                 sum += c4[i][j] * c4[i][j];
380             }
381         }
382     }
383
384     /* butterfly */
385     s = 0.5 * sqrt(2.0);
386     for(i=0;i<4;i++) {
387         for(j=0;j<8;j++) {
388             block1[8*(2*i)+j] = (block[8*(2*i)+j] + block[8*(2*i+1)+j]) * s;
389             block1[8*(2*i+1)+j] = (block[8*(2*i)+j] - block[8*(2*i+1)+j]) * s;
390         }
391     }
392
393     /* idct8 on lines */
394     for(i=0;i<8;i++) {
395         for(j=0;j<8;j++) {
396             sum = 0;
397             for(k=0;k<8;k++)
398                 sum += c8[k][j] * block1[8*i+k];
399             block2[8*i+j] = sum;
400         }
401     }
402
403     /* idct4 */
404     for(i=0;i<8;i++) {
405         for(j=0;j<4;j++) {
406             /* top */
407             sum = 0;
408             for(k=0;k<4;k++)
409                 sum += c4[k][j] * block2[8*(2*k)+i];
410             block3[8*(2*j)+i] = sum;
411
412             /* bottom */
413             sum = 0;
414             for(k=0;k<4;k++)
415                 sum += c4[k][j] * block2[8*(2*k+1)+i];
416             block3[8*(2*j+1)+i] = sum;
417         }
418     }
419
420     /* clamp and store the result */
421     for(i=0;i<8;i++) {
422         for(j=0;j<8;j++) {
423             v = block3[8*i+j];
424             if (v < 0)
425                 v = 0;
426             else if (v > 255)
427                 v = 255;
428             dest[i * linesize + j] = (int)rint(v);
429         }
430     }
431 }
432
433 void idct248_error(const char *name,
434                     void (*idct248_put)(uint8_t *dest, int line_size, int16_t *block))
435 {
436     int it, i, it1, ti, ti1, err_max, v;
437
438     srandom(0);
439
440     /* just one test to see if code is correct (precision is less
441        important here) */
442     err_max = 0;
443     for(it=0;it<NB_ITS;it++) {
444
445         /* XXX: use forward transform to generate values */
446         for(i=0;i<64;i++)
447             block1[i] = (random() % 256) - 128;
448         block1[0] += 1024;
449
450         for(i=0; i<64; i++)
451             block[i]= block1[i];
452         idct248_ref(img_dest1, 8, block);
453
454         for(i=0; i<64; i++)
455             block[i]= block1[i];
456         idct248_put(img_dest, 8, block);
457
458         for(i=0;i<64;i++) {
459             v = abs((int)img_dest[i] - (int)img_dest1[i]);
460             if (v == 255)
461                 printf("%d %d\n", img_dest[i], img_dest1[i]);
462             if (v > err_max)
463                 err_max = v;
464         }
465 #if 0
466         printf("ref=\n");
467         for(i=0;i<8;i++) {
468             int j;
469             for(j=0;j<8;j++) {
470                 printf(" %3d", img_dest1[i*8+j]);
471             }
472             printf("\n");
473         }
474
475         printf("out=\n");
476         for(i=0;i<8;i++) {
477             int j;
478             for(j=0;j<8;j++) {
479                 printf(" %3d", img_dest[i*8+j]);
480             }
481             printf("\n");
482         }
483 #endif
484     }
485     printf("%s %s: err_inf=%d\n",
486            1 ? "IDCT248" : "DCT248",
487            name, err_max);
488
489     ti = gettime();
490     it1 = 0;
491     do {
492         for(it=0;it<NB_ITS_SPEED;it++) {
493             for(i=0; i<64; i++)
494                 block[i]= block1[i];
495 //            memcpy(block, block1, sizeof(DCTELEM) * 64);
496 // dont memcpy especially not fastmemcpy because it does movntq !!!
497             idct248_put(img_dest, 8, block);
498         }
499         it1 += NB_ITS_SPEED;
500         ti1 = gettime() - ti;
501     } while (ti1 < 1000000);
502     emms_c();
503
504     printf("%s %s: %0.1f kdct/s\n",
505            1 ? "IDCT248" : "DCT248",
506            name, (double)it1 * 1000.0 / (double)ti1);
507 }
508
509 void help(void)
510 {
511     printf("dct-test [-i] [<test-number>]\n"
512            "test-number 0 -> test with random matrixes\n"
513            "            1 -> test with random sparse matrixes\n"
514            "            2 -> do 3. test from mpeg4 std\n"
515            "-i          test IDCT implementations\n"
516            "-4          test IDCT248 implementations\n");
517 }
518
519 int main(int argc, char **argv)
520 {
521     int test_idct = 0, test_248_dct = 0;
522     int c,i;
523     int test=1;
524
525     init_fdct();
526     idct_mmx_init();
527
528     for(i=0;i<256;i++) cropTbl[i + MAX_NEG_CROP] = i;
529     for(i=0;i<MAX_NEG_CROP;i++) {
530         cropTbl[i] = 0;
531         cropTbl[i + MAX_NEG_CROP + 256] = 255;
532     }
533
534     for(;;) {
535         c = getopt(argc, argv, "ih4");
536         if (c == -1)
537             break;
538         switch(c) {
539         case 'i':
540             test_idct = 1;
541             break;
542         case '4':
543             test_248_dct = 1;
544             break;
545         default :
546         case 'h':
547             help();
548             return 0;
549         }
550     }
551
552     if(optind <argc) test= atoi(argv[optind]);
553
554     printf("ffmpeg DCT/IDCT test\n");
555
556     if (test_248_dct) {
557         idct248_error("SIMPLE-C", simple_idct248_put);
558     } else {
559       for (i=0;algos[i].name;i++)
560         if (algos[i].is_idct == test_idct) {
561           dct_error (algos[i].name, algos[i].is_idct, algos[i].func, algos[i].ref, algos[i].format, test);
562         }
563     }
564     return 0;
565 }