quartic.txt

求解一元三次议程的源代码,用于求解一元三次方程

#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,x2,x3,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()/*求一元四次方程a0x^4+a1x^3+a2x^2+a3x+a4=0的根x1,x2,x3,x4*/
{
  double a0=2.0,a1=3.0,a2=-5.0,a3=-4.0,a4=4.0,b1,b2,b3,p1,p2,p3,p4,x_21r,x_22r,r_2,s_1,s_2,d1,d2,d3;
  p1=a1/a0;p2=a2/a0;p3=a3/a0;p4=a4/a0;
  b1=-3*p1*p1/8.0+p2;
  b2=p1*p1*p1/8.0-p1*p2/2.0+p3;
  b3=-3*p1*p1*p1*p1/256.0+p1*p1*p2/16.0-p1*p3/4.0+p4;
  if(b2==0)
       {if(b1*b1-4.0*b3>=0)
                 { x_21r=-0.5*b1+sqrt(0.25*b1*b1-b3);
                   x_22r=-0.5*b1-sqrt(0.25*b1*b1-b3);
                   if(x_21r>=0)
                          {r4[0][0]=sqrt(x_21r)-0.25*p1;      r4[0][1]=0.0;
                           r4[1][0]=-1*sqrt(x_21r)-0.25*p1;   r4[1][1]=0.0;
                          };
                   if(x_21r<0)
                          {r4[0][0]=-0.25*p1;      r4[0][1]=sqrt(-1*x_21r);
                           r4[1][0]=-0.25*p1;      r4[1][1]=-1*sqrt(-1*x_21r);
                          };
                   if(x_22r>=0)
                          {r4[2][0]=sqrt(x_22r)-0.25*p1;      r4[2][1]=0.0;
                           r4[3][0]=-1*sqrt(x_22r)-0.25*p1;   r4[3][1]=0.0;
                          };
                   if(x_22r<0)
                          {r4[2][0]=-0.25*p1;      r4[2][1]=sqrt(-1*x_22r);
                           r4[3][0]=-0.25*p1;      r4[3][1]=-1*sqrt(-1*x_22r);
                          };
                 };
         if(b1*b1-4.0*b3<0)
                 { r_2=sqrt(b3);
                   s_1=atan(-2.0*sqrt(b3-0.25*b1*b1)/b1);
                   s_2=atan(2.0*sqrt(b3-0.25*b1*b1)/b1);
                   r4[0][0]=sqrt(r_2)*cos(s_1/2.0)-0.25*p1;         r4[0][1]=sqrt(r_2)*sin(s_1/2.0);
                   r4[1][0]=sqrt(r_2)*cos(s_1/2.0+pi)-0.25*p1;      r4[1][1]=sqrt(r_2)*sin(s_1/2.0+pi);
                   r4[2][0]=sqrt(r_2)*cos(s_2/2.0)-0.25*p1;         r4[2][1]=sqrt(r_2)*sin(s_2/2.0);
                   r4[3][0]=sqrt(r_2)*cos(s_2/2.0+pi)-0.25*p1;      r4[3][1]=sqrt(r_2)*sin(s_2/2.0+pi); 
                 };
        };
  if(b2!=0)
       {cubic(1.0,b1*2.0,(b1*b1-4.0*b3),-1*b2*b2);
        if(r3[0][0]!=0&&r3[0][1]==0)d1=sqrt(r3[0][0]);
        if(r3[1][0]!=0&&r3[1][1]==0)d1=sqrt(r3[1][0]);
        if(r3[2][0]!=0&&r3[2][1]==0)d1=sqrt(r3[2][0]);
        d2=(b1+d1*d1-b2/d1)/2.0;
        d3=(b1+d1*d1+b2/d1)/2.0;
        if(d1*d1-4.0*d2>=0)
                 {r4[0][0]=-0.5*d1+sqrt(0.25*d1*d1-d2)-0.25*p1;         r4[0][1]=0;
                  r4[1][0]=-0.5*d1-sqrt(0.25*d1*d1-d2)-0.25*p1;         r4[1][1]=0;
                 };
        if(d1*d1-4.0*d2<0)
                 {r4[0][0]=-0.5*d1-0.25*p1;         r4[0][1]=sqrt(d2-0.25*d1*d1);
                  r4[1][0]=-0.5*d1-0.25*p1;         r4[1][1]=-1*sqrt(d2-0.25*d1*d1);
                 };
        if(d1*d1-4.0*d3>=0)
                 {r4[2][0]=0.5*d1+sqrt(0.25*d1*d1-d3)-0.25*p1;         r4[2][1]=0;
                  r4[3][0]=0.5*d1-sqrt(0.25*d1*d1-d3)-0.25*p1;         r4[3][1]=0;
                 };
        if(d1*d1-4.0*d3<0)
                 {r4[2][0]=0.5*d1-0.25*p1;         r4[2][1]=sqrt(d3-0.25*d1*d1);
                  r4[3][0]=0.5*d1-0.25*p1;         r4[3][1]=-1*sqrt(d3-0.25*d1*d1);
                 };    
       };
  printf("x1的实部为:%f\t",r4[0][0]);
  printf("x1的虚部为:%f\n",r4[0][1]);
  printf("x2的实部为:%f\t",r4[1][0]);
  printf("x2的虚部为:%f\n",r4[1][1]);
  printf("x3的实部为:%f\t",r4[2][0]);
  printf("x3的虚部为:%f\n",r4[2][1]);
  printf("x4的实部为:%f\t",r4[3][0]);
  printf("x4的虚部为:%f\n",r4[3][1]);
}