header.h

求解一元高次方程

/*! @file
********************************************************************************
<PRE>
模块名       : 函数模块
文件名       : header.h
相关文件     : 
文件实现功能 : 包含了一系列函数，主要有输出函数，计算界的函数，二分法函数，strum函数
，等
作者         : <xxx>
版本         : 1.0
--------------------------------------------------------------------------------
备注         : 
--------------------------------------------------------------------------------
修改记录 : 
日 期						版本				修改人					修改内容 
YYYY/MM/DD					1.0					<xxx>					  创建
2008/12/19					1.0					<Name>					  创建
 </PRE>
********************************************************************************

* 版权所有(c) YYYY, <xxx>, 保留所有权利

*******************************************************************************/

#ifndef HEADER_INCLUDED
#define HEADER_INCLUDED
/*! @function
********************************************************************************
<PRE>
函数名   : ArrayDisplay
功能     : 输出一个double数组
参数     : double*dp_array,int i_degree
返回值   : 
抛出异常 : 
--------------------------------------------------------------------------------
备注     : 
典型用法 : ArrayDisplay(arrayname,degree_of_array);
--------------------------------------------------------------------------------
作者     : <xxx>
</PRE>
*******************************************************************************/
void ArrayDisplay(double*dp_array,int i_degree)
{
	cout<<endl<<"---------------------------------------------------------------------"<<endl;
	cout<<"数组元素分别为:"<<endl;
	for (int i=0;i<=i_degree;i++)
	{
		cout<<dp_array[i]<<"      ";
	}
	cout<<endl;
	cout<<"---------------------------------------------------------------------"<<endl<<endl;

}

