/************************************************* * libr263: fast H.263 encoder library * * Copyright (C) 1996, Roalt Aalmoes, Twente University * SPA multimedia group * * Based on Telenor TMN 1.6 encoder (Copyright (C) 1995, Telenor R&D) * created by Karl Lillevold  * * Author encoder: Roalt Aalmoes, <aalmoes@huygens.nl> *  * Date: 31-07-96 **************************************************//***************************************************************** * Some routines are translated from Gisle Bj鴑tegaards's FORTRAN * routines by Robert.Danielsen@nta.no * *****************************************************************//***************************************************************** * This source includes sources from Berkeley's MPEG-1 encoder * which are copyright of Berkeley University, California, USA *****************************************************************/#include "macros.h"#include "sim.h"#include <math.h>#ifndef PI# ifdef M_PI#  define PI M_PI# else#  define PI 3.14159265358979323846# endif#endifint	zigzag[8][8] = {  {0, 1, 5, 6,14,15,27,28},  {2, 4, 7,13,16,26,29,42},  {3, 8,12,17,25,30,41,43},  {9,11,18,24,31,40,44,53},  {10,19,23,32,39,45,52,54},  {20,22,33,38,46,51,55,60},  {21,34,37,47,50,56,59,61},  {35,36,48,49,57,58,62,63},};#ifndef FASTDCT/********************************************************************** * *	Name:		Dct *	Description:	Does dct on an 8x8 block, does zigzag-scanning of *			coefficients * *	Input:		64 pixels in a 1D array *	Returns:	64 coefficients in a 1D array *	Side effects:	 * *	Date: 930128	Author: Robert.Danielsen@nta.no * **********************************************************************/int Dct( int *block, int *coeff){  int        j1, i, j, k;  float	b[8];  float 	b1[8];  float 	d[8][8];  float 	f0=.7071068,f1=.4903926,f2=.4619398,f3=.4157348,f4=.3535534;  float 	f5=.2777851,f6=.1913417,f7=.0975452;  for (i = 0, k = 0; i < 8; i++, k += 8) {    for (j = 0; j < 8; j++) {      b[j] = block[k+j];    }    /* Horizontal transform */    for (j = 0; j < 4; j++) {      j1 = 7 - j;      b1[j] = b[j] + b[j1];      b1[j1] = b[j] - b[j1];    }    b[0] = b1[0] + b1[3];    b[1] = b1[1] + b1[2];    b[2] = b1[1] - b1[2];    b[3] = b1[0] - b1[3];    b[4] = b1[4];    b[5] = (b1[6] - b1[5]) * f0;    b[6] = (b1[6] + b1[5]) * f0;    b[7] = b1[7];    d[i][0] = (b[0] + b[1]) * f4;    d[i][4] = (b[0] - b[1]) * f4;    d[i][2] = b[2] * f6 + b[3] * f2;    d[i][6] = b[3] * f6 - b[2] * f2;    b1[4] = b[4] + b[5];    b1[7] = b[7] + b[6];    b1[5] = b[4] - b[5];    b1[6] = b[7] - b[6];    d[i][1] = b1[4] * f7 + b1[7] * f1;    d[i][5] = b1[5] * f3 + b1[6] * f5;    d[i][7] = b1[7] * f7 - b1[4] * f1;    d[i][3] = b1[6] * f3 - b1[5] * f5;  }  /* Vertical transform */  for (i = 0; i < 8; i++) {    for (j = 0; j < 4; j++) {      j1 = 7 - j;      b1[j] = d[j][i] + d[j1][i];      b1[j1] = d[j][i] - d[j1][i];    }    b[0] = b1[0] + b1[3];    b[1] = b1[1] + b1[2];    b[2] = b1[1] - b1[2];    b[3] = b1[0] - b1[3];    b[4] = b1[4];    b[5] = (b1[6] - b1[5]) * f0;    b[6] = (b1[6] + b1[5]) * f0;    b[7] = b1[7];    d[0][i] = (b[0] + b[1]) * f4;    d[4][i] = (b[0] - b[1]) * f4;    d[2][i] = b[2] * f6 + b[3] * f2;    d[6][i] = b[3] * f6 - b[2] * f2;    b1[4] = b[4] + b[5];    b1[7] = b[7] + b[6];    b1[5] = b[4] - b[5];    b1[6] = b[7] - b[6];    d[1][i] = b1[4] * f7 + b1[7] * f1;    d[5][i] = b1[5] * f3 + b1[6] * f5;    d[7][i] = b1[7] * f7 - b1[4] * f1;    d[3][i] = b1[6] * f3 - b1[5] * f5;  }  /* Zigzag - scanning */  for (i = 0; i < 8; i++) {    for (j = 0; j < 8; j++) {      *(coeff + zigzag[i][j]) = (int)(d[i][j]);    }  }  return 0;}#else				/* Start of MPEG DCT operation */typedef unsigned char uint8;typedef char int8;typedef unsigned short uint16;typedef short int16;#ifdef LONG_32		typedef unsigned long uint32;typedef long int32;#elsetypedef unsigned int uint32;typedef int int32;#endif/* this is ansi.h */#undef _ANSI_ARGS_#undef const#ifdef NON_ANSI_COMPILER#define _ANSI_ARGS_(x)       ()#define CONST#else#define _ANSI_ARGS_(x)   x#define CONST const#ifdef __cplusplus#define VARARGS (...)#else#define VARARGS ()#endif#endif#define DCTSIZE     8   /* you really don't want to change this */#define DCTSIZE_SQ 64   /* you really don't want to change this */#define DCTSIZE2    DCTSIZE*DCTSIZEtypedef short DCTELEM;typedef DCTELEM DCTBLOCK[DCTSIZE2];typedef DCTELEM DCTBLOCK_2D[DCTSIZE][DCTSIZE];/* We assume that right shift corresponds to signed division by 2 with * rounding towards minus infinity.  This is correct for typical "arithmetic * shift" instructions that shift in copies of the sign bit.  But some * C compilers implement >> with an unsigned shift.  For these machines you * must define RIGHT_SHIFT_IS_UNSIGNED. * RIGHT_SHIFT provides a proper signed right shift of an int32 quantity. * It is only applied with constant shift counts.  SHIFT_TEMPS must be * included in the variables of any routine using RIGHT_SHIFT. */#ifdef RIGHT_SHIFT_IS_UNSIGNED#define SHIFT_TEMPS     int32 shift_temp/*#define RIGHT_SHIFT(x,shft)  ((shift_temp = (x)) < 0 ? (shift_temp >> (shft)) |((~((int32) 0)) << (32-(shft))) : (shift_temp >> (shft)))     	 */   #else#define SHIFT_TEMPS#define RIGHT_SHIFT(x,shft)     ((x) >> (shft))#endif          #define LG2_DCT_SCALE 16#define ONE	((int32) 1)#define DCT_SCALE (ONE << LG2_DCT_SCALE)/* In some places we shift the inputs left by a couple more bits, *//* so that they can be added to fractional results without too much */     /* loss of precision. */#define LG2_OVERSCALE 2#define OVERSCALE  (ONE << LG2_OVERSCALE)#define OVERSHIFT(x)  ((x) <<= LG2_OVERSCALE)/* Scale a fractional constant by DCT_SCALE */#define FIX(x)	((int32) ((x) * DCT_SCALE + 0.5))/* Scale a fractional constant by DCT_SCALE/OVERSCALE *//* Such a constant can be multiplied with an overscaled input *//* to produce something that's scaled by DCT_SCALE */#define FIXO(x)  ((int32) ((x) * DCT_SCALE / OVERSCALE + 0.5))/* Descale and correctly round a value that's scaled by DCT_SCALE */#define UNFIX(x)   RIGHT_SHIFT((x) + (ONE << (LG2_DCT_SCALE-1)), LG2_DCT_SCALE)/* Same with an additional division by 2, ie, correctly rounded UNFIX(x/2) */#define UNFIXH(x)  RIGHT_SHIFT((x) + (ONE << LG2_DCT_SCALE), LG2_DCT_SCALE+1)/* Take a value scaled by DCT_SCALE and round to integer scaled by OVERSCALE */#define UNFIXO(x)  RIGHT_SHIFT((x) + (ONE << (LG2_DCT_SCALE-1-LG2_OVERSCALE)),LG2_DCT_SCALE-LG2_OVERSCALE)/* Here are the constants we need *//* SIN_i_j is sine of i*pi/j, scaled by DCT_SCALE *//* COS_i_j is cosine of i*pi/j, scaled by DCT_SCALE */#define SIN_1_4 FIX(0.707106781)#define COS_1_4 SIN_1_4#define SIN_1_8 FIX(0.382683432)#define COS_1_8 FIX(0.923879533)#define SIN_3_8 COS_1_8#define COS_3_8 SIN_1_8#define SIN_1_16 FIX(0.195090322)#define COS_1_16 FIX(0.980785280)#define SIN_7_16 COS_1_16#define COS_7_16 SIN_1_16#define SIN_3_16 FIX(0.555570233)#define COS_3_16 FIX(0.831469612)#define SIN_5_16 COS_3_16#define COS_5_16 SIN_3_16/* OSIN_i_j is sine of i*pi/j, scaled by DCT_SCALE/OVERSCALE *//* OCOS_i_j is cosine of i*pi/j, scaled by DCT_SCALE/OVERSCALE */#define OSIN_1_4 FIXO(0.707106781)#define OCOS_1_4 OSIN_1_4#define OSIN_1_8 FIXO(0.382683432)#define OCOS_1_8 FIXO(0.923879533)#define OSIN_3_8 OCOS_1_8#define OCOS_3_8 OSIN_1_8#define OSIN_1_16 FIXO(0.195090322)#define OCOS_1_16 FIXO(0.980785280)#define OSIN_7_16 OCOS_1_16#define OCOS_7_16 OSIN_1_16#define OSIN_3_16 FIXO(0.555570233)#define OCOS_3_16 FIXO(0.831469612)#define OSIN_5_16 OCOS_3_16#define OCOS_5_16 OSIN_3_16     /*==================* * TYPE DEFINITIONS * *==================*//*   *  your basic Block type */typedef int32 Block[DCTSIZE][DCTSIZE];typedef int32 FlatBlock[DCTSIZE_SQ];typedef	    int32   LumBlock[2*DCTSIZE][2*DCTSIZE];typedef	    int32   ChromBlock[DCTSIZE][DCTSIZE];/* Prototypes */void reference_fwd_dct(Block block, Block dest);void mp_fwd_dct_fast(Block data2d, Block dest2d);/*void mp_fwd_dct_fast _ANSI_ARGS_((int16 *data2d, int16 *dest2d)); */void init_fdct(void);/* * -------------------------------------------------------------- * * mp_fwd_dct_fast -- * * Perform the forward DCT on one block of samples. * * A 2-D DCT can be done by 1-D DCT on each row followed by 1-D DCT on each * column. * * Results: None * * Side effects: Overwrites the input data * * -------------------------------------------------------------- *//*__inline__ voidmp_fwd_dct_fast(data2d, dest2d)    Block data2d, dest2d;    */static inline void mp_fwd_dct_fast(Block data2d, Block dest2d){    int32 *data = (int32 *) data2d;	/* this algorithm wants					 * a 1-d array */    int32 *dest = (int32 *) dest2d;    int rowctr, columncounter;    register int32 *inptr, *outptr;    int32 workspace[DCTSIZE_SQ];    int32 tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7;    int32 tmp10, tmp11, tmp12, tmp13;    int32 tmp14, tmp15, tmp16, tmp17;    int32 tmp25, tmp26;    SHIFT_TEMPS;    /*     * Each iteration of the inner loop performs one 8-point 1-D DCT. It     * reads from a *row* of the input matrix and stores into a *column*     * of the output matrix.  In the first pass, we read from the data[]     * array and store into the local workspace[].  In the second pass,     * we read from the workspace[] array and store into data[], thus     * performing the equivalent of a columnar DCT pass with no variable     * array indexing.     */    inptr = data;		/* initialize pointers for first pass */    outptr = workspace;       /* PASS ONE */	    for (rowctr = DCTSIZE - 1; rowctr >= 0; rowctr--) {      /*       * many tmps have nonoverlapping lifetime -- flashy       * register colourers should be able to do this lot       * very well       */      /* SHIFT_TEMPS */            /* temp0 through tmp7:  -512 to +512 */      /* if I-block, then -256 to +256 */      tmp0 = inptr[7] + inptr[0];      tmp1 = inptr[6] + inptr[1];      tmp2 = inptr[5] + inptr[2];      tmp3 = inptr[4] + inptr[3];      tmp4 = inptr[3] - inptr[4];      tmp5 = inptr[2] - inptr[5];      tmp6 = inptr[1] - inptr[6];      tmp7 = inptr[0] - inptr[7];            /* tmp10 through tmp13:  -1024 to +1024 */      /* if I-block, then -512 to +512 */      tmp10 = tmp3 + tmp0;      tmp11 = tmp2 + tmp1;      tmp12 = tmp1 - tmp2;            tmp13 = tmp0 - tmp3;            outptr[0] = (int32) UNFIXH((tmp10 + tmp11) * SIN_1_4);      outptr[DCTSIZE * 4] = (int32) UNFIXH((tmp10 - tmp11) * COS_1_4);            outptr[DCTSIZE * 2] = (int32) UNFIXH(tmp13 * COS_1_8 + tmp12 * SIN_1_8);      outptr[DCTSIZE * 6] = (int32) UNFIXH(tmp13 * SIN_1_8 - tmp12 * COS_1_8);            tmp16 = UNFIXO((tmp6 + tmp5) * SIN_1_4);      tmp15 = UNFIXO((tmp6 - tmp5) * COS_1_4);            OVERSHIFT(tmp4);      OVERSHIFT(tmp7);            /*       * tmp4, tmp7, tmp15, tmp16 are overscaled by       * OVERSCALE       */            tmp14 = tmp4 + tmp15;      tmp25 = tmp4 - tmp15;      tmp26 = tmp7 - tmp16;      tmp17 = tmp7 + tmp16;            outptr[DCTSIZE] = (int32) UNFIXH(tmp17 * OCOS_1_16 + tmp14 * OSIN_1_16);      outptr[DCTSIZE * 7] = (int32) UNFIXH(tmp17 * OCOS_7_16 - tmp14 * OSIN_7_16);      outptr[DCTSIZE * 5] = (int32) UNFIXH(tmp26 * OCOS_5_16 + tmp25 * OSIN_5_16);      outptr[DCTSIZE * 3] = (int32) UNFIXH(tmp26 * OCOS_3_16 - tmp25 * OSIN_3_16);            inptr += DCTSIZE;	/* advance inptr to next row */      outptr++;		/* advance outptr to next column */        }    /* end of pass; in case it was pass 1, set up for pass 2 */    inptr = workspace;    outptr = dest;    columncounter = 0;    /* PASS TWO */    for (rowctr = DCTSIZE - 1; rowctr >= 0; rowctr--) {      /*       * many tmps have nonoverlapping lifetime -- flashy       * register colourers should be able to do this lot       * very well       */      /* SHIFT_TEMPS */            /* temp0 through tmp7:  -512 to +512 */      /* if I-block, then -256 to +256 */      tmp0 = inptr[7] + inptr[0];      tmp1 = inptr[6] + inptr[1];      tmp2 = inptr[5] + inptr[2];      tmp3 = inptr[4] + inptr[3];      tmp4 = inptr[3] - inptr[4];      tmp5 = inptr[2] - inptr[5];      tmp6 = inptr[1] - inptr[6];      tmp7 = inptr[0] - inptr[7];