header.h

求解一元高次方程

--------------------------------------------------------------------------------
作者     : <xxx>
</PRE>
*******************************************************************************/
void Sturm(double* dp_equation,int i_degree,double** &dp_sturm,int &i_count)
{
	i_count=2;
	int i = 0;

	for (i=0;i<=i_degree;i++)
	{
		dp_sturm[0][i]=dp_equation[i];
		if (0==i)
		{
			dp_sturm[1][0]=0.0;
		}
		else
		{
			dp_sturm[1][i]=(i_degree-i+1)*dp_equation[i-1];//次数对齐了
		}
	}

	ArrayDisplay(dp_sturm[0],i_degree);
	ArrayDisplay(dp_sturm[1],i_degree);


// 	int flag=i_degree;
// 	int i_lowerDegree; 
	while (true)
	{
		dp_sturm[i_count]=Divide(dp_sturm[i_count-2],dp_sturm[i_count-1],i_degree);

		int flag = GetDegree(dp_sturm[i_count],i_degree);

		if (-1 == flag)
		{
			return;
		}
// 		ArrayDisplay(dp_sturm[i_count],i_degree);
// 		i_lowerDegree	= i_degree;
// 		for (i = 0;i<=i_degree;i++)
// 		{
// 			if (dp_sturm[i_count][i]<=1e-6)				//??????????
// 			{
// 				i_lowerDegree--;
// 			}
// 		}                         /*注------------------------*/
// 		if (i_lowerDegree<=0)   ///1、对于重根来说，fs+1(x)(余式)必=0，即i_lowerDegree=-1停止，fs(x)必为上次的值，一次项系数！=0
// 		{                       ///2、对于单根来说，fs(x)(余式)为常数，即i_lowerDegree=0停止，但是此时余式恰好是sturm组的最后一项fs(x)
// 			break;              ///   如果这时停止，返回的i_count是上一次的值，这时sturm序列是否少了一个fs(x)?
// 			                    ///3、所以令i_lowerDegree<0时停止，这时让重根和单根其最终余式为0时停止，这时有两种情况
// 		}                              ///1、如果上一次余式是非0次多项式(i_lowerDegree>=1),  则有重根 
//                                        ///2、如果上一次余式是0次多项式(i_lowerDegree=0)，常数项，则是单根
		i_count++;
	}
	return;
}                     //如果是0次常数项则无重根?
/*! @function
********************************************************************************
<PRE>
函数名   : int CalculateTransform(int i_Count,double * dp_array)
功能     : 计算数组的变号数
参数     : i_Count 数组元素个数，dp_array数组
返回值   : 返回变号数
抛出异常 : 
--------------------------------------------------------------------------------
备注     : 
典型用法 : 
--------------------------------------------------------------------------------
作者     : <xxx>
</PRE>
*******************************************************************************/
int CalculateTransform(int i_Count,double * dp_array)
{
	int i = 0;
	int count = 0;
	double temp = dp_array[0];
	for (i=1;i<i_Count;i++)
	{
		if (fabs(dp_array[i])<1e-6) //if (fabs(dp_array[i])<1e-6)?
		{
			continue;
		}
		else
		{
			if (temp*dp_array[i]<0)
			{
				count++;
			}
			temp=dp_array[i];
		}
	}
	return count;

}
/*! @function
********************************************************************************
<PRE>
函数名   : int RootCount(double d_Left,double d_Right,doule ** dp_sturm,int i_Count)
功能     : 计算d_Left 和d_Right之间的根的个数
参数     : double d_Left 下界，double d_Right 上界,double **dp_sturm sturm数组,i_Count，sturm数组的个数
返回值   : int 返回界限内的根的个数
抛出异常 : 
--------------------------------------------------------------------------------
备注     : 
典型用法 : 
--------------------------------------------------------------------------------
作者     : <xxx>
</PRE>
*******************************************************************************/
int RootCount(double d_Left,double d_Right,double ** dp_sturm,int i_Count,int i_degree)
{
	double * dp_LeftSturm = new double[i_Count];
	double * dp_RightSturm = new double[i_Count];
	int i = 0;
	for (i = 0; i< i_Count; i++)
	{
		dp_LeftSturm[i] = CalculateEquation(dp_sturm[i], i_degree, d_Left);
		dp_RightSturm[i] = CalculateEquation(dp_sturm[i], i_degree, d_Right);
	}
//	ArrayDisplay(dp_LeftSturm,i_Count);
//	ArrayDisplay(dp_RightSturm,i_Count);
	int i_LeftCount = CalculateTransform(i_Count,dp_LeftSturm);
	int i_RightCount = CalculateTransform(i_Count,dp_RightSturm);
	return i_RightCount-i_LeftCount;
}
/*! @function
********************************************************************************
<PRE>
函数名   : double * RealEquation(double* dp_equation, int i_degree, int & i_realDegree)
功能     : 返回去除高次项为零的方程
参数     : double * dp_equation 输入的方程，i_degree 输入方程次数，i_realDegree 输出最高次项不为零的方程的次数
返回值   : double * 返回最高次项不为零的方程
抛出异常 : 
--------------------------------------------------------------------------------
备注     : 
典型用法 : 
--------------------------------------------------------------------------------
作者     : <xxx>
</PRE>
*******************************************************************************/
double * RealEquation(double* dp_equation, int i_degree, int & i_realDegree)
{
	i_realDegree = GetDegree(dp_equation,i_degree);
	if (i_realDegree == i_degree)
	{
		return dp_equation;
	}
	else
	{
		double * dp_newEquation = new double[i_realDegree+1];
		int i = i_degree-i_realDegree;
		int k = 0;
		for (;i<=i_degree;i++)
		{
			dp_newEquation[k++] = dp_equation[i]; 
		}
		return dp_newEquation;
	}



}
/*! @function
********************************************************************************
<PRE>
函数名   : void BiRoot(int i_degree, double** dp_sturm,int i_sturmCount,double& d_left,double& d_right)
功能     : 
参数     : 
返回值   : 
抛出异常 : 
--------------------------------------------------------------------------------
备注     : 
典型用法 : 
--------------------------------------------------------------------------------
作者     : <xxx>
</PRE>
*******************************************************************************/
bool BiRoot(int i_degree, double** dp_sturm,int i_sturmCount,double& d_left,double& d_right,double d_precision)
{

	int i_rootCount=RootCount(d_left,d_right,dp_sturm,i_sturmCount,i_degree);
	if (i_rootCount!=0)
	{                                          /*注------------------------*/
		if (1 == i_rootCount||-1==i_rootCount) ///怎么会有-1？不懂？
		{
			return true;							
		}
		else
		{
			if (RootCount(d_left,(d_left+d_right)/2,dp_sturm,i_sturmCount,i_degree)!=0)
			{
				d_right = (d_right+d_left)/2;
				return	BiRoot(i_degree,dp_sturm,i_sturmCount,d_left,d_right,d_precision);
			}
			else
			{
				d_left = (d_left+d_right)/2;
				return	BiRoot(i_degree,dp_sturm,i_sturmCount,d_left,d_right,d_precision);
			}
		}
	}
	else
	{
		return false;								//返回状态0，无根
	}
}
/*! @function
********************************************************************************
<PRE>
函数名   : double SolveEquation(double* dp_equation, int i_degree, double** dp_sturm,int i_sturmCount)
功能     : 求解方程
参数     : dp_equation 方程式子，i_degree 方程次数， dp_sturm 方程sturm数组， i_sturmCount sturm数组个数
返回值   : 
抛出异常 : 
--------------------------------------------------------------------------------
备注     : 
典型用法 : 
--------------------------------------------------------------------------------
作者     : <xxx>
</PRE>
*******************************************************************************/
bool SolveEquation(double* dp_equation, int i_degree, double d_precision,double& d_result)
{

	
	int i = 0;

	//dp_equation,i_degree是输入的方程，dp_realEquation,i_realDegree是实际处理的方程
	//////////////////////////////////////////////////////////////////////////
	//对方程处理，使最高次项不为零
	int i_realDegree = GetDegree(dp_equation,i_degree);
	double * dp_realEquation = new double[i_realDegree];
	dp_realEquation = RealEquation(dp_equation, i_degree, i_realDegree);

	/************************************************************************/
	/*//计算sturm数组                                                       */
	/************************************************************************/
	//////////////////////////////////////////////////////////////////////////
	//构造sturm数组
	double ** dp_sturm = new double*[100];
	int i_sturmCount = 100;
	for (i = 0;i<100;i++)
	{
		dp_sturm[i] = new double[i_realDegree+1];
	}
	//////////////////////////////////////////////////////////////////////////
	//计算
	Sturm(dp_realEquation,i_realDegree,dp_sturm,i_sturmCount);

	for (i=0;i<i_sturmCount;i++)
	{
		ArrayDisplay(dp_sturm[i],i_realDegree);
	}
	int temp=i_realDegree;
	//有重根
	if (GetDegree(dp_sturm[i_sturmCount-1],i_realDegree) != 0)
	{
		delete dp_realEquation;
		dp_realEquation = NULL;
		i_realDegree = GetDegree(dp_sturm[0],i_realDegree)-GetDegree(dp_sturm[i_sturmCount-1],i_realDegree);
		
		dp_realEquation = new double[i_realDegree+1];
		dp_realEquation = NewEquation(dp_sturm[0],dp_sturm[i_sturmCount-1],temp);
		dp_realEquation = RealEquation(dp_realEquation, i_realDegree, i_realDegree);
		delete dp_sturm;
		dp_sturm = new double*[100];
		i_sturmCount = 100;
		for (i = 0;i<100;i++)
		{
			dp_sturm[i] = new double[i_realDegree+1];
		}
		Sturm(dp_realEquation,i_degree,dp_sturm,i_sturmCount);
		for (i=0;i<i_sturmCount;i++)
		{
			ArrayDisplay(dp_sturm[i],i_realDegree);
		}
	}


	//////////////////////////////////////////////////////////////////////////
	//求方程解的界限
	double	d_right=CalculateBoundary(dp_realEquation,i_degree);
	double	d_left=-d_right;

	//判断方程是否有根
	bool b_status = BiRoot(i_realDegree,dp_sturm,i_sturmCount,d_left,d_right,d_precision);

	if (b_status )
	{
		cout<<d_left<<"       "<<d_right<<endl;
		d_result = Bisection(dp_realEquation,i_realDegree,d_left,d_right,d_precision);
		cout<<"The solution of this equation is "<<d_result<<endl;
		cout<<"The value of this equation is "<<CalculateEquation(dp_realEquation,i_realDegree,d_result);
		return true;
	}

	else
	{
		cout<<"no solution";
		return false;
	}

// 	if (fabs(CalculateEquation(dp_equation,i_degree,d_right))<=d_precision)
// 	{
// 		return d_right;
// 	}
// 	
// 	if (fabs(CalculateEquation(dp_equation,i_degree,d_left))<=d_precision)
// 	{
// 		return d_left;
// 	}
// 	bool status = BiRoot(i_degree,dp_sturm,i_sturmCount,d_left,d_right,d_precision);
// 	if (status)
// 	{
// 		cout<<d_left<<"       "<<d_right<<endl;
// 		d_result = Bisection(dp_equation,i_degree,d_left,d_right,d_precision);
// 		cout<<"The solution of this equation is "<<d_result<<endl;
// 		cout<<"The value of this equation is "<<CalculateEquation(dp_equation,i_degree,d_result);
// 		return true;
// 	}
// 	else
// 	{
// 		cout<<"no solution";
// 		return false;
// 	}

}
/*! @function
********************************************************************************
<PRE>
函数名   : 
功能     : 
参数     : 
返回值   : 
抛出异常 : 
--------------------------------------------------------------------------------
备注     : 
典型用法 : 
--------------------------------------------------------------------------------
作者     : <xxx>
</PRE>
*******************************************************************************/
bool SixSolve(double d_parm0, double d_parm1, double d_parm2, double d_parm3, double d_parm4, double d_parm5, double d_parm6, double dp_precision,double& d_result)
{
	double * dp_equation = new double[7];
	dp_equation[0] = d_parm0;
	dp_equation[1] = d_parm1;
	dp_equation[2] = d_parm2;
	dp_equation[3] = d_parm3;
	dp_equation[4] = d_parm4;
	dp_equation[5] = d_parm5;
	dp_equation[6] = d_parm6;

	int i_degree = 6;
	
	bool b_haveSolution = SolveEquation(dp_equation, i_degree, dp_precision, d_result);
	cout <<d_result<<"                        "<<CalculateEquation(dp_equation,i_degree,d_result)<<endl;
	return i_degree;

}
#endif