/* This Digital Cosine Transform (DCT) algorithm implements a 2-dimensional */
/* 64-element transform.                                                    */
/* It is taken from T. Kientzle, "Implementing Fast DCTs",                  */
/* Dr Dobb's Journal, #297, March 1999, p. 119.                             */
/* Translated into PL/I from C language by R. A. Vowels, 8 May 1999.        */
(FIXEDOVERFLOW):
DCT2DTEST:
	PROCEDURE (DCTBLOCK) OPTIONS (REORDER);
	DECLARE DCTBLOCK(*,*) FIXED BINARY (16);
	DECLARE PI   FIXED BINARY(31,25) STATIC INITIAL (3.1415928);
	DECLARE C1   FIXED BINARY(31,10) INITIAL (COS(PI/16)),
		S1   FIXED BINARY(31,10) INITIAL (SIN(PI/16)),
		C3   FIXED BINARY(31,10) INITIAL (COS(3*PI/16)),
		S3   FIXED BINARY(31,10) INITIAL (SIN(3*PI/16)),
		R2C6 FIXED BINARY(31,10) INITIAL ((SQRT(2.00000)*COS(6*PI/16)) ),
		R2S6 FIXED BINARY(31,10) INITIAL ((SQRT(2.00000)*SIN(6*PI/16)) ),
		R2   FIXED BINARY(31, 7) STATIC INITIAL (SQRT(2.000000));
	DECLARE (Row, Col) FIXED BINARY(7);
	DECLARE (X0, X1, X2, X3, X4, X5, X6, X7) FIXED BINARY;
	DECLARE (Y0, Y1, Y2, Y3, Y4, Y5, Y6, Y7, Y8) FIXED BINARY (31,10);
	DECLARE X8 FIXED BINARY;
	DECLARE Z FIXED BINARY (31,3);

	DO Row = 0 TO 7;
		X0 = DCTBLOCK(Row,0); X1 = DCTBLOCK(Row,1); X2 = DCTBLOCK(Row,2);
		X3 = DCTBLOCK(Row,3); X4 = DCTBLOCK(Row,4); X5 = DCTBLOCK(Row,5);
		X6 = DCTBLOCK(Row,6); X7 = DCTBLOCK(Row,7);

		/* Stage 1 */
		X8 = X7 + X0;  X0 = X0 - X7;  X7 = X1 + X6;  X1 = X1 - X6;
		X6 = X2 + X5;  X2 = X2 - X5;  X5 = X3 + X4;  X3 = X3 - X4;

		/* Stage 2 */
		X4 = X8 + X5;  X8 = X8 - X5;  X5 = X7 + X6;  X7 = X7 - X6;
		Y6 = C1*(X1 + X2);  Y2 = (-S1 - C1)*X2 + Y6;  Y1 = (S1 - C1)*X1 + Y6;
		Y6 = C3*(X0 + X3);  Y3 = (-S3 - C3)*X3 + Y6;  Y0 = (S3 - C3)*X0 + Y6;

		/* Stage 3 */
		X6 = X4 + X5;  X4 = X4 - X5;
		Y5 = R2C6*(X7 + X8);  Y7 = (-R2S6 - R2C6)*X7 + Y5;  Y8 = (R2S6 - R2C6)*X8 + Y5;
		Y5 = Y0 + Y2;  Y0 = Y0 - Y2;  Y2 = Y3 + Y1;  Y3 = Y3 - Y1;

		/* Stage 4: Round and return result. */
		DCTBLOCK(Row,0) = X6;
		DCTBLOCK(Row,4) = X4;
		DCTBLOCK(Row,2) = ROUND(Y8, 0);
		DCTBLOCK(Row,6) = ROUND(Y7, 0);
		DCTBLOCK(Row,7) = ROUND(Y2 - Y5, 0);
		DCTBLOCK(Row,1) = ROUND(Y2 + Y5, 0);
		DCTBLOCK(Row,3) = ROUND(Y3*R2, 0);
		DCTBLOCK(Row,5) = ROUND(Y0*R2, 0);
	END;

	DO Col = 0 TO 7;
		X0 = DCTBLOCK(0,Col); X1 = DCTBLOCK(1,Col); X2 = DCTBLOCK(2,Col);
		X3 = DCTBLOCK(3,Col); X4 = DCTBLOCK(4,Col); X5 = DCTBLOCK(5,Col);
		X6 = DCTBLOCK(6,Col); X7 = DCTBLOCK(7,Col);

		/* Stage 1 */
		X8 = X7 + X0;  X0 = X0 - X7;  X7 = X1 + X6;  X1 = X1 - X6;
		X6 = X2 + X5;  X2 = X2 - X5;  X5 = X3 + X4;  X3 = X3 - X4;

		/* Stage 2 */
		X4 = X8 + X5;  X8 = X8 - X5;  X5 = X7 + X6;  X7 = X7 - X6;
		Y6 = C1*(X1 + X2);  Y2 = (-S1 - C1)*X2 + Y6;  Y1 = (S1 - C1)*X1 + Y6;
		Y6 = C3*(X0 + X3);  Y3 = (-S3 - C3)*X3 + Y6;  Y0 = (S3 - C3)*X0 + Y6;

		/* Stage 3 */
		X6 = X4 + X5;  X4 = X4 - X5;
		Y5 = R2C6*(X7 + X8);  Y7 = (-R2S6 - R2C6)*X7 + Y5;  Y8 = (R2S6 - R2C6)*X8 + Y5;
		Y5 = Y0 + Y2;  Y0 = Y0 - Y2;  Y2 = Y3 + Y1;  Y3 = Y3 - Y1;

		/* Stage 4: Round and return result. */

		DCTBLOCK(0,Col) = ROUND(X6, -4)/32;
		DCTBLOCK(4,Col) = ROUND(X4, -4)/32;
		DCTBLOCK(2,Col) = ROUND(Y8, -14)/32;
		DCTBLOCK(6,Col) = ROUND(Y7, -14)/32;
		DCTBLOCK(7,Col) = ROUND(Y2 - Y5, -14)/32;
		DCTBLOCK(1,Col) = ROUND(Y2 + Y5, -14)/32;
		Z = TRUNC(Y3/8192)*R2;
		DCTBLOCK(3,Col) = ROUND(Z*8, -2);
		Z = TRUNC(Y0/8192)*R2;
		DCTBLOCK(5,Col) = ROUND(Z*8, -2);
	END;


END  DCT2DTEST;