数学计算

开发平台：
C++ Builder

idct.c：源码内容
							/* idct.c, inverse fast discrete cosine transform                           */
/* Copyright (C) 1996, MPEG Software Simulation Group. All Rights Reserved. */
/*
 * Disclaimer of Warranty
 *
 * These software programs are available to the user without any license fee or
 * royalty on an "as is" basis.  The MPEG Software Simulation Group disclaims
 * any and all warranties, whether express, implied, or statuary, including any
 * implied warranties or merchantability or of fitness for a particular
 * purpose.  In no event shall the copyright-holder be liable for any
 * incidental, punitive, or consequential damages of any kind whatsoever
 * arising from the use of these programs.
 *
 * This disclaimer of warranty extends to the user of these programs and user's
 * customers, employees, agents, transferees, successors, and assigns.
 *
 * The MPEG Software Simulation Group does not represent or warrant that the
 * programs furnished hereunder are free of infringement of any third-party
 * patents.
 *
 * Commercial implementations of MPEG-1 and MPEG-2 video, including shareware,
 * are subject to royalty fees to patent holders.  Many of these patents are
 * general enough such that they are unavoidable regardless of implementation
 * design.
 *
 */
/**********************************************************/
/* inverse two dimensional DCT, Chen-Wang algorithm       */
/* (cf. IEEE ASSP-32, pp. 803-816, Aug. 1984)             */
/* 32-bit integer arithmetic (8 bit coefficients)         */
/* 11 mults, 29 adds per DCT                              */
/*                                      sE, 18.8.91       */
/**********************************************************/
/* coefficients extended to 12 bit for IEEE1180-1990      */
/* compliance                           sE,  2.1.94       */
/**********************************************************/
/* floating point conversion by Miha Peternel             */
/*                                      27-29.11.2000     */
/**********************************************************/
//#define ModelX 123 // enable old code
//#define DoTest 123 // enable error logging to idct_log.txt
#ifdef DoTest
#include <stdio.h>
static FILE *logf; /* log file */
static int iter = 0;
static int isum = 0;
static int imax = 0;
static double fsum = 0;
static double fmax = 0;
static double c[8][8]; /* cosine transform matrix for 8x1 IDCT */
extern void idct_mmx_32( short *blk );
#include <math.h>
#ifndef PI
# ifdef M_PI
#  define PI M_PI
# else
#  define PI 3.14159265358979323846
# endif
#endif
#endif
#ifdef ModelX
/////////////////////////////////////////////////////
//
// old integer stuff is disabled... scroll down !!!
//
/////////////////////////////////////////////////////
/* this code assumes >> to be a two's-complement arithmetic */
/* right shift: (-2)>>1 == -1 , (-3)>>1 == -2               */
#include "config.h"
#define W1 2841 /* 2048*sqrt(2)*cos(1*pi/16) */
#define W2 2676 /* 2048*sqrt(2)*cos(2*pi/16) */
#define W3 2408 /* 2048*sqrt(2)*cos(3*pi/16) */
#define W5 1609 /* 2048*sqrt(2)*cos(5*pi/16) */
#define W6 1108 /* 2048*sqrt(2)*cos(6*pi/16) */
#define W7 565  /* 2048*sqrt(2)*cos(7*pi/16) */
/* global declarations */
void Initialize_Fast_IDCT _ANSI_ARGS_((void));
void Fast_IDCT _ANSI_ARGS_((short *block));
/* private data */
static short iclip[1024]; /* clipping table */
static short *iclp;
/* private prototypes */
static void idctrow _ANSI_ARGS_((short *blk));
static void idctcol _ANSI_ARGS_((short *blk));
/* row (horizontal) IDCT
 *
 *           7                       pi         1
 * dst[k] = sum c[l] * src[l] * cos( -- * ( k + - ) * l )
 *          l=0                      8          2
 *
 * where: c[0]    = 128
 *        c[1..7] = 128*sqrt(2)
 */
