[libav.git] / libavcodec / dct-test.c

/*
 * (c) 2001 Fabrice Bellard
 *     2007 Marc Hoffman <marc.hoffman@analog.com>
 *
 * This file is part of FFmpeg.
 *
 * FFmpeg is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 * FFmpeg is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with FFmpeg; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
 */

/**
 * @file dct-test.c
 * DCT test. (c) 2001 Fabrice Bellard.
 * Started from sample code by Juan J. Sierralta P.
 */

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <sys/time.h>
#include <unistd.h>
#include <math.h>

#include "libavutil/common.h"
#include "dsputil.h"

#include "simple_idct.h"
#include "faandct.h"
#include "faanidct.h"
#include "i386/idct_xvid.h"

#undef printf
#undef random

void *fast_memcpy(void *a, const void *b, size_t c){return memcpy(a,b,c);};

/* reference fdct/idct */
extern void fdct(DCTELEM *block);
extern void idct(DCTELEM *block);
extern void init_fdct();

extern void ff_mmx_idct(DCTELEM *data);
extern void ff_mmxext_idct(DCTELEM *data);

extern void odivx_idct_c (short *block);

// BFIN
extern void ff_bfin_idct (DCTELEM *block)  ;
extern void ff_bfin_fdct (DCTELEM *block) ;

// ALTIVEC
extern void fdct_altivec (DCTELEM *block);
//extern void idct_altivec (DCTELEM *block);?? no routine


struct algo {
  const char *name;
  enum { FDCT, IDCT } is_idct;
  void (* func) (DCTELEM *block);
  void (* ref)  (DCTELEM *block);
  enum formattag { NO_PERM,MMX_PERM, MMX_SIMPLE_PERM, SCALE_PERM, SSE2_PERM } format;
  int  mm_support;
};

#ifndef FAAN_POSTSCALE
#define FAAN_SCALE SCALE_PERM
#else
#define FAAN_SCALE NO_PERM
#endif

struct algo algos[] = {
  {"REF-DBL",         0, fdct,               fdct, NO_PERM},
  {"FAAN",            0, ff_faandct,         fdct, FAAN_SCALE},
  {"FAANI",           1, ff_faanidct,        idct, NO_PERM},
  {"IJG-AAN-INT",     0, fdct_ifast,         fdct, SCALE_PERM},
  {"IJG-LLM-INT",     0, ff_jpeg_fdct_islow, fdct, NO_PERM},
  {"REF-DBL",         1, idct,               idct, NO_PERM},
  {"INT",             1, j_rev_dct,          idct, MMX_PERM},
  {"SIMPLE-C",        1, ff_simple_idct,     idct, NO_PERM},

#ifdef HAVE_MMX
  {"MMX",             0, ff_fdct_mmx,        fdct, NO_PERM, MM_MMX},
#ifdef HAVE_MMX2
  {"MMX2",            0, ff_fdct_mmx2,       fdct, NO_PERM, MM_MMXEXT},
#endif

#ifdef CONFIG_GPL
  {"LIBMPEG2-MMX",    1, ff_mmx_idct,        idct, MMX_PERM, MM_MMX},
  {"LIBMPEG2-MMXEXT", 1, ff_mmxext_idct,     idct, MMX_PERM, MM_MMXEXT},
#endif
  {"SIMPLE-MMX",      1, ff_simple_idct_mmx, idct, MMX_SIMPLE_PERM, MM_MMX},
  {"XVID-MMX",        1, ff_idct_xvid_mmx,   idct, NO_PERM, MM_MMX},
  {"XVID-MMX2",       1, ff_idct_xvid_mmx2,  idct, NO_PERM, MM_MMXEXT},
  {"XVID-SSE2",       1, ff_idct_xvid_sse2,  idct, SSE2_PERM, MM_SSE2},
#endif

#ifdef HAVE_ALTIVEC
  {"altivecfdct",     0, fdct_altivec,       fdct, NO_PERM, MM_ALTIVEC},
#endif

#ifdef ARCH_BFIN
  {"BFINfdct",        0, ff_bfin_fdct,       fdct, NO_PERM},
  {"BFINidct",        1, ff_bfin_idct,       idct, NO_PERM},
#endif

  { 0 }
};

#define AANSCALE_BITS 12
static const unsigned short aanscales[64] = {
    /* precomputed values scaled up by 14 bits */
    16384, 22725, 21407, 19266, 16384, 12873,  8867,  4520,
    22725, 31521, 29692, 26722, 22725, 17855, 12299,  6270,
    21407, 29692, 27969, 25172, 21407, 16819, 11585,  5906,
    19266, 26722, 25172, 22654, 19266, 15137, 10426,  5315,
    16384, 22725, 21407, 19266, 16384, 12873,  8867,  4520,
    12873, 17855, 16819, 15137, 12873, 10114,  6967,  3552,
    8867, 12299, 11585, 10426,  8867,  6967,  4799,  2446,
    4520,  6270,  5906,  5315,  4520,  3552,  2446,  1247
};

