horner계산을 이용한 테일러급수 계산 및 loop unrolling

void Eval_Poly_Naive(double y[], double x[], int n_x, double a[], int deg)
{
//	p = c[0] + c[1] * x + c[2] * pow(x, 2.0) + c[3] *
//		pow(x, 3.0) + c[4] * pow(x, 4.0);

	for (int j = 0; j < deg; j++) {
		y[j] = a[0];
		for (int i = 1; i < deg; i++) {
			y[j] += a[i] * pow(x[j], i);
		}

		printf("naive ret %lf\n", y[j]);
	}
}


void Eval_Poly_Horner(double y[], double x[], int n_x, double a[], int deg)
{
	for (int j = 0; j < deg; j++) {
		y[j] = a[deg - 1];  // Initialize result 

		// Evaluate value of polynomial using Horner's method 
		for (int i = deg - 2; i >= 0; --i) {
	//		printf("%lf %lf %lf\n", y[j], x[j], a[i]);
			y[j] = y[j] * x[j] + a[i];
		}

		printf("horner ret %lf\n", y[j]);
	}
}​

double Taylor_series_ex(double x, int n)
{
	//a5 x + a4 : a4는 영원히 1 ㅁ5는 1/n!
	double e = 1 + x/n;
	
	for (int i = n-1; i >= 1; --i) {
		//coeffiecient
		e = e * x/i + 1;
	}

	return e;
}


float Taylor_series(double x, int n)
{
	float e = 1 + x / n;

	for (int i = n - 1; i >= 1; --i) {
		//coeffiecient
		e = e * x / i + 1;
	}

	return e;
}

void MultiplySquareMatrices_1(double *pDestMatrix, double *pLeftMatrix ,double *pRightMatrix, int MatSize)
{
	int i, j, k;
	for (i = 0; i < MatSize; i++) {
		for (j = 0; j < MatSize; j++) {
			pDestMatrix[i * MatSize + j] = 0;
			for (k = 0; k < MatSize; k++) {
				pDestMatrix[i * MatSize + j] = pLeftMatrix[i * MatSize + k]
					* pRightMatrix[k * MatSize + j];
			}
		}
	}
}

void MultiplySquareMatrices_2(double *pDestMatrix, double *pLeftMatrix, double *pRightMatrix, int MatSize)
{
	int i, j, k;
	for (j = 0; j < MatSize; j++) {
		for (i = 0; i < MatSize; i++) {
			pDestMatrix[i * MatSize + j] = 0;
			for (k = 0; k < MatSize; k++) {
				pDestMatrix[i * MatSize + j] = pLeftMatrix[i * MatSize + k]
					* pRightMatrix[k * MatSize + j];
			}
		}
	}
}

void MultiplySquareMatrices_3(double *pDestMatrix, double *pLeftMatrix, double *pRightMatrix, int MatSize)
{
	int i, j, k;
	for (j = 0; j < MatSize; j++) {
		for (i = 0; i < MatSize; i++) {
			pDestMatrix[i * MatSize + j] = 0;

			//m개단위로 풀기
			for (k = 0; k < MatSize; k+=8) {
				pDestMatrix[i * MatSize + j] = pLeftMatrix[i * MatSize + k]
					* pRightMatrix[k * MatSize + j];
				pDestMatrix[i * MatSize + j] = pLeftMatrix[i * MatSize + (k + 1)]
					* pRightMatrix[(k + 1) * MatSize + j];
				pDestMatrix[i * MatSize + j] = pLeftMatrix[i * MatSize + (k + 2)]
					* pRightMatrix[(k + 2) * MatSize + j];
				pDestMatrix[i * MatSize + j] = pLeftMatrix[i * MatSize + (k + 3)]
					* pRightMatrix[(k + 3) * MatSize + j];
				pDestMatrix[i * MatSize + j] = pLeftMatrix[i * MatSize + (k + 4)]
					* pRightMatrix[(k + 4) * MatSize + j];
				pDestMatrix[i * MatSize + j] = pLeftMatrix[i * MatSize + (k + 5)]
					* pRightMatrix[(k + 5) * MatSize + j];
				pDestMatrix[i * MatSize + j] = pLeftMatrix[i * MatSize + (k + 6)]
					* pRightMatrix[(k + 6) * MatSize + j];
				pDestMatrix[i * MatSize + j] = pLeftMatrix[i * MatSize + (k + 7)]
					* pRightMatrix[(k + 7) * MatSize + j];
			}
		}
	}
}

void MultiplySquareMatrices_4(double *pDestMatrix, double *pLeftMatrix, double *pRightMatrix, int MatSize)
{
	int i, j, k;
	for (j = 0; j < MatSize; j++) {
		for (i = 0; i < MatSize; i++) {
			pDestMatrix[i * MatSize + j] = 0;

			//m개단위로 풀기
			for (k = 0; k < MatSize; k += 16) {
				pDestMatrix[i * MatSize + j] = pLeftMatrix[i * MatSize + k]
					* pRightMatrix[k * MatSize + j];
				pDestMatrix[i * MatSize + j] = pLeftMatrix[i * MatSize + (k + 1)]
					* pRightMatrix[(k + 1) * MatSize + j];
				pDestMatrix[i * MatSize + j] = pLeftMatrix[i * MatSize + (k + 2)]
					* pRightMatrix[(k + 2) * MatSize + j];
				pDestMatrix[i * MatSize + j] = pLeftMatrix[i * MatSize + (k + 3)]
					* pRightMatrix[(k + 3) * MatSize + j];
				pDestMatrix[i * MatSize + j] = pLeftMatrix[i * MatSize + (k + 4)]
					* pRightMatrix[(k + 4) * MatSize + j];
				pDestMatrix[i * MatSize + j] = pLeftMatrix[i * MatSize + (k + 5)]
					* pRightMatrix[(k + 5) * MatSize + j];
				pDestMatrix[i * MatSize + j] = pLeftMatrix[i * MatSize + (k + 6)]
					* pRightMatrix[(k + 6) * MatSize + j];
				pDestMatrix[i * MatSize + j] = pLeftMatrix[i * MatSize + (k + 7)]
					* pRightMatrix[(k + 7) * MatSize + j];
				pDestMatrix[i * MatSize + j] = pLeftMatrix[i * MatSize + (k + 8)]
					* pRightMatrix[(k + 8) * MatSize + j];
				pDestMatrix[i * MatSize + j] = pLeftMatrix[i * MatSize + (k + 9)]
					* pRightMatrix[(k + 9) * MatSize + j];
				pDestMatrix[i * MatSize + j] = pLeftMatrix[i * MatSize + (k + 10)]
					* pRightMatrix[(k + 10) * MatSize + j];
				pDestMatrix[i * MatSize + j] = pLeftMatrix[i * MatSize + (k + 11)]
					* pRightMatrix[(k + 11) * MatSize + j];
				pDestMatrix[i * MatSize + j] = pLeftMatrix[i * MatSize + (k + 12)]
					* pRightMatrix[(k + 12) * MatSize + j];
				pDestMatrix[i * MatSize + j] = pLeftMatrix[i * MatSize + (k + 13)]
					* pRightMatrix[(k + 13) * MatSize + j];
				pDestMatrix[i * MatSize + j] = pLeftMatrix[i * MatSize + (k + 14)]
					* pRightMatrix[(k + 14) * MatSize + j];
				pDestMatrix[i * MatSize + j] = pLeftMatrix[i * MatSize + (k + 15)]
					* pRightMatrix[(k + 15) * MatSize + j];
			}
		}
	}
}