static void idctrow(blk)
short *blk;
{
  int x0, x1, x2, x3, x4, x5, x6, x7, x8;
  /* shortcut */
  if (!((x1 = blk[4]<<11) | (x2 = blk[6]) | (x3 = blk[2]) |
        (x4 = blk[1]) | (x5 = blk[7]) | (x6 = blk[5]) | (x7 = blk[3])))
  {
    blk[0]=blk[1]=blk[2]=blk[3]=blk[4]=blk[5]=blk[6]=blk[7]=blk[0]<<3;
    return;
  }
  x0 = (blk[0]<<11) + 128; /* for proper rounding in the fourth stage */
  /* first stage */
  x8 = W7*(x4+x5);
  x4 = x8 + (W1-W7)*x4;
  x5 = x8 - (W1+W7)*x5;
  x8 = W3*(x6+x7);
  x6 = x8 - (W3-W5)*x6;
  x7 = x8 - (W3+W5)*x7;
  
  /* second stage */
  x8 = x0 + x1;
  x0 -= x1;
  x1 = W6*(x3+x2);
  x2 = x1 - (W2+W6)*x2;
  x3 = x1 + (W2-W6)*x3;
  x1 = x4 + x6;
  x4 -= x6;
  x6 = x5 + x7;
  x5 -= x7;
  
  /* third stage */
  x7 = x8 + x3;
  x8 -= x3;
  x3 = x0 + x2;
  x0 -= x2;
  x2 = (181*(x4+x5)+128)>>8;
  x4 = (181*(x4-x5)+128)>>8;
  
  /* fourth stage */
  blk[0] = (x7+x1)>>8;
  blk[1] = (x3+x2)>>8;
  blk[2] = (x0+x4)>>8;
  blk[3] = (x8+x6)>>8;
  blk[4] = (x8-x6)>>8;
  blk[5] = (x0-x4)>>8;
  blk[6] = (x3-x2)>>8;
  blk[7] = (x7-x1)>>8;
}
/* column (vertical) IDCT
 *
 *             7                         pi         1
 * dst[8*k] = sum c[l] * src[8*l] * cos( -- * ( k + - ) * l )
 *            l=0                        8          2
 *
 * where: c[0]    = 1/1024
 *        c[1..7] = (1/1024)*sqrt(2)
 */