uint8_t cropTbl[256 + 2 * MAX_NEG_CROP];

int64_t gettime(void)
{
    struct timeval tv;
    gettimeofday(&tv,NULL);
    return (int64_t)tv.tv_sec * 1000000 + tv.tv_usec;
}

#define NB_ITS 20000
#define NB_ITS_SPEED 50000

static short idct_mmx_perm[64];

static short idct_simple_mmx_perm[64]={
        0x00, 0x08, 0x04, 0x09, 0x01, 0x0C, 0x05, 0x0D,
        0x10, 0x18, 0x14, 0x19, 0x11, 0x1C, 0x15, 0x1D,
        0x20, 0x28, 0x24, 0x29, 0x21, 0x2C, 0x25, 0x2D,
        0x12, 0x1A, 0x16, 0x1B, 0x13, 0x1E, 0x17, 0x1F,
        0x02, 0x0A, 0x06, 0x0B, 0x03, 0x0E, 0x07, 0x0F,
        0x30, 0x38, 0x34, 0x39, 0x31, 0x3C, 0x35, 0x3D,
        0x22, 0x2A, 0x26, 0x2B, 0x23, 0x2E, 0x27, 0x2F,
        0x32, 0x3A, 0x36, 0x3B, 0x33, 0x3E, 0x37, 0x3F,
};

static const uint8_t idct_sse2_row_perm[8] = {0, 4, 1, 5, 2, 6, 3, 7};

void idct_mmx_init(void)
{
    int i;

    /* the mmx/mmxext idct uses a reordered input, so we patch scan tables */
    for (i = 0; i < 64; i++) {
        idct_mmx_perm[i] = (i & 0x38) | ((i & 6) >> 1) | ((i & 1) << 2);
//        idct_simple_mmx_perm[i] = simple_block_permute_op(i);
    }
}

static DCTELEM block[64] __attribute__ ((aligned (16)));
static DCTELEM block1[64] __attribute__ ((aligned (8)));
static DCTELEM block_org[64] __attribute__ ((aligned (8)));

