一元三次和四次方程的根式解法
http://www.qthome.org/bbs/pascal_view.asp?id=28一元三次和四次方程的根式解法
基本信息：
提交人：tianyuan_0 提交日期：2006-1-17 23:44:40 点击数：2075
编辑 删除
程序简介：
前面看到了有人用二分法解一元三次方程.我正在编的程序也需要求解一元三次方程.
代数学早已解决了一元三次和一元四次方程的求解问题,因而再用二分法求解无疑是效率不高的.
下面的程序省略了主函数.
平台是Visual C++ 2005.
由于我没有在VC++2005的<complex>中找到操作符(这也是令人奇怪的一件事),我重写了complex类.
程序内容：
namespace mathapp{
const double PAI=3.141592653589793626;
public ref struct mathappEx:public System::Exception
{
public:
int i;
System::String^ Desp;
mathappEx(int n,System::String^ s)
{
i=n;
Desp=s;
}
};
///complex
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
value class complex
{
public:double _im,_re;
public:
complex(double re,double im)
{
_im=im;
_re=re;
}
//
//
//property
//
//
property double mo
{
double get()
{
return System::Math::Sqrt(_im*_im+_re*_re);
}
}
property double arg
{
double get()
{
if(_re==0)
{
if(_im>0)
return PAI/2;
else
return -PAI/2;
}
else
{
if(_re>0)
return System::Math::Atan(_im/_re);
else
return PAI+System::Math::Atan(_im/_re);
}
}
}
//
//operator
//
static complex operator ^(complex c,double up)
{
double t=System::Math::Pow(c.mo,up);
return complex::complex(t*System::Math::Cos(c.arg*up),t*System::Math::Sin(c.arg*up));
}
static complex operator +(complex left,complex right)
{
return complex::complex(left._re+right._re,left._im+right._im);
}
static complex operator +(complex left,double d)
{
return complex::complex(left._re+d,left._im);
}
static complex operator +(double d,complex left)
{
return complex::complex(d+left._re,left._im);
}
static complex operator -(complex left,complex right)
{
return complex::complex(left._re-right._re,left._im-right._im);
}
static complex operator-(double left,complex right)
{
return complex::complex(left-right._re,-right._im);
}
static complex operator-(complex right,double left)
{
return complex::complex(right._re-left,right._im);
}
static complex operator *(complex left,complex right)
{
return complex::complex(left._re*right._re-left._im*right._im,left._re*right._im+left._im*right._re);
}
static complex operator *(double left,complex right)
{
return complex::complex(left*right._re,left*right._im);
}
static complex operator *(complex left,double right)
{
return complex::complex(left._re*right,left._im*right);
}
static complex operator /(complex left,complex right)
{
double t=right._im*right._im+right._re*right._re;
if (t==0)
{
mathappEx^ e=gcnew mathappEx(0,"除数为0");
throw e;
}
else
{
complex c=complex::complex(right._re,-right._im);
return t*c*left;
}
}
static complex operator /(complex left,double right)
{
if(right==0)
{
mathappEx^e=gcnew mathappEx(0,"除数为0");
throw e;
}
else
{
return complex::complex(left._re/right,left._im/right);
}
}
static complex operator /(double left,complex right)
{
double t=right._im*right._im+right._re*right._re;
if (t==0)
{
mathappEx^ e=gcnew mathappEx(0,"除数为0");
throw e;
}
else
{
complex c=complex::complex(right._re,-right._im);
return c*left/t;
}
}
static intoperator==(complex__z1, complex__z2)
{
return __z1._re == __z2._re && __z1._im == __z2._im;
}
static int operator==(complex c,double d)
{
return (c._re==d)&&(c._im==0.0);
}
static int operator==(double d,complex c)
{
return (c._re==d)&&(c._im==0.0);
}
static int operator!=(complex__z1, complex__z2)
{
return __z1._re != __z2._re || __z1._im != __z2._im;
}
static int operator!=(complexc,doubled)
{
return c._im!=0||c._re!=d;
}
static int operator!=(double d,complex c)
{
return c!=d;
}
complex operator +()
{
return *this;
}
complex operator -()
{
return complex::complex(-_re,-_im);
}
//
//convert
//
virtual System::String^ ToString()override
{
if(_im>=0)
return this->_re.ToString()+"+"+this->_im.ToString()+"i";
else
return this->_re.ToString()+"-"+this->_im.ToString()+"i";
}
};
const complex omiga(-0.5,0.86602540378443864676372317075294);
const complex omiga2(-0.5,-0.86602540378443864676372317075294);
static int Sign(complex c)
{
if(c._re>0)
return 1;
if(c._re=0)
return 0;
return -1;
}
////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////
ref class MathApp
{
public:
MathApp(void);
static int OneDegree(double x,double c,double &jie)
{
if(x==0) return 0;
jie=-c/x;
return 1;
}
static int QuadricOne(double xx,double x,double c,cli::array<complex>^ root)
{
if(xx==0)
{
double t;
if(OneDegree(x,c,t)!=0)
{
root._re=t;
root._im=0;
return 1;
}
else
return 0;
}
complex cc(x*x-4*xx*c,0);
root=(-x+(cc^0.5))/(2*xx);
root=(-x-(cc^0.5))/(2*xx);
if(cc==0.0) return 1;
return 2;
}
static int CubicOne(double xxx,double xx,double x,double c,cli::array<complex, 1>^ root)//一元三次方程求解
{
/*1.x^3+px+q=0的解：
xk=w^k*(-q/2+((q/2)^2+(p/3)^3))^(1/2))^(1/3)+
w^(3-k)*(-q/2-((q/2)^2+(p/3)^3))^(1/2))^(1/3)
其中：w=cos(2*pi/3)+i*sin(2*pi/3),k=1,2,3
2.a*x^3+b*x^2+c*x+d=0
令：y=x-b/(3*a)化为1.*/
if(xxx==0)
return QuadricOne(xx,x,c,root);
xx=xx/xxx;
x=x/xxx;
c=c/xxx;
double p=x*0.33333333333333333333333-xx*xx*0.1111111111111111111111111111111;
double q=xx*xx*xx/27.0 -xx*x*0.166666666666666666666666666 +c*0.5;
complex d(p*p*p+q*q,0.0);
complex d1=d^0.5;
d=(d1-q)^(0.33333333333333333333333333333);
if(d.mo!=0)
d1=-p/d;
else
d1=(-q-d1)^(0.333333333333333333333333333);
root=d+d1-xx*0.3333333333333333333333333333;
root=d*omiga+d1*omiga2-xx*0.3333333333333333333333333333;
root=d*omiga2+d1*omiga-xx*0.3333333333333333333333333333;
return 3;
}
static int Quartic(double xxxx,double xxx,double xx,double x,double c,cli::array<complex>^ root)//一元四次方程求解
{
if(xxxx==0)
return CubicOne(xxx,xx,x,c,root);
xxx=xxx/xxxx;
xx=xx/xxxx;
x=x/xxxx;
c=c/xxxx;
double p=-0.375000000*xxx*xxx+xx;
double q=0.12500000*xxx*xxx*xxx-0.50*xxx*xx+x;
double r=-0.01171875*xxx*xxx*xxx*xxx+0.0625*xxx*xxx*xx-0.25*xxx*x+c;
cli::array<complex>^ carray=gcnew array<complex>(3);
int j=CubicOne(1.0,-2*p,p*p-4*r,q*q,carray);
for(j=0;j<3;j++)
carray=(-carray)^0.5;
if(Sign(carray*carray*carray)*System::Math::Sign(q)<0)
carray=-carray;
xxx=0.250*xxx;
root=0.50*(carray-carray+carray)-xxx;
root=0.50*(carray+carray-carray)-xxx;
root=0.50*(-carray+carray+carray)-xxx;
root=0.50*(-carray-carray-carray)-xxx;
return 1;
}
};
}
//yaos注： 未验证正确性 其实用牛顿迭代法计算（3次和4次方程）的数值解效率已经很高了。 呵呵
实用的不优美
优美的不实用啊
对了mathe
知道5次以上代数方程的三角函数解法么？ 5次以上代数方程的三角函数解法
:Q:
google没找到，可否详之 不知道5次方程的三角函数解法。
不过我知道五次方程的解可以用椭圆函数来表示，这个可以google到的 5次以下的根式解法有现成的公式，
5次以上的用群论的知识根据所给的系数也是可以解出的
页:
[1]