static void idctcol(blk)
short *blk;
{
  int x0, x1, x2, x3, x4, x5, x6, x7, x8;
  /* shortcut */
  if (!((x1 = (blk[8*4]<<8)) | (x2 = blk[8*6]) | (x3 = blk[8*2]) |
        (x4 = blk[8*1]) | (x5 = blk[8*7]) | (x6 = blk[8*5]) | (x7 = blk[8*3])))
  {
    blk[8*0]=blk[8*1]=blk[8*2]=blk[8*3]=blk[8*4]=blk[8*5]=blk[8*6]=blk[8*7]=
      iclp[(blk[8*0]+32)>>6];
    return;
  }
  x0 = (blk[8*0]<<8) + 8192;
  /* first stage */
  x8 = W7*(x4+x5) + 4;
  x4 = (x8+(W1-W7)*x4)>>3;
  x5 = (x8-(W1+W7)*x5)>>3;
  x8 = W3*(x6+x7) + 4;
  x6 = (x8-(W3-W5)*x6)>>3;
  x7 = (x8-(W3+W5)*x7)>>3;
  
  /* second stage */
  x8 = x0 + x1;
  x0 -= x1;
  x1 = W6*(x3+x2) + 4;
  x2 = (x1-(W2+W6)*x2)>>3;
  x3 = (x1+(W2-W6)*x3)>>3;
  x1 = x4 + x6;
  x4 -= x6;
  x6 = x5 + x7;
  x5 -= x7;
  
  /* third stage */
  x7 = x8 + x3;
  x8 -= x3;
  x3 = x0 + x2;
  x0 -= x2;
  x2 = (181*(x4+x5)+128)>>8;
  x4 = (181*(x4-x5)+128)>>8;
  
  /* fourth stage */
  blk[8*0] = iclp[(x7+x1)>>14];
  blk[8*1] = iclp[(x3+x2)>>14];
  blk[8*2] = iclp[(x0+x4)>>14];
  blk[8*3] = iclp[(x8+x6)>>14];
  blk[8*4] = iclp[(x8-x6)>>14];
  blk[8*5] = iclp[(x0-x4)>>14];
  blk[8*6] = iclp[(x3-x2)>>14];
  blk[8*7] = iclp[(x7-x1)>>14];
}
/* two dimensional inverse discrete cosine transform */
void Fast_IDCT(block)
short *block;
{
  int i;
#ifdef DoTest
  int j, k, v;
  double partial_product;
  double tmp[64];
	double hfsum = 0;
	double hfmax = 0;
	int hisum = 0;
	int himax = 0;
  for (i=0; i<8; i++)
    for (j=0; j<8; j++)
    {
      partial_product = 0.0;
      for (k=0; k<8; k++)
        partial_product+= c[k][j]*block[8*i+k];
      tmp[8*i+j] = partial_product;
    }
#endif
  for (i=0; i<8; i++)
    idctrow(block+8*i);
  for (i=0; i<8; i++)
    idctcol(block+i);
#ifdef DoTest
  /* Transpose operation is integrated into address mapping by switching 
     loop order of i and j */
  for (j=0; j<8; j++)
    for (i=0; i<8; i++)
    {
		  int ierr;
			double ferr;
      partial_product = 0.0;
      for (k=0; k<8; k++)
        partial_product+= c[k][i]*tmp[8*k+j];
      v = (int) floor(partial_product+0.5);
      v = (v<-256) ? -256 : ((v>255) ? 255 : v);
      ierr = abs(v-block[8*i+j]);
			//ferr = fabs(partial_product-fblock[8*i+j]);
			ferr = fabs(partial_product-block[8*i+j]); // integer test only
			block[8*i+j] = v;
			hisum += ierr; if( ierr > himax ) himax = ierr;
			hfsum += ferr; if( ferr > hfmax ) hfmax = ferr;
    }
	iter++;
	isum += hisum;
	if( himax > imax ) imax = himax;
	fsum += hfsum;
	if( hfmax > fmax ) fmax = hfmax;
	if( himax > 0 || (iter & ((1<<10)-1)) == 0 )
	fprintf( logf, "%6d | %lg/%d %lg/%lg | %lg/%d %lg/%lg | %d %lgn",
	  iter,
	  hisum/64.0, himax,
		hfsum/64.0, hfmax,
	  isum/(64.0*iter), imax,
		fsum/(64.0*iter), fmax,
		isum, fsum );
#endif
}
void Initialize_Fast_IDCT()
{
  int i;
#ifdef DoTest
  int freq, time;
  double scale;
#endif
  iclp = iclip+512;
  for (i= -512; i<512; i++)
    iclp[i] = (i<-256) ? -256 : ((i>255) ? 255 : i);
#ifdef DoTest
  for (freq=0; freq < 8; freq++)
  {
    scale = (freq == 0) ? sqrt(0.125) : 0.5;
    for (time=0; time<8; time++)
      c[freq][time] = scale*cos((PI/8.0)*freq*(time + 0.5));
  }
	logf = fopen( "idct_log.txt", "wt" );
	// warning: we'll never close the file !!!!
	iter = 0;
	isum = 0;
	imax = 0;
	fsum = 0;
	fmax = 0;
#endif
}
//////////////////////////////////////////////////////////////
#else // ModelX
/////////////////////////////////////////////////////
//
// here comes the real thing - 64-bit floating point
// converted from MPEG sources by Miha Peternel
//
// TODO:
// - loops can be easily vectorized for SIMD
//
/////////////////////////////////////////////////////
#include "config.h"
#include <math.h>
#ifndef PI
# ifdef M_PI
#  define PI M_PI
# else
#  define PI 3.14159265358979323846
# endif
#endif
#define FLOAT double
///*
#define SCALE 1.0 // 2048.0
#define FACT1 1.0
#define FACT2 (1.0/8) // 256.0
/*
#define SCALE 2048.0
#define FACT1 16384.0
#define FACT2 256.0
/**/
const static double RC = 1.0*1024*1024*1024*1024*256*16 + 1024; // magic + clip center
static FLOAT W1; // 2841 /* 2048*sqrt(2)*cos(1*pi/16) */
static FLOAT W2; // 2676 /* 2048*sqrt(2)*cos(2*pi/16) */
static FLOAT W5; // 1609 /* 2048*sqrt(2)*cos(5*pi/16) */
static FLOAT W1_8;
static FLOAT W2_8;
static FLOAT W5_8;
static FLOAT W7; // 565  /* 2048*sqrt(2)*cos(7*pi/16) */
static FLOAT W1mW7; // W1-W7
static FLOAT W1pW7; // W1+W7
static FLOAT W3; // 2408 /* 2048*sqrt(2)*cos(3*pi/16) */
static FLOAT W3mW5; // W3-W5
static FLOAT W3pW5; // W3+W5
static FLOAT W6; // 1108 /* 2048*sqrt(2)*cos(6*pi/16) */
static FLOAT W2mW6; // W2-W6
static FLOAT W2pW6; // W2+W6
static FLOAT S2; // 1/sqrt(2)
static FLOAT D8 = FACT2; // 1.0/8
static FLOAT W7_8;
static FLOAT W1mW7_8;
static FLOAT W1pW7_8;
static FLOAT W3_8;
static FLOAT W3mW5_8;
static FLOAT W3pW5_8;
static FLOAT W6_8;
static FLOAT W2mW6_8;
static FLOAT W2pW6_8;
//static FLOAT fblock[8*8];
/* global declarations */
void Initialize_Fast_IDCT _ANSI_ARGS_((void));
void Fast_IDCT _ANSI_ARGS_((short *block));
/* private data */
static short iclip[1024+1024]; /* clipping table */
static short *iclp;
/* private prototypes */
//static void idctrow _ANSI_ARGS_((short *blk));
//static void idctcol _ANSI_ARGS_((short *blk));
/* row (horizontal) IDCT
 *
 *           7                       pi         1
 * dst[k] = sum c[l] * src[l] * cos( -- * ( k + - ) * l )
 *          l=0                      8          2
 *
 * where: c[0]    = 128
 *        c[1..7] = 128*sqrt(2)
 */
