s.c

c++解一员三次方程源码

#include<stdio.h>
#include<math.h>
#define pi 3.14159265
static double r3[3][2],r4[4][2];
void cubic(double a0,double a1,double a2,double a3)
/*求一元三次方程a0x^3+a1x^2+a2x+a3=0的根x1,x2,x3*/
{
  double b1,b2,b3,p,q,u1,v1,u2,v2,u,v,s,x1r,x1i,x2r,x2i,x3r,x3i;
  b1=a1/a0;         b2=a2/a0;       b3=a3/a0;
  p=b2/3.0-b1*b1/9.0;
  q=2*pow(b1,3)/27.0-b1*b2/3.0+b3;
  if(q*q/4.0+p*p*p>=0)
          { u1=-0.5*q+sqrt(0.25*q*q+p*p*p);
            v1=-0.5*q-sqrt(0.25*q*q+p*p*p);
            u=(2.0*u1+1.0/u1)/3.0;/*求u1的立方根u*/
            u2=u1;
            while(fabs((u2-u)/u2)>=1E-6){
                u2=u;
                u=(2.0*u2+u1/(u2*u2))/3.0;
            }
            v=(2.0*v1+1.0/v1)/3.0;/*求v1的立方根v*/
            v2=v1;
            while(fabs((v2-v)/v2)>=1E-6){
                v2=v;
                v=(2.0*v2+v1/(v2*v2))/3.0;
            }
            x1r=u+v-b1/3;               x1i=0;
            x2r=-1*(u/2.0+v/2.0+b1/3.0);x2i=(u-v)*pow(3,0.5)/2.0;
            x3r=-1*(u/2.0+v/2.0+b1/3.0);x3i=-1*(u-v)*pow(3,0.5)/2.0;
          };
  if(q*q/4.0+p*p*p<0)
          { if(q==0)s=pi/2;
            if(q<0)s=atan(-1.0*sqrt(-1*q*q-4*p*p*p)/q);
            if(q>0)s=pi+atan(sqrt(-q*q-4*p*p*p)/q);
            x1r=2*sqrt(-1.0*p)*cos(s/3.0)-b1/3;       x1i=0.0;
            x2r=2*sqrt(-1.0*p)*cos((s+2*pi)/3.0)-b1/3;x2i=0.0;
            x3r=2*sqrt(-1.0*p)*cos((s+4*pi)/3.0)-b1/3;x3i=0.0;
          };
  r3[0][0]=x1r;r3[0][1]=x1i;r3[1][0]=x2r;r3[1][1]=x2i;r3[2][0]=x3r;r3[2][1]=x3i;
}
main()
{
  float a0,a1,a2,a3;
  printf("输入参数:");
  scanf("%f%f%f%f",&a0,&a1,&a2,&a3);
  cubic(a0,a1,a2,a3);
  printf("x1的实部为:%f\t",r3[0][0]);
  printf("x1的虚部为:%f\n",r3[0][1]);
  printf("x2的实部为:%f\t",r3[1][0]);
  printf("x2的虚部为:%f\n",r3[1][1]);
  printf("x3的实部为:%f\t",r3[2][0]);
  printf("x3的虚部为:%f\n",r3[2][1]);
 
}