复件 quartic1.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)

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]);
}