dct-test: split table in two for idct and fdct
[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 Libav.
6  *
7  * Libav 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  * Libav 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 Libav; 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
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 #include <math.h>
34
35 #include "libavutil/cpu.h"
36 #include "libavutil/common.h"
37 #include "libavutil/lfg.h"
38
39 #include "simple_idct.h"
40 #include "aandcttab.h"
41 #include "faandct.h"
42 #include "faanidct.h"
43 #include "x86/idct_xvid.h"
44 #include "dctref.h"
45
46 #undef printf
47
48 void ff_mmx_idct(DCTELEM *data);
49 void ff_mmxext_idct(DCTELEM *data);
50
51 void odivx_idct_c(short *block);
52
53 // BFIN
54 void ff_bfin_idct(DCTELEM *block);
55 void ff_bfin_fdct(DCTELEM *block);
56
57 // ALTIVEC
58 void fdct_altivec(DCTELEM *block);
59 //void idct_altivec(DCTELEM *block);?? no routine
60
61 // ARM
62 void ff_j_rev_dct_arm(DCTELEM *data);
63 void ff_simple_idct_arm(DCTELEM *data);
64 void ff_simple_idct_armv5te(DCTELEM *data);
65 void ff_simple_idct_armv6(DCTELEM *data);
66 void ff_simple_idct_neon(DCTELEM *data);
67
68 void ff_simple_idct_axp(DCTELEM *data);
69
70 struct algo {
71     const char *name;
72     enum { FDCT, IDCT } is_idct;
73     void (*func)(DCTELEM *block);
74     void (*ref) (DCTELEM *block);
75     enum formattag { NO_PERM, MMX_PERM, MMX_SIMPLE_PERM, SCALE_PERM,
76                      SSE2_PERM, PARTTRANS_PERM } format;
77     int mm_support;
78 };
79
80 #ifndef FAAN_POSTSCALE
81 #define FAAN_SCALE SCALE_PERM
82 #else
83 #define FAAN_SCALE NO_PERM
84 #endif
85
86 static int cpu_flags;
87
88 static const struct algo fdct_tab[] = {
89     {"REF-DBL",         0, ff_ref_fdct,        ff_ref_fdct, NO_PERM},
90     {"FAAN",            0, ff_faandct,         ff_ref_fdct, FAAN_SCALE},
91     {"IJG-AAN-INT",     0, fdct_ifast,         ff_ref_fdct, SCALE_PERM},
92     {"IJG-LLM-INT",     0, ff_jpeg_fdct_islow, ff_ref_fdct, NO_PERM},
93
94 #if HAVE_MMX
95     {"MMX",             0, ff_fdct_mmx,        ff_ref_fdct, NO_PERM, AV_CPU_FLAG_MMX},
96     {"MMX2",            0, ff_fdct_mmx2,       ff_ref_fdct, NO_PERM, AV_CPU_FLAG_MMX2},
97     {"SSE2",            0, ff_fdct_sse2,       ff_ref_fdct, NO_PERM, AV_CPU_FLAG_SSE2},
98 #endif
99
100 #if HAVE_ALTIVEC
101     {"altivecfdct",     0, fdct_altivec,       ff_ref_fdct, NO_PERM, AV_CPU_FLAG_ALTIVEC},
102 #endif
103
104 #if ARCH_BFIN
105     {"BFINfdct",        0, ff_bfin_fdct,       ff_ref_fdct, NO_PERM},
106 #endif
107
108     { 0 }
109 };
110
111 static const struct algo idct_tab[] = {
112     {"FAANI",           1, ff_faanidct,        ff_ref_idct, NO_PERM},
113     {"REF-DBL",         1, ff_ref_idct,        ff_ref_idct, NO_PERM},
114     {"INT",             1, j_rev_dct,          ff_ref_idct, MMX_PERM},
115     {"SIMPLE-C",        1, ff_simple_idct,     ff_ref_idct, NO_PERM},
116
117 #if HAVE_MMX
118 #if CONFIG_GPL
119     {"LIBMPEG2-MMX",    1, ff_mmx_idct,        ff_ref_idct, MMX_PERM, AV_CPU_FLAG_MMX},
120     {"LIBMPEG2-MMX2",   1, ff_mmxext_idct,     ff_ref_idct, MMX_PERM, AV_CPU_FLAG_MMX2},
121 #endif
122     {"SIMPLE-MMX",      1, ff_simple_idct_mmx, ff_ref_idct, MMX_SIMPLE_PERM, AV_CPU_FLAG_MMX},
123     {"XVID-MMX",        1, ff_idct_xvid_mmx,   ff_ref_idct, NO_PERM, AV_CPU_FLAG_MMX},
124     {"XVID-MMX2",       1, ff_idct_xvid_mmx2,  ff_ref_idct, NO_PERM, AV_CPU_FLAG_MMX2},
125     {"XVID-SSE2",       1, ff_idct_xvid_sse2,  ff_ref_idct, SSE2_PERM, AV_CPU_FLAG_SSE2},
126 #endif
127
128 #if ARCH_BFIN
129     {"BFINidct",        1, ff_bfin_idct,       ff_ref_idct, NO_PERM},
130 #endif
131
132 #if ARCH_ARM
133     {"SIMPLE-ARM",      1, ff_simple_idct_arm, ff_ref_idct, NO_PERM },
134     {"INT-ARM",         1, ff_j_rev_dct_arm,   ff_ref_idct, MMX_PERM },
135 #endif
136 #if HAVE_ARMV5TE
137     {"SIMPLE-ARMV5TE",  1, ff_simple_idct_armv5te, ff_ref_idct, NO_PERM },
138 #endif
139 #if HAVE_ARMV6
140     {"SIMPLE-ARMV6",    1, ff_simple_idct_armv6, ff_ref_idct, MMX_PERM },
141 #endif
142 #if HAVE_NEON
143     {"SIMPLE-NEON",     1, ff_simple_idct_neon, ff_ref_idct, PARTTRANS_PERM },
144 #endif
145
146 #if ARCH_ALPHA
147     {"SIMPLE-ALPHA",    1, ff_simple_idct_axp,  ff_ref_idct, NO_PERM },
148 #endif
149
150     { 0 }
151 };
152
153 #define AANSCALE_BITS 12
154
155 uint8_t cropTbl[256 + 2 * MAX_NEG_CROP];
156
157 static int64_t gettime(void)
158 {
159     struct timeval tv;
160     gettimeofday(&tv, NULL);
161     return (int64_t)tv.tv_sec * 1000000 + tv.tv_usec;
162 }
163
164 #define NB_ITS 20000
165 #define NB_ITS_SPEED 50000
166
167 static short idct_mmx_perm[64];
168
169 static short idct_simple_mmx_perm[64] = {
170     0x00, 0x08, 0x04, 0x09, 0x01, 0x0C, 0x05, 0x0D,
171     0x10, 0x18, 0x14, 0x19, 0x11, 0x1C, 0x15, 0x1D,
172     0x20, 0x28, 0x24, 0x29, 0x21, 0x2C, 0x25, 0x2D,
173     0x12, 0x1A, 0x16, 0x1B, 0x13, 0x1E, 0x17, 0x1F,
174     0x02, 0x0A, 0x06, 0x0B, 0x03, 0x0E, 0x07, 0x0F,
175     0x30, 0x38, 0x34, 0x39, 0x31, 0x3C, 0x35, 0x3D,
176     0x22, 0x2A, 0x26, 0x2B, 0x23, 0x2E, 0x27, 0x2F,
177     0x32, 0x3A, 0x36, 0x3B, 0x33, 0x3E, 0x37, 0x3F,
178 };
179
180 static const uint8_t idct_sse2_row_perm[8] = { 0, 4, 1, 5, 2, 6, 3, 7 };
181
182 static void idct_mmx_init(void)
183 {
184     int i;
185
186     /* the mmx/mmxext idct uses a reordered input, so we patch scan tables */
187     for (i = 0; i < 64; i++) {
188         idct_mmx_perm[i] = (i & 0x38) | ((i & 6) >> 1) | ((i & 1) << 2);
189     }
190 }
191
192 DECLARE_ALIGNED(16, static DCTELEM, block)[64];
193 DECLARE_ALIGNED(8,  static DCTELEM, block1)[64];
194 DECLARE_ALIGNED(8,  static DCTELEM, block_org)[64];
195
196 static inline void mmx_emms(void)
197 {
198 #if HAVE_MMX
199     if (cpu_flags & AV_CPU_FLAG_MMX)
200         __asm__ volatile ("emms\n\t");
201 #endif
202 }
203
204 static void dct_error(const struct algo *dct, int test)
205 {
206     int it, i, scale;
207     int err_inf, v;
208     int64_t err2, ti, ti1, it1;
209     int64_t sysErr[64], sysErrMax = 0;
210     int maxout = 0;
211     int blockSumErrMax = 0, blockSumErr;
212     AVLFG prng;
213
214     av_lfg_init(&prng, 1);
215
216     err_inf = 0;
217     err2 = 0;
218     for (i = 0; i < 64; i++)
219         sysErr[i] = 0;
220     for (it = 0; it < NB_ITS; it++) {
221         for (i = 0; i < 64; i++)
222             block1[i] = 0;
223         switch (test) {
224         case 0:
225             for (i = 0; i < 64; i++)
226                 block1[i] = (av_lfg_get(&prng) % 512) - 256;
227             if (dct->is_idct) {
228                 ff_ref_fdct(block1);
229                 for (i = 0; i < 64; i++)
230                     block1[i] >>= 3;
231             }
232             break;
233         case 1: {
234                 int num = av_lfg_get(&prng) % 10 + 1;
235                 for (i = 0; i < num; i++)
236                     block1[av_lfg_get(&prng) % 64] =
237                         av_lfg_get(&prng) % 512 - 256;
238             }
239             break;
240         case 2:
241             block1[0] = av_lfg_get(&prng) % 4096 - 2048;
242             block1[63] = (block1[0] & 1) ^ 1;
243             break;
244         }
245
246         for (i = 0; i < 64; i++)
247             block_org[i] = block1[i];
248
249         if (dct->format == MMX_PERM) {
250             for (i = 0; i < 64; i++)
251                 block[idct_mmx_perm[i]] = block1[i];
252         } else if (dct->format == MMX_SIMPLE_PERM) {
253             for (i = 0; i < 64; i++)
254                 block[idct_simple_mmx_perm[i]] = block1[i];
255         } else if (dct->format == SSE2_PERM) {
256             for (i = 0; i < 64; i++)
257                 block[(i & 0x38) | idct_sse2_row_perm[i & 7]] = block1[i];
258         } else if (dct->format == PARTTRANS_PERM) {
259             for (i = 0; i < 64; i++)
260                 block[(i & 0x24) | ((i & 3) << 3) | ((i >> 3) & 3)] = block1[i];
261         } else {
262             for (i = 0; i < 64; i++)
263                 block[i] = block1[i];
264         }
265
266         dct->func(block);
267         mmx_emms();
268
269         if (dct->format == SCALE_PERM) {
270             for (i = 0; i < 64; i++) {
271                 scale = 8 * (1 << (AANSCALE_BITS + 11)) / ff_aanscales[i];
272                 block[i] = (block[i] * scale) >> AANSCALE_BITS;
273             }
274         }
275
276         dct->ref(block1);
277
278         blockSumErr = 0;
279         for (i = 0; i < 64; i++) {
280             v = abs(block[i] - block1[i]);
281             if (v > err_inf)
282                 err_inf = v;
283             err2 += v * v;
284             sysErr[i] += block[i] - block1[i];
285             blockSumErr += v;
286             if (abs(block[i]) > maxout)
287                 maxout = abs(block[i]);
288         }
289         if (blockSumErrMax < blockSumErr)
290             blockSumErrMax = blockSumErr;
291     }
292     for (i = 0; i < 64; i++)
293         sysErrMax = FFMAX(sysErrMax, FFABS(sysErr[i]));
294
295     for (i = 0; i < 64; i++) {
296         if (i % 8 == 0)
297             printf("\n");
298         printf("%7d ", (int) sysErr[i]);
299     }
300     printf("\n");
301
302     printf("%s %s: err_inf=%d err2=%0.8f syserr=%0.8f maxout=%d blockSumErr=%d\n",
303            dct->is_idct ? "IDCT" : "DCT", dct->name, err_inf,
304            (double) err2 / NB_ITS / 64.0, (double) sysErrMax / NB_ITS,
305            maxout, blockSumErrMax);
306
307     /* speed test */
308     for (i = 0; i < 64; i++)
309         block1[i] = 0;
310
311     switch (test) {
312     case 0:
313         for (i = 0; i < 64; i++)
314             block1[i] = av_lfg_get(&prng) % 512 - 256;
315         if (dct->is_idct) {
316             ff_ref_fdct(block1);
317             for (i = 0; i < 64; i++)
318                 block1[i] >>= 3;
319         }
320         break;
321     case 1:
322     case 2:
323         block1[0] = av_lfg_get(&prng) % 512 - 256;
324         block1[1] = av_lfg_get(&prng) % 512 - 256;
325         block1[2] = av_lfg_get(&prng) % 512 - 256;
326         block1[3] = av_lfg_get(&prng) % 512 - 256;
327         break;
328     }
329
330     if (dct->format == MMX_PERM) {
331         for (i = 0; i < 64; i++)
332             block[idct_mmx_perm[i]] = block1[i];
333     } else if (dct->format == MMX_SIMPLE_PERM) {
334         for (i = 0; i < 64; i++)
335             block[idct_simple_mmx_perm[i]] = block1[i];
336     } else {
337         for (i = 0; i < 64; i++)
338             block[i] = block1[i];
339     }
340
341     ti = gettime();
342     it1 = 0;
343     do {
344         for (it = 0; it < NB_ITS_SPEED; it++) {
345             for (i = 0; i < 64; i++)
346                 block[i] = block1[i];
347             dct->func(block);
348         }
349         it1 += NB_ITS_SPEED;
350         ti1 = gettime() - ti;
351     } while (ti1 < 1000000);
352     mmx_emms();
353
354     printf("%s %s: %0.1f kdct/s\n", dct->is_idct ? "IDCT" : "DCT", dct->name,
355            (double) it1 * 1000.0 / (double) ti1);
356 }
357
358 DECLARE_ALIGNED(8, static uint8_t, img_dest)[64];
359 DECLARE_ALIGNED(8, static uint8_t, img_dest1)[64];
360
361 static void idct248_ref(uint8_t *dest, int linesize, int16_t *block)
362 {
363     static int init;
364     static double c8[8][8];
365     static double c4[4][4];
366     double block1[64], block2[64], block3[64];
367     double s, sum, v;
368     int i, j, k;
369
370     if (!init) {
371         init = 1;
372
373         for (i = 0; i < 8; i++) {
374             sum = 0;
375             for (j = 0; j < 8; j++) {
376                 s = (i == 0) ? sqrt(1.0 / 8.0) : sqrt(1.0 / 4.0);
377                 c8[i][j] = s * cos(M_PI * i * (j + 0.5) / 8.0);
378                 sum += c8[i][j] * c8[i][j];
379             }
380         }
381
382         for (i = 0; i < 4; i++) {
383             sum = 0;
384             for (j = 0; j < 4; j++) {
385                 s = (i == 0) ? sqrt(1.0 / 4.0) : sqrt(1.0 / 2.0);
386                 c4[i][j] = s * cos(M_PI * i * (j + 0.5) / 4.0);
387                 sum += c4[i][j] * c4[i][j];
388             }
389         }
390     }
391
392     /* butterfly */
393     s = 0.5 * sqrt(2.0);
394     for (i = 0; i < 4; i++) {
395         for (j = 0; j < 8; j++) {
396             block1[8 * (2 * i) + j] =
397                 (block[8 * (2 * i) + j] + block[8 * (2 * i + 1) + j]) * s;
398             block1[8 * (2 * i + 1) + j] =
399                 (block[8 * (2 * i) + j] - block[8 * (2 * i + 1) + j]) * s;
400         }
401     }
402
403     /* idct8 on lines */
404     for (i = 0; i < 8; i++) {
405         for (j = 0; j < 8; j++) {
406             sum = 0;
407             for (k = 0; k < 8; k++)
408                 sum += c8[k][j] * block1[8 * i + k];
409             block2[8 * i + j] = sum;
410         }
411     }
412
413     /* idct4 */
414     for (i = 0; i < 8; i++) {
415         for (j = 0; j < 4; j++) {
416             /* top */
417             sum = 0;
418             for (k = 0; k < 4; k++)
419                 sum += c4[k][j] * block2[8 * (2 * k) + i];
420             block3[8 * (2 * j) + i] = sum;
421
422             /* bottom */
423             sum = 0;
424             for (k = 0; k < 4; k++)
425                 sum += c4[k][j] * block2[8 * (2 * k + 1) + i];
426             block3[8 * (2 * j + 1) + i] = sum;
427         }
428     }
429
430     /* clamp and store the result */
431     for (i = 0; i < 8; i++) {
432         for (j = 0; j < 8; j++) {
433             v = block3[8 * i + j];
434             if      (v < 0)   v = 0;
435             else if (v > 255) v = 255;
436             dest[i * linesize + j] = (int) rint(v);
437         }
438     }
439 }
440
441 static void idct248_error(const char *name,
442                           void (*idct248_put)(uint8_t *dest, int line_size,
443                                               int16_t *block))
444 {
445     int it, i, it1, ti, ti1, err_max, v;
446     AVLFG prng;
447
448     av_lfg_init(&prng, 1);
449
450     /* just one test to see if code is correct (precision is less
451        important here) */
452     err_max = 0;
453     for (it = 0; it < NB_ITS; it++) {
454         /* XXX: use forward transform to generate values */
455         for (i = 0; i < 64; i++)
456             block1[i] = av_lfg_get(&prng) % 256 - 128;
457         block1[0] += 1024;
458
459         for (i = 0; i < 64; i++)
460             block[i] = block1[i];
461         idct248_ref(img_dest1, 8, block);
462
463         for (i = 0; i < 64; i++)
464             block[i] = block1[i];
465         idct248_put(img_dest, 8, block);
466
467         for (i = 0; i < 64; i++) {
468             v = abs((int) img_dest[i] - (int) img_dest1[i]);
469             if (v == 255)
470                 printf("%d %d\n", img_dest[i], img_dest1[i]);
471             if (v > err_max)
472                 err_max = v;
473         }
474     }
475     printf("%s %s: err_inf=%d\n", 1 ? "IDCT248" : "DCT248", name, err_max);
476
477     ti = gettime();
478     it1 = 0;
479     do {
480         for (it = 0; it < NB_ITS_SPEED; it++) {
481             for (i = 0; i < 64; i++)
482                 block[i] = block1[i];
483             idct248_put(img_dest, 8, block);
484         }
485         it1 += NB_ITS_SPEED;
486         ti1 = gettime() - ti;
487     } while (ti1 < 1000000);
488     mmx_emms();
489
490     printf("%s %s: %0.1f kdct/s\n", 1 ? "IDCT248" : "DCT248", name,
491            (double) it1 * 1000.0 / (double) ti1);
492 }
493
494 static void help(void)
495 {
496     printf("dct-test [-i] [<test-number>]\n"
497            "test-number 0 -> test with random matrixes\n"
498            "            1 -> test with random sparse matrixes\n"
499            "            2 -> do 3. test from mpeg4 std\n"
500            "-i          test IDCT implementations\n"
501            "-4          test IDCT248 implementations\n");
502 }
503
504 int main(int argc, char **argv)
505 {
506     int test_idct = 0, test_248_dct = 0;
507     int c, i;
508     int test = 1;
509
510     cpu_flags = av_get_cpu_flags();
511
512     ff_ref_dct_init();
513     idct_mmx_init();
514
515     for (i = 0; i < 256; i++)
516         cropTbl[i + MAX_NEG_CROP] = i;
517     for (i = 0; i < MAX_NEG_CROP; i++) {
518         cropTbl[i] = 0;
519         cropTbl[i + MAX_NEG_CROP + 256] = 255;
520     }
521
522     for (;;) {
523         c = getopt(argc, argv, "ih4");
524         if (c == -1)
525             break;
526         switch (c) {
527         case 'i':
528             test_idct = 1;
529             break;
530         case '4':
531             test_248_dct = 1;
532             break;
533         default:
534         case 'h':
535             help();
536             return 0;
537         }
538     }
539
540     if (optind < argc)
541         test = atoi(argv[optind]);
542
543     printf("ffmpeg DCT/IDCT test\n");
544
545     if (test_248_dct) {
546         idct248_error("SIMPLE-C", ff_simple_idct248_put);
547     } else {
548         const struct algo *algos = test_idct ? idct_tab : fdct_tab;
549         for (i = 0; algos[i].name; i++)
550             if (!(~cpu_flags & algos[i].mm_support)) {
551                 dct_error(&algos[i], test);
552             }
553     }
554     return 0;
555 }