static void idctrow( short *blk, FLOAT *fblk )
{
  FLOAT x0, x1, x2, x3, x4, x5, x6, x7, x8;
  /* shortcut */
///*
	int *i = (int *)blk;
  if( !( blk[1] | i[1] | i[2] | i[3] ))
  {
    fblk[0]=fblk[1]=fblk[2]=fblk[3]=fblk[4]=fblk[5]=fblk[6]=fblk[7]= blk[0]*SCALE;//*8.0; //8=SCALE/FACT2
    return;
  }
//*/
  /* shortcut */
	/*
  if (!((x1 = blk[4]<<11) | (x2 = blk[6]) | (x3 = blk[2]) |
        (x4 = blk[1]) | (x5 = blk[7]) | (x6 = blk[5]) | (x7 = blk[3])))
  {
    blk[0]=blk[1]=blk[2]=blk[3]=blk[4]=blk[5]=blk[6]=blk[7]=blk[0]<<3;
    return;
  }
	*/
	x1 = blk[4]*SCALE;//*2048;
	x2 = blk[6];
	x3 = blk[2];
  x4 = blk[1];
	x5 = blk[7];
	x6 = blk[5];
	x7 = blk[3];
  //x0 = (blk[0]*SCALE) + 128; /* for proper rounding in the fourth stage */
  x0 = blk[0]*SCALE;//*2048;
  /* first stage */
  x8 = W7*(x4+x5);
  x4 = x8 + (W1mW7)*x4;
  x5 = x8 - (W1pW7)*x5;
  x8 = W3*(x6+x7);
  x6 = x8 - (W3mW5)*x6;
  x7 = x8 - (W3pW5)*x7;
  
  /* second stage */
  x8 = x0 + x1;
  x0 -= x1;
  x1 = W6*(x3+x2);
  x2 = x1 - (W2pW6)*x2;
  x3 = x1 + (W2mW6)*x3;
  x1 = x4 + x6;
  x4 -= x6;
  x6 = x5 + x7;
  x5 -= x7;
  
  /* third stage */
  x7 = x8 + x3;
  x8 -= x3;
  x3 = x0 + x2;
  x0 -= x2;
  x2 = S2*(x4+x5);
  x4 = S2*(x4-x5);
  
  /* fourth stage */
  fblk[0] = (x7+x1); //*(1.0/256);
  fblk[7] = (x7-x1); //*(1.0/256);
  fblk[3] = (x8+x6); //*(1.0/256);
  fblk[4] = (x8-x6); //*(1.0/256);
  fblk[1] = (x3+x2); //*(1.0/256);
  fblk[6] = (x3-x2); //*(1.0/256);
  fblk[2] = (x0+x4); //*(1.0/256);
  fblk[5] = (x0-x4); //*(1.0/256);
}
/* column (vertical) IDCT
 *
 *             7                         pi         1
 * dst[8*k] = sum c[l] * src[8*l] * cos( -- * ( k + - ) * l )
 *            l=0                        8          2
 *
 * where: c[0]    = 1/1024
 *        c[1..7] = (1/1024)*sqrt(2)
 */
