一元三次和四次方程的根式解法

无心人

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:ow(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 int  operator==(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!=(complex  c,double  d)
{
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[0]._re=t;
root[0]._im=0;
return 1;
}
else
return 0;
}
complex cc(x*x-4*xx*c,0);
root[0]=(-x+(cc^0.5))/(2*xx);
root[1]=(-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[0]=d+d1-xx*0.3333333333333333333333333333;
root[1]=d*omiga+d1*omiga2-xx*0.3333333333333333333333333333;
root[2]=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[j]=(-carray[j])^0.5;

if(Sign(carray[0]*carray[1]*carray[2])*System::Math::Sign(q)<0)
carray[0]=-carray[0];
            
xxx=0.250*xxx;
root[0]=0.50*(carray[0]-carray[1]+carray[2])-xxx;
root[1]=0.50*(carray[0]+carray[1]-carray[2])-xxx;
            root[2]=0.50*(-carray[0]+carray[1]+carray[2])-xxx;
root[3]=0.50*(-carray[0]-carray[1]-carray[2])-xxx;
return 1;
}
};
}


//yaos注： 未验证正确性

mathe

其实用牛顿迭代法计算（3次和4次方程）的数值解效率已经很高了。

无心人

呵呵
实用的不优美
优美的不实用啊

对了mathe
知道5次以上代数方程的三角函数解法么？

wayne

5次以上代数方程的三角函数解法

google没找到，可否详之

mathe

不知道5次方程的三角函数解法。
不过我知道五次方程的解可以用椭圆函数来表示，这个可以google到的

plp626

5次以下的根式解法有现成的公式，
5次以上的用群论的知识根据所给的系数也是可以解出的