/**********************************************************************
 * 預處理定義
**********************************************************************/
#include <stdio.h>
#include <math.h>

/**********************************************************************
 * 函數聲明部分
**********************************************************************/
int ChangeMaxToTop(double **a,int row,int line,int m,int n);
int LieZhuYuan(double **a,int m,int n);
double QiuRoot(int num,double **a,int m,int n);
double Array(double **a,int i,int j,int n);
double DoubleAbs(double x);

/**********************************************************************
 * 函 數 名:getRoot
 * 功能描述:求根
 * 輸    入:a -- 待求根方程系數矩陣
 *			 m -- 二維數組行數
 *			 n -- 二維數組列數
 * 輸    出:x -- 根數組
 * 返    回:無
***********************************************************************/
void getRoot(double x[],double **a,int m,int n)
{
	int i;
	double **p;
	p = a;
	LieZhuYuan(a,m,n);

	for(i = 0;i < m;i ++)
	{
		x[i] = QiuRoot(i + 1,a,m,n);
		printf("x%d = %f\n",i + 1,x[i]);
	}
	/*
	x[1] = QiuRoot(2,p,3,4);
	printf("x2 = %f\n",x[1]);

	x[0] = QiuRoot(1,p,3,4);
	printf("x1 = %f\n",x[0]);
	*/
}
/**********************************************************************
 * 函 數 名:Array
 * 功能描述:計算在一維數組中的二維數據
 * 輸    入:a -- 一維數組 
 *			 i -- 二維數組行下標
 *           j -- 二維數組列下標
 *			 n -- 二維數組列數
 * 輸    出:無
 * 返    回:a[i][j]的值
**********************************************************************/
double Array(double **a,int i,int j,int n)
{
	double temp = *((double*)a + n * i + j);
	return temp;
}

/**********************************************************************
 * 函 數 名:ChangeMaxToTop
 * 功能描述:求第row行第line列中的最大的行，并與第row行交換位置
 * 輸    入:a -- 待交換數組
 *			 row -- 需要交換的數組行下標
 *           line -- 需要交換的數組列下標
 *			 m -- 二維數組行數
 *			 n -- 二維數組列數
 * 輸    出:交換后的數組
 * 返    回:1 -- 成功
 *			 0 -- 失敗
**********************************************************************/
int ChangeMaxToTop(double **a,int row,int line,int m,int n)
{
	int i,j;
	double temp;
	temp = Array(a,row,line,n);
	j = row;
	for(i = row + 1;i < m;i++) //找到絕對值比當前主元絕對值更大的數
		if(DoubleAbs(temp) < DoubleAbs(Array(a,i,line,n)))	
		{
			temp = Array(a,i,line,n);
			j = i;
		}
	if(j < row || j > m)
		return 0;
	else if(j == row)	//當前主元最大，不需要調換
		return 1;
	for(i = 0;i < n;i++)
	{
		temp = Array(a,j,i,n);
		*((double*)a + n * j + i) = Array(a,row,i,n);
		*((double*)a + n * row + i) = temp;
	}
	return 1;
}

/**********************************************************************
 * 函 數 名:LieZhuYuan
 * 功能描述:列主元法消元
 * 輸    入:a -- 系數矩陣
 *			 m -- 矩陣行數
 *			 n -- 矩陣列數
 * 輸    出:a -- 銷元后系數矩陣
 * 返    回:1 -- 成功
 *			 0 -- 失敗
**********************************************************************/
int LieZhuYuan(double **a,int m,int n)
{
	int i,j,k,flag;
	double temp;
	for(i = 0;i < m - 1;i++)
	{
		flag = ChangeMaxToTop(a,i,i,m,n);	//交換最大主元
		if (flag == 0)
			return 0;
		for(k = i + 1;k < m;k++)
		{
			temp = Array(a,k,i,n) / (0 - Array(a,i,i,n));
			for(j = 0;j < n;j++)	//置主元下面的系數為0
				*((double*)a + n * k + j) = Array(a,k,j,n) + Array(a,i,j,n) * temp;
		}
	}
	return 1;
}

/**********************************************************************
 * 函 數 名:QiuRoot
 * 功能描述:求根
 * 輸    入:num -- 要求的根
 *			 a -- 系數矩陣
 *			 m -- 矩陣行數
 *			 n -- 矩陣列數
 * 輸    出:無
 * 返    回:第num個根值
**********************************************************************/
double QiuRoot(int num,double **a,int m,int n)
{
	int j;
	double temp = 0;		//臨時變量，用來求和
	if(num == m)
		return Array(a,m - 1,n - 1,n) / Array(a,m - 1,m - 1,n);		//第m個根
	else
	{
		for(j = num;j < m;j++)
			temp = Array(a,num - 1,j,n) * QiuRoot(j + 1,a,m,n) + temp;	//遞歸調用
		return (Array(a,num - 1,n - 1,n) - temp) / Array(a,num - 1,num - 1,n);
	}
}

/**********************************************************************
 * 函 數 名:DoubleAbs
 * 功能描述:求double型數據的絕對值
 * 輸    入:x -- double型數據
 * 輸    出:無
 * 返    回:x的絕對值
**********************************************************************/
double DoubleAbs(double x)
{
	if (x < 0)
		return 0 - x;
	else
		return x;
}