/* idctref_miha.c, Inverse Discrete Fourier Transform, double precision */

/*************************************************************/
/*                                                           */
/* x87 hand-optimized assembly by Miha Peternel              */
/*                                     27.11. - 20.1.2001    */
/*                                                           */
/* You are free to use this code in your project if:         */
/* - no changes are made to this message                     */
/* - any changes to this code are publicly available         */
/* - your project documentation contains the following text: */
/*   "This software contains fast high-quality IDCT decoder  */
/*    by Miha Peternel."                                     */
/*                                                           */
/*************************************************************/

/*  Perform IEEE 1180 reference (64-bit floating point, separable 8x1
 *  direct matrix multiply) Inverse Discrete Cosine Transform
*/

/* Here we use math.h to generate constants.  Compiler results may vary a little */

#include <math.h>

#define M_PI	3.1415926535897932384626433832795
const static double HALF = 0.5;

/* cosine transform matrix for 8x1 IDCT */
static double c[8][8];

/* initialize DCT coefficient matrix */
void Initialize_REF_IDCT()
{
  int freq, time;
  double scale;

  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((M_PI/8.0)*freq*(time + 0.5));
  }
}

void REF_IDCT(short *block)
{
  double tmp[64];
	double rnd[64];
	int int0, int1, int2, int3, int4, int5, int6, int7;
	unsigned short fpold;
	unsigned short fpnew;

	int *b = (int *) block;


  if( !(b[0]|(b[31]&~0x10000)) )
	{
	  if( b[ 1]|b[ 2]|b[ 3]|b[ 4]|b[ 5]|b[ 6] )
		  goto normal;
	  if( b[ 7]|b[ 8]|b[ 9]|b[10]|b[11]|b[12] )
		  goto normal;
	  if( b[13]|b[14]|b[15]|b[16]|b[17]|b[18] )
		  goto normal;
	  if( b[19]|b[20]|b[21]|b[22]|b[23]|b[24] )
		  goto normal;
	  if( b[25]|b[26]|b[27]|b[28]|b[29]|b[30] )
		  goto normal;
		b[31]=0;
		return;
	}
normal:

	__asm
	{
		// do the IDCT
		mov esi,[block]
		lea eax,[c]
		lea edi,[tmp]
		//mov ebx,8
		mov ebx,8 // 0x77000000 // 8
		align 16
	__col1:
	    movzx edx,[esi+1*2]
			mov   ecx,[esi+2*2]
			or    edx,[esi+4*2]
			or    ecx,[esi+6*2]
			or edx,ecx
			//mov ecx,8
			mov ecx,8/2 // 0x77000000 // 8

			jnz __row1
				fild  word ptr [esi+0*2]
				fmul qword ptr [eax+0*8*8]
				fst  qword ptr [edi+0*8]
				fst  qword ptr [edi+1*8]
				fst  qword ptr [edi+2*8]
				fst  qword ptr [edi+3*8]
				fst  qword ptr [edi+4*8]
				fst  qword ptr [edi+5*8]
				fst  qword ptr [edi+6*8]
				fstp qword ptr [edi+7*8]
				add edi,8*8
				jmp __next1
			align 16
		__row1:
				fild  word ptr [esi+0*2]
				fmul qword ptr [eax+0*8*8]
				fild  word ptr [esi+1*2]
				fmul qword ptr [eax+1*8*8]
				fadd
				fild  word ptr [esi+2*2]
				fmul qword ptr [eax+2*8*8]
				fadd
				fild  word ptr [esi+3*2]
				fmul qword ptr [eax+3*8*8]
				fadd
				fild  word ptr [esi+4*2]
				fmul qword ptr [eax+4*8*8]
				fadd
				fild  word ptr [esi+5*2]
				fmul qword ptr [eax+5*8*8]
				fadd
				fild  word ptr [esi+6*2]
				fmul qword ptr [eax+6*8*8]
				fadd
				fild  word ptr [esi+7*2]
				fmul qword ptr [eax+7*8*8]
				fadd

				fild  word ptr [esi+0*2]
				fmul qword ptr [eax+0*8*8+8]
				fild  word ptr [esi+1*2]
				fmul qword ptr [eax+1*8*8+8]
				fadd
				fild  word ptr [esi+2*2]
				fmul qword ptr [eax+2*8*8+8]
				fadd
				fild  word ptr [esi+3*2]
				fmul qword ptr [eax+3*8*8+8]
				fadd
				fild  word ptr [esi+4*2]
				fmul qword ptr [eax+4*8*8+8]
				fadd
				fild  word ptr [esi+5*2]
				fmul qword ptr [eax+5*8*8+8]
				fadd
				fild  word ptr [esi+6*2]
				fmul qword ptr [eax+6*8*8+8]
				fadd
				fild  word ptr [esi+7*2]
				fmul qword ptr [eax+7*8*8+8]
				fadd
				add eax,8*2
				fxch st(1)
				fstp qword ptr [edi]//
				fstp qword ptr [edi+8]
				add edi,8*2

			dec ecx

			jnz __row1
			add eax,-8*8
			  //align 16
		__next1:
			add esi,+8*2

		sub ebx,0x80000001 // add ebx,ebx 
		js  __col1
			//align 16
			test ebx,ebx // align jump &| redo flags
		jnz __col1

		lea esi,[tmp]
		lea eax,[c]
		lea edi,[rnd]
		//mov edi,[block]
    fld qword ptr [HALF]
		mov ebx,8
	__row2:
			mov ecx,8/2
			align 16
			__col2:
				fld  qword ptr [esi+0*8*8]
				fmul qword ptr [eax+0*8*8]
				fld  qword ptr [esi+1*8*8]
				fmul qword ptr [eax+1*8*8]
				fadd
				fld  qword ptr [esi+2*8*8]
				fmul qword ptr [eax+2*8*8]
				fadd
				fld  qword ptr [esi+3*8*8]
				fmul qword ptr [eax+3*8*8]
				fadd
				fld  qword ptr [esi+4*8*8]
				fmul qword ptr [eax+4*8*8]
				fadd
				fld  qword ptr [esi+5*8*8]
				fmul qword ptr [eax+5*8*8]
				fadd
				fld  qword ptr [esi+6*8*8]
				fmul qword ptr [eax+6*8*8]
				fadd
				fld  qword ptr [esi+7*8*8]
				fmul qword ptr [eax+7*8*8]
				fadd
				fadd st(0),st(1)

				fxch st(1)

				fld  qword ptr [esi+0*8*8]
				fmul qword ptr [eax+0*8*8+8]
				fld  qword ptr [esi+1*8*8]
				fmul qword ptr [eax+1*8*8+8]
				fadd
				fld  qword ptr [esi+2*8*8]
				fmul qword ptr [eax+2*8*8+8]
				fadd
				fld  qword ptr [esi+3*8*8]
				fmul qword ptr [eax+3*8*8+8]
				fadd
				fld  qword ptr [esi+4*8*8]
				fmul qword ptr [eax+4*8*8+8]
				fadd
				fld  qword ptr [esi+5*8*8]
				fmul qword ptr [eax+5*8*8+8]
				fadd
				fld  qword ptr [esi+6*8*8]
				fmul qword ptr [eax+6*8*8+8]
				fadd
				fld  qword ptr [esi+7*8*8]
				fmul qword ptr [eax+7*8*8+8]
				fadd
				fadd st(0),st(1)
				add eax,8*2

				fxch st(2)
				fstp qword ptr [edi]
				fxch st(1)
				fstp qword ptr [edi+8*8]
				add edi,8*8*2

			dec ecx

			jnz __col2
			add eax,-8*8
			add esi,+8
			add edi,8-8*8*8

		sub ebx,0x80000001
		js  __row2
			  //align 16
				test ebx,ebx // align jump &| redo flags
		jnz __row2
		ffree st(0) // bye bye 0.5

	  // set x87 to floor mode
		fstcw [fpold]
		movzx eax, [fpold]

		or eax, 0x0400 // round down - floor
		mov [fpnew], ax
		fldcw [fpnew]

		// now floor the damn array
		lea esi, [rnd]
		mov edi, [block]
		mov ebx, -256 // clip min
		mov edx, +255 // clip max
		mov ecx, 8
		align 16
	__floor:
		  fld   qword ptr [esi+0*8]
			fistp dword ptr [int0]
			  mov eax,[int0]
				cmp   eax,ebx
				cmovl eax,ebx
				cmp   eax,edx
				cmovg eax,edx
		  fld   qword ptr [esi+1*8]
			fistp dword ptr [int1]
				mov word ptr [edi+0*2],ax
			  mov eax,[int1]
				cmp   eax,ebx
				cmovl eax,ebx
				cmp   eax,edx
				cmovg eax,edx
		  fld   qword ptr [esi+2*8]
			fistp dword ptr [int2]
				mov word ptr [edi+1*2],ax
			  mov eax,[int2]
				cmp   eax,ebx
				cmovl eax,ebx
				cmp   eax,edx
				cmovg eax,edx
		  fld   qword ptr [esi+3*8]
			fistp dword ptr [int3]
				mov word ptr [edi+2*2],ax
			  mov eax,[int3]
				cmp   eax,ebx
				cmovl eax,ebx
				cmp   eax,edx
				cmovg eax,edx
		  fld   qword ptr [esi+4*8]
			fistp dword ptr [int4]
				mov word ptr [edi+3*2],ax
			  mov eax,[int4]
				cmp   eax,ebx
				cmovl eax,ebx
				cmp   eax,edx
				cmovg eax,edx
		  fld   qword ptr [esi+5*8]
			fistp dword ptr [int5]
				mov word ptr [edi+4*2],ax
			  mov eax,[int5]
				cmp   eax,ebx
				cmovl eax,ebx
				cmp   eax,edx
				cmovg eax,edx
		  fld   qword ptr [esi+6*8]
			fistp dword ptr [int6]
				mov word ptr [edi+5*2],ax
			  mov eax,[int6]
				cmp   eax,ebx
				cmovl eax,ebx
				cmp   eax,edx
				cmovg eax,edx
		  fld   qword ptr [esi+7*8]
			fistp dword ptr [int7]
				mov word ptr [edi+6*2],ax
			  mov eax,[int7]
				cmp   eax,ebx
				cmovl eax,ebx
				cmp   eax,edx
				cmovg eax,edx
				mov word ptr [edi+7*2],ax

			add esi, 8*8
			add edi, 8*2

		sub ecx,0x80000001
		js  __floor
			  //align 16
				test ecx,ecx // align jump &| redo flags
		jnz __floor

		// set x87 to default mode
		fldcw [fpold]
	};

}