void dct_error(const char *name, int is_idct,
               void (*fdct_func)(DCTELEM *block),
               void (*fdct_ref)(DCTELEM *block), int form, int test)
{
    int it, i, scale;
    int err_inf, v;
    int64_t err2, ti, ti1, it1;
    int64_t sysErr[64], sysErrMax=0;
    int maxout=0;
    int blockSumErrMax=0, blockSumErr;

    srandom(0);

    err_inf = 0;
    err2 = 0;
    for(i=0; i<64; i++) sysErr[i]=0;
    for(it=0;it<NB_ITS;it++) {
        for(i=0;i<64;i++)
            block1[i] = 0;
        switch(test){
        case 0:
            for(i=0;i<64;i++)
                block1[i] = (random() % 512) -256;
            if (is_idct){
                fdct(block1);

                for(i=0;i<64;i++)
                    block1[i]>>=3;
            }
        break;
        case 1:{
            int num= (random()%10)+1;
            for(i=0;i<num;i++)
                block1[random()%64] = (random() % 512) -256;
        }break;
        case 2:
            block1[0]= (random()%4096)-2048;
            block1[63]= (block1[0]&1)^1;
        break;
        }

#if 0 // simulate mismatch control
{ int sum=0;
        for(i=0;i<64;i++)
           sum+=block1[i];

        if((sum&1)==0) block1[63]^=1;
}
#endif

        for(i=0; i<64; i++)
            block_org[i]= block1[i];

        if (form == MMX_PERM) {
            for(i=0;i<64;i++)
                block[idct_mmx_perm[i]] = block1[i];
            } else if (form == MMX_SIMPLE_PERM) {
            for(i=0;i<64;i++)
                block[idct_simple_mmx_perm[i]] = block1[i];

        } else if (form == SSE2_PERM) {
            for(i=0; i<64; i++)
                block[(i&0x38) | idct_sse2_row_perm[i&7]] = block1[i];
        } else {
            for(i=0; i<64; i++)
                block[i]= block1[i];
        }
#if 0 // simulate mismatch control for tested IDCT but not the ref
{ int sum=0;
        for(i=0;i<64;i++)
           sum+=block[i];

        if((sum&1)==0) block[63]^=1;
}
#endif

        fdct_func(block);
        emms_c(); /* for ff_mmx_idct */

        if (form == SCALE_PERM) {
            for(i=0; i<64; i++) {
                scale = 8*(1 << (AANSCALE_BITS + 11)) / aanscales[i];
                block[i] = (block[i] * scale /*+ (1<<(AANSCALE_BITS-1))*/) >> AANSCALE_BITS;
            }
        }

        fdct_ref(block1);

        blockSumErr=0;
        for(i=0;i<64;i++) {
            v = abs(block[i] - block1[i]);
            if (v > err_inf)
                err_inf = v;
            err2 += v * v;
            sysErr[i] += block[i] - block1[i];
            blockSumErr += v;
            if( abs(block[i])>maxout) maxout=abs(block[i]);
        }
        if(blockSumErrMax < blockSumErr) blockSumErrMax= blockSumErr;
#if 0 // print different matrix pairs
        if(blockSumErr){
            printf("\n");
            for(i=0; i<64; i++){
                if((i&7)==0) printf("\n");
                printf("%4d ", block_org[i]);
            }
            for(i=0; i<64; i++){
                if((i&7)==0) printf("\n");
                printf("%4d ", block[i] - block1[i]);
            }
        }
#endif
    }
    for(i=0; i<64; i++) sysErrMax= FFMAX(sysErrMax, FFABS(sysErr[i]));

#if 1 // dump systematic errors
    for(i=0; i<64; i++){
        if(i%8==0) printf("\n");
        printf("%5d ", (int)sysErr[i]);
    }
    printf("\n");
#endif

    printf("%s %s: err_inf=%d err2=%0.8f syserr=%0.8f maxout=%d blockSumErr=%d\n",
           is_idct ? "IDCT" : "DCT",
           name, err_inf, (double)err2 / NB_ITS / 64.0, (double)sysErrMax / NB_ITS, maxout, blockSumErrMax);
#if 1 //Speed test
    /* speed test */
    for(i=0;i<64;i++)
        block1[i] = 0;
    switch(test){
    case 0:
        for(i=0;i<64;i++)
            block1[i] = (random() % 512) -256;
        if (is_idct){
            fdct(block1);

            for(i=0;i<64;i++)
                block1[i]>>=3;
        }
    break;
    case 1:{
    case 2:
        block1[0] = (random() % 512) -256;
        block1[1] = (random() % 512) -256;
        block1[2] = (random() % 512) -256;
        block1[3] = (random() % 512) -256;
    }break;
    }

    if (form == MMX_PERM) {
        for(i=0;i<64;i++)
            block[idct_mmx_perm[i]] = block1[i];
    } else if(form == MMX_SIMPLE_PERM) {
        for(i=0;i<64;i++)
            block[idct_simple_mmx_perm[i]] = block1[i];
    } else {
        for(i=0; i<64; i++)
            block[i]= block1[i];
    }

    ti = gettime();
    it1 = 0;
    do {
        for(it=0;it<NB_ITS_SPEED;it++) {
            for(i=0; i<64; i++)
                block[i]= block1[i];
//            memcpy(block, block1, sizeof(DCTELEM) * 64);
// do not memcpy especially not fastmemcpy because it does movntq !!!
            fdct_func(block);
        }
        it1 += NB_ITS_SPEED;
        ti1 = gettime() - ti;
    } while (ti1 < 1000000);
    emms_c();

    printf("%s %s: %0.1f kdct/s\n",
           is_idct ? "IDCT" : "DCT",
           name, (double)it1 * 1000.0 / (double)ti1);
#endif
}

static uint8_t img_dest[64] __attribute__ ((aligned (8)));
static uint8_t img_dest1[64] __attribute__ ((aligned (8)));

void idct248_ref(uint8_t *dest, int linesize, int16_t *block)
{
    static int init;
    static double c8[8][8];
    static double c4[4][4];
    double block1[64], block2[64], block3[64];
    double s, sum, v;
    int i, j, k;

    if (!init) {
        init = 1;

        for(i=0;i<8;i++) {
            sum = 0;
            for(j=0;j<8;j++) {
                s = (i==0) ? sqrt(1.0/8.0) : sqrt(1.0/4.0);
                c8[i][j] = s * cos(M_PI * i * (j + 0.5) / 8.0);
                sum += c8[i][j] * c8[i][j];
            }
        }

        for(i=0;i<4;i++) {
            sum = 0;
            for(j=0;j<4;j++) {
                s = (i==0) ? sqrt(1.0/4.0) : sqrt(1.0/2.0);
                c4[i][j] = s * cos(M_PI * i * (j + 0.5) / 4.0);
                sum += c4[i][j] * c4[i][j];
            }
        }
    }

    /* butterfly */
    s = 0.5 * sqrt(2.0);
    for(i=0;i<4;i++) {
        for(j=0;j<8;j++) {
            block1[8*(2*i)+j] = (block[8*(2*i)+j] + block[8*(2*i+1)+j]) * s;
            block1[8*(2*i+1)+j] = (block[8*(2*i)+j] - block[8*(2*i+1)+j]) * s;
        }
    }

    /* idct8 on lines */
    for(i=0;i<8;i++) {
        for(j=0;j<8;j++) {
            sum = 0;
            for(k=0;k<8;k++)
                sum += c8[k][j] * block1[8*i+k];
            block2[8*i+j] = sum;
        }
    }

    /* idct4 */
    for(i=0;i<8;i++) {
        for(j=0;j<4;j++) {
            /* top */
            sum = 0;
            for(k=0;k<4;k++)
                sum += c4[k][j] * block2[8*(2*k)+i];
            block3[8*(2*j)+i] = sum;

            /* bottom */
            sum = 0;
            for(k=0;k<4;k++)
                sum += c4[k][j] * block2[8*(2*k+1)+i];
            block3[8*(2*j+1)+i] = sum;
        }
    }

    /* clamp and store the result */
    for(i=0;i<8;i++) {
        for(j=0;j<8;j++) {
            v = block3[8*i+j];
            if (v < 0)
                v = 0;
            else if (v > 255)
                v = 255;
            dest[i * linesize + j] = (int)rint(v);
        }
    }
}