/*! @function
********************************************************************************
<PRE>
函数名   : CalculateEquation
功能     : 计算一个函数值
参数     : double*dp_array,int i_degree,double d_value
返回值   : double 
抛出异常 : 
--------------------------------------------------------------------------------
备注     : 
典型用法 : double test=CalculateEquation(dp_array,i_degree,d_value);
--------------------------------------------------------------------------------
作者     : <xxx>
</PRE>
*******************************************************************************/
double CalculateEquation(double*dp_array,int i_degree,double d_value)
{
	double result=0;
	for (int i=0;i<=i_degree;i++)
	{
		double temp=pow(d_value,i_degree-i);
		result+=dp_array[i]*temp;
	}
	return result;
}
/*! @function
********************************************************************************
<PRE>
函数名   : int GetDegree(double * dp_equation,int i_degree)
功能     : 返回一个方程的次数
参数     : double * dp_equation 方程式，i_degree 方程的原次数
返回值   : int 返回方程的实际次数
抛出异常 : 
--------------------------------------------------------------------------------
备注     : 
典型用法 : 
--------------------------------------------------------------------------------
作者     : <xxx>
</PRE>
*******************************************************************************/
int GetDegree(double * dp_equation,int i_degree)
{
	int i = 0;
	int i_realDegree = i_degree;
	for (i = 0;i<= i_degree;i++)
	{
		if (fabs(dp_equation[i])<= 1e-6)
		{
			i_realDegree--;
		}
		else
		{
			break;
		}
	}
	return i_realDegree;

}
/*! @function
********************************************************************************
<PRE>
函数名   : double CalculateBoundary(double*dpArray,int i_degree)
功能     : 求一元高次方程解的界限
参数     : double*dpArray,int i_degree
返回值   : 解的上界值
抛出异常 : 
--------------------------------------------------------------------------------
备注     : 
典型用法 : CalculateBoundary(dpArray,i_degree);
--------------------------------------------------------------------------------
作者     : <xxx>
</PRE>
*******************************************************************************/
double CalculateBoundary(double*dpArray,int i_degree)
{
	double a0=fabs(dpArray[0]);
	double dTemp=0.0;
	double dAbsoluteValue=0.0;      /*注------------------------*/
	for (int i=1;i<=i_degree;i++)  //是否应该是i=1？根的界限中dTemp好像不包括a0？
	{
		dAbsoluteValue=fabs(dpArray[i]);
		if (dAbsoluteValue>dTemp)
		{
			dTemp=dAbsoluteValue;
		}
	}
	return 1.0+(dTemp/a0);
}
/*! @function
********************************************************************************
<PRE>
函数名   : double Bisection(double*dp_array,int i_degree,double d_left,
			double d_right,double d_precision)
功能     : 二分法求解方程的一个解
参数     : double*dp_array,int i_degree,double d_left,double d_right,double d_precision
返回值   : 返回方程的一个解
抛出异常 : 
--------------------------------------------------------------------------------
备注     : 用d_precision控制精度
典型用法 : cout<<Bisection(dp_coefficient,i_degree,0.0,3.0,0.0005);
--------------------------------------------------------------------------------
作者     : <xxx>
</PRE>
*******************************************************************************/
double Bisection(double*dp_array,int i_degree,double d_left,double d_right,double d_precision)
{
	double d_centerValue=CalculateEquation(dp_array,i_degree,(d_right+d_left)/2);

	//若在精度范围内返回
	if (fabs(d_centerValue)<=d_precision)
	{
		return (d_left+d_right)/2;
	}
	//精度不够递归求解
	else
	{
		double d_leftValue=CalculateEquation(dp_array,i_degree,d_left);
		double d_rightValue=CalculateEquation(dp_array,i_degree,d_right);
		if (d_leftValue*d_centerValue<0)
		{
			return	Bisection(dp_array,i_degree,d_left,(d_left+d_right)/2,d_precision);
		}                                                                              
		else                                                                            
		{
			return Bisection(dp_array,i_degree,(d_left+d_right)/2,d_right,d_precision);
		}	
	}
}
/*! @function
********************************************************************************
<PRE>
函数名   : double* Divide(const double* dp_Upper,const double* dp_Lower , int i_degree)
功能     : 输入dp_Upper为被除式,dp_Lower为除式，i_degree为输入的方程次数
参数     : 
返回值   : double* 返回余式
抛出异常 : 
--------------------------------------------------------------------------------
备注     : 
典型用法 : 
--------------------------------------------------------------------------------
作者     : <xxx>
</PRE>
*******************************************************************************/
double* Divide(const double* dp_UpperOriginal,const double* dp_Lower , int i_degree)
{
	int i=0;
	int i_UpperDegree=i_degree;
	int i_LowerDegree=i_degree;

	double*dp_Upper=new double[i_degree+1];
	for (i=0;i<=i_degree;i++)
	{
		dp_Upper[i]=dp_UpperOriginal[i];
	}
	//获取输出两个式子的次数
	for (i=0;i<=i_degree;i++)
	{
		if (fabs(dp_Upper[i])<=1e-6)  // if(fabs(dp_Upper[i])<1e-6)?有必要没？
		{                        
			i_UpperDegree--;     
		}                        
		else                      
		{
			break;
		}
	}
	for (i=0;i<=i_degree;i++)
	{
		if (fabs(dp_Lower[i])<=1e-6)  // if(fabs(dp_Lower[i])<1e-6)?有必要没？
		{
			i_LowerDegree--;
		}
		else
		{
			break;
		}
	}
	int j=0;
	double temp=0.0;
	while (true)
	{
		j=i_degree-i_UpperDegree;
		temp=dp_Upper[i_degree-i_UpperDegree];
		for (i=i_degree-i_LowerDegree;i<=i_degree;i++)
		{
			dp_Upper[j]=dp_Upper[j]-(temp*dp_Lower[i])/dp_Lower[i_degree-i_LowerDegree];
			j++;
		}
		i_UpperDegree=i_degree;
		for (i=0;i<=i_degree;i++)
		{
			if (fabs(dp_Upper[i])<=1e-6)// if(fabs(dp_Upper[i])<1e-6)?dp_Upper[0]必是0.0，但dp_Upper[1]不一定0.0，这时近似是0？
			{                          //如：// 3 3.000001   
				i_UpperDegree--;            //- 3 3
			}                               //-------------
			else                             // 0  0.000001
			{
				break;
			}
		}
		if (i_UpperDegree<i_LowerDegree)
		{
			break;
		}
	}

	for (i=0;i<=i_degree;i++)
	{
		dp_Upper[i]=-dp_Upper[i];
	}
	return dp_Upper;
}