static void idctcol( short *blk, FLOAT *fblk )
{
  double rc[8];
  FLOAT x0, x1, x2, x3, x4, x5, x6, x7, x8;
  /* shortcut */
/*
  if (!((x1 = (blk[8*4]<<8)) | (x2 = blk[8*6]) | (x3 = blk[8*2]) |
        (x4 = blk[8*1]) | (x5 = blk[8*7]) | (x6 = blk[8*5]) | (x7 = blk[8*3])))
  {
    blk[8*0]=blk[8*1]=blk[8*2]=blk[8*3]=blk[8*4]=blk[8*5]=blk[8*6]=blk[8*7]=
      iclp[(blk[8*0]+32)>>6];
    return;
  }
*/
	x1 = fblk[8*4]*D8;
	x2 = fblk[8*6];
	x3 = fblk[8*2];
	x4 = fblk[8*1];
	x5 = fblk[8*7];
	x6 = fblk[8*5];
	x7 = fblk[8*3];
  x0 = fblk[8*0]*D8;
  /* first stage */
  x8 = W7_8*(x4+x5);
  x4 = (x8+(W1mW7_8)*x4);
  x5 = (x8-(W1pW7_8)*x5);
  x8 = W3_8*(x6+x7);
  x6 = (x8-(W3mW5_8)*x6);
  x7 = (x8-(W3pW5_8)*x7);
  
  /* second stage */
  x8 = x0 + x1;
  x0 -= x1;
  x1 = W6_8*(x3+x2);
  x2 = (x1-(W2pW6_8)*x2);
  x3 = (x1+(W2mW6_8)*x3);
  x1 = x4 + x6;
  x4 -= x6;
  x6 = x5 + x7;
  x5 -= x7;
  
  /* third stage */
  x7 = x8 + x3;
  x8 -= x3;
  x3 = x0 + x2;
  x0 -= x2;
  x2 = S2*(x4+x5);
  x4 = S2*(x4-x5);
  
  /* fourth stage */
  //blk[8*0] = iclp[(int)((x7+x1)*(1.0/FACT1/FACT2)+0.5)];//*(1.0/16384/256)];
  //blk[8*7] = iclp[(int)((x7-x1)*(1.0/FACT1/FACT2)+0.5)];//*(1.0/16384/256)];
  //blk[8*3] = iclp[(int)((x8+x6)*(1.0/FACT1/FACT2)+0.5)];//*(1.0/16384/256)];
  //blk[8*4] = iclp[(int)((x8-x6)*(1.0/FACT1/FACT2)+0.5)];//*(1.0/16384/256)];
  //blk[8*1] = iclp[(int)((x3+x2)*(1.0/FACT1/FACT2)+0.5)];//*(1.0/16384/256)];
  //blk[8*6] = iclp[(int)((x3-x2)*(1.0/FACT1/FACT2)+0.5)];//*(1.0/16384/256)];
  //blk[8*2] = iclp[(int)((x0+x4)*(1.0/FACT1/FACT2)+0.5)];//*(1.0/16384/256)];
  //blk[8*5] = iclp[(int)((x0-x4)*(1.0/FACT1/FACT2)+0.5)];//*(1.0/16384/256)];
#ifdef DoTest
  fblk[8*0] = (x7+x1);//*(1.0/16384/256)];
  fblk[8*7] = (x7-x1);//*(1.0/16384/256)];
  fblk[8*3] = (x8+x6);//*(1.0/16384/256)];
  fblk[8*4] = (x8-x6);//*(1.0/16384/256)];
  fblk[8*1] = (x3+x2);//*(1.0/16384/256)];
  fblk[8*6] = (x3-x2);//*(1.0/16384/256)];
  fblk[8*2] = (x0+x4);//*(1.0/16384/256)];
  fblk[8*5] = (x0-x4);//*(1.0/16384/256)];
#endif
  rc[0] = (x7+x1)+RC;
  rc[7] = (x7-x1)+RC;
  rc[3] = (x8+x6)+RC;
  rc[4] = (x8-x6)+RC;
  rc[1] = (x3+x2)+RC;
  rc[6] = (x3-x2)+RC;
  rc[2] = (x0+x4)+RC;
  rc[5] = (x0-x4)+RC;
	blk[8*0] = iclip[*((int*)(&rc[0]))];
	blk[8*7] = iclip[*((int*)(&rc[7]))];
	blk[8*3] = iclip[*((int*)(&rc[3]))];
	blk[8*4] = iclip[*((int*)(&rc[4]))];
	blk[8*1] = iclip[*((int*)(&rc[1]))];
	blk[8*6] = iclip[*((int*)(&rc[6]))];
	blk[8*2] = iclip[*((int*)(&rc[2]))];
	blk[8*5] = iclip[*((int*)(&rc[5]))];
}
/* two dimensional inverse discrete cosine transform */
static FLOAT fblock[8*8];
void Fast_IDCT(block)
short *block;
{
  int i;
  //FLOAT fblock[8*8];
#ifdef DoTest
  int j, k, v;
  double partial_product;
  double tmp[64];
	double hfsum = 0;
	double hfmax = 0;
	int hisum = 0;
	int himax = 0;
  for (i=0; i<8; i++)
    for (j=0; j<8; j++)
    {
      partial_product = 0.0;
      for (k=0; k<8; k++)
        partial_product+= c[k][j]*block[8*i+k];
      tmp[8*i+j] = partial_product;
    }
#endif
  for (i=0; i<8; i++)
    idctrow(block+8*i,fblock+8*i);
  for (i=0; i<8; i++)
    idctcol(block+i,fblock+i);
  //idct_mmx_32( block ); // test only
#ifdef DoTest
  /* Transpose operation is integrated into address mapping by switching 
     loop order of i and j */
  for (j=0; j<8; j++)
    for (i=0; i<8; i++)
    {
		  int ierr;
			double ferr;
      partial_product = 0.0;
      for (k=0; k<8; k++)
        partial_product+= c[k][i]*tmp[8*k+j];
      v = (int) floor(partial_product+0.5);
      v = (v<-256) ? -256 : ((v>255) ? 255 : v);
      ierr = abs(v-block[8*i+j]);
			ferr = fabs(partial_product-fblock[8*i+j]);
			//ferr = fabs(partial_product-block[8*i+j]); // integer test only
			block[8*i+j] = v;
			hisum += ierr; if( ierr > himax ) himax = ierr;
			hfsum += ferr; if( ferr > hfmax ) hfmax = ferr;
    }
	iter++;
	isum += hisum;
	if( himax > imax ) imax = himax;
	fsum += hfsum;
	if( hfmax > fmax ) fmax = hfmax;
	if( himax > 0 || (iter & ((1<<10)-1)) == 0 )
	fprintf( logf, "%6d | %lg/%d %lg/%lg | %lg/%d %lg/%lg | %d %lgn",
	  iter,
	  hisum/64.0, himax,
		hfsum/64.0, hfmax,
	  isum/(64.0*iter), imax,
		fsum/(64.0*iter), fmax,
		isum, fsum );
#endif
}
void Initialize_Fast_IDCT()
{
  int i;
#ifdef DoTest
  int freq, time;
  double scale;
#endif
  S2 = sqrt(0.5); // 1.0/sqrt(2);
  W1 = SCALE*sqrt(2)*cos(PI*(1.0/16)); 
	W1_8 = W1/8;
  W2 = SCALE*sqrt(2)*cos(PI*(2.0/16)); 
	W2_8 = W2/8;
  W3 = SCALE*sqrt(2)*cos(PI*(3.0/16)); 
	W3_8 = W3/8;
  W5 = SCALE*sqrt(2)*cos(PI*(5.0/16)); 
	W5_8 = W5/8;
  W6 = SCALE*sqrt(2)*cos(PI*(6.0/16)); 
	W6_8 = W6/8;
  W7 = SCALE*sqrt(2)*cos(PI*(7.0/16));
	W7_8 = W7/8;
  //W1 = 16*sqrt(2)*cos(PI*(1.0/16)); 
  //W2 = 16*sqrt(2)*cos(PI*(2.0/16)); 
  //W3 = 16*sqrt(2)*cos(PI*(3.0/16)); 
  //W5 = 16*sqrt(2)*cos(PI*(5.0/16)); 
  //W6 = 16*sqrt(2)*cos(PI*(6.0/16)); 
  //W7 = 16*sqrt(2)*cos(PI*(7.0/16));
  W1mW7 = W1-W7;  W1mW7_8 = W1mW7/8;
  W1pW7 = W1+W7;  W1pW7_8 = W1pW7/8;
  W3mW5 = W3-W5;  W3mW5_8 = W3mW5/8;
  W3pW5 = W3+W5;  W3pW5_8 = W3pW5/8;
  W2mW6 = W2-W6;  W2mW6_8 = W2mW6/8;
  W2pW6 = W2+W6;  W2pW6_8 = W2pW6/8;
  iclp = iclip+1024;
  for (i= -1024; i<1024; i++)
    iclp[i] = (i<-256) ? -256 : ((i>255) ? 255 : i);
#ifdef DoTest
  for (freq=0; freq < 8; freq++)
  {
    scale = (freq == 0) ? sqrt(0.125) : 0.5;
    for (time=0; time<8; time++)
      c[freq][time] = scale*cos((PI/8.0)*freq*(time + 0.5));
  }
	logf = fopen( "idct_log.txt", "wt" );
	// warning: we'll never close the file !!!!
	iter = 0;
	isum = 0;
	imax = 0;
	fsum = 0;
	fmax = 0;
#endif
}
#endif // ModelX

						