void idct248_error(const char *name,
                    void (*idct248_put)(uint8_t *dest, int line_size, int16_t *block))
{
    int it, i, it1, ti, ti1, err_max, v;

    srandom(0);

    /* just one test to see if code is correct (precision is less
       important here) */
    err_max = 0;
    for(it=0;it<NB_ITS;it++) {

        /* XXX: use forward transform to generate values */
        for(i=0;i<64;i++)
            block1[i] = (random() % 256) - 128;
        block1[0] += 1024;

        for(i=0; i<64; i++)
            block[i]= block1[i];
        idct248_ref(img_dest1, 8, block);

        for(i=0; i<64; i++)
            block[i]= block1[i];
        idct248_put(img_dest, 8, block);

        for(i=0;i<64;i++) {
            v = abs((int)img_dest[i] - (int)img_dest1[i]);
            if (v == 255)
                printf("%d %d\n", img_dest[i], img_dest1[i]);
            if (v > err_max)
                err_max = v;
        }
#if 0
        printf("ref=\n");
        for(i=0;i<8;i++) {
            int j;
            for(j=0;j<8;j++) {
                printf(" %3d", img_dest1[i*8+j]);
            }
            printf("\n");
        }

        printf("out=\n");
        for(i=0;i<8;i++) {
            int j;
            for(j=0;j<8;j++) {
                printf(" %3d", img_dest[i*8+j]);
            }
            printf("\n");
        }
#endif
    }
    printf("%s %s: err_inf=%d\n",
           1 ? "IDCT248" : "DCT248",
           name, err_max);

    ti = gettime();
    it1 = 0;
    do {
        for(it=0;it<NB_ITS_SPEED;it++) {
            for(i=0; i<64; i++)
                block[i]= block1[i];
//            memcpy(block, block1, sizeof(DCTELEM) * 64);
// do not memcpy especially not fastmemcpy because it does movntq !!!
            idct248_put(img_dest, 8, block);
        }
        it1 += NB_ITS_SPEED;
        ti1 = gettime() - ti;
    } while (ti1 < 1000000);
    emms_c();

    printf("%s %s: %0.1f kdct/s\n",
           1 ? "IDCT248" : "DCT248",
           name, (double)it1 * 1000.0 / (double)ti1);
}

void help(void)
{
    printf("dct-test [-i] [<test-number>]\n"
           "test-number 0 -> test with random matrixes\n"
           "            1 -> test with random sparse matrixes\n"
           "            2 -> do 3. test from mpeg4 std\n"
           "-i          test IDCT implementations\n"
           "-4          test IDCT248 implementations\n");
}

int main(int argc, char **argv)
{
    int test_idct = 0, test_248_dct = 0;
    int c,i;
    int test=1;
    int cpu_flags = mm_support();

    init_fdct();
    idct_mmx_init();

    for(i=0;i<256;i++) cropTbl[i + MAX_NEG_CROP] = i;
    for(i=0;i<MAX_NEG_CROP;i++) {
        cropTbl[i] = 0;
        cropTbl[i + MAX_NEG_CROP + 256] = 255;
    }

    for(;;) {
        c = getopt(argc, argv, "ih4");
        if (c == -1)
            break;
        switch(c) {
        case 'i':
            test_idct = 1;
            break;
        case '4':
            test_248_dct = 1;
            break;
        default :
        case 'h':
            help();
            return 0;
        }
    }

    if(optind <argc) test= atoi(argv[optind]);

    printf("ffmpeg DCT/IDCT test\n");

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