/*! @function
********************************************************************************
<PRE>
函数名   : 
功能     : 
参数     : 
返回值   : 
抛出异常 : 
--------------------------------------------------------------------------------
备注     : 
典型用法 : 
--------------------------------------------------------------------------------
作者     : <xxx>
</PRE>
*******************************************************************************/
double* NewEquation(const double* dp_UpperOriginal,double* dp_Lower , int i_degree)
{
	int i=0;
	int i_UpperDegree=i_degree;
	int i_LowerDegree=i_degree;
	
	double*dp_Upper=new double[i_degree+1];
	for (i=0;i<=i_degree;i++)
	{
		dp_Upper[i]=dp_UpperOriginal[i];
	}
	//获取输出两个式子的次数
	i_UpperDegree = GetDegree(dp_Upper,i_degree);
// 	for (i=0;i<=i_degree;i++)
// 	{
// 		if (fabs(dp_Upper[i])<=1e-6)  // if(fabs(dp_Upper[i])<1e-6)?有必要没？
// 		{                        
// 			i_UpperDegree--;     
// 		}                        
// 		else                      
// 		{
// 			break;
// 		}
// 	}
	i_LowerDegree = GetDegree(dp_Lower, i_degree);
// 	for (i=0;i<=i_degree;i++)
// 	{
// 		if (fabs(dp_Lower[i])<=1e-6)  // if(fabs(dp_Lower[i])<1e-6)?有必要没？
// 		{
// 			i_LowerDegree--;
// 		}
// 		else
// 		{
// 			break;
// 		}
// 	}
	int j=0;
	double temp=0.0;
	double * dp_equation;
	int i_newDegree = i_UpperDegree-i_LowerDegree;
	dp_equation = new double [i_newDegree+1];
	for (int k=0;k<=i_newDegree;k++)
	{
		j=i_degree-i_UpperDegree;
		temp=dp_Upper[i_degree-i_UpperDegree];
		for (i=i_degree-i_LowerDegree;i<=i_degree;i++)
		{
			dp_Upper[j]=dp_Upper[j]-(temp*dp_Lower[i])/dp_Lower[i_degree-i_LowerDegree];
			j++;
 		}
		i_UpperDegree--;
// 		i_UpperDegree=i_degree;
// 		for (i=0;i<=i_degree;i++)
// 		{
// 			if (fabs(dp_Upper[i])<=1e-6)// if(fabs(dp_Upper[i])<1e-6)?dp_Upper[0]必是0.0，但dp_Upper[1]不一定0.0，这时近似是0？
// 			{                          //如：// 3 3.000001   
// 				i_UpperDegree--;            //- 3 3
// 			}                               //-------------
// 			else                             // 0  0.000001
// 			{
// 				break;
// 			}
//  		}
		dp_equation[k] = temp/dp_Lower[i_degree-i_LowerDegree];
	}                                                          
	i_degree = i_newDegree;
	return dp_equation;
}

/*! @function
********************************************************************************
<PRE>
函数名   : void Sturm(double* dp_equation,int i_degree,double** &dp_sturm)
功能     : 计算sturm数组
参数     : double* dp_equation [in],int i_degree [in],double** &dp_sturm [out]
返回值   : void
抛出异常 : 
--------------------------------------------------------------------------------
备注     : 
典型用法 : Sturm(dp_equation,i_realDegree,dp_sturm);