/********************************************* * Solving algebraic equations * * by Bairstow's Method * * ------------------------------------------ * * Ref.: "Mathématiques par l'informatique * * individuelle - Programmes en basic * * By H. Lehning and D. Jakubowicz, * * MASSON Paris, 1982" [BIBLI 06]. * * * * SAMPLE RUN: * * (Find real and complex roots of equation: * * x^6 -127x^5 +215x^4 +28x^3 -39x^2 +20x * * -15 = 0) * * * * Input degree of equation: 6 * * A(6) = 1 * * A(5) = -127 * * A(4) = 215 * * A(3) = 28 * * A(2) = -39 * * A(1) = 20 * * A(0) = -15 * * * * ROOTS ERROR * * 0.03989619 0.44667179 1.874776e-011 * * 0.03989619 -0.44667179 1.874776e-011 * * 0.52383379 0.00000000 1.297699e-006 * * -0.64575255 0.00000000 3.497111e-006 * * 125.28210889 0.00000000 3.889075e-009 * * 1.76001412 0.00000000 1.438089e-006 * * * * NOTE: if a better precision is desired, * * you can set e to a smaller value or * * use Newton's method for a real root. * * * * C++ version by J-P Moreau. * * (www.jpmoreau.fr) * ********************************************** See explanation file Bairstow.txt --------------------------------- */ #include #include //Labels: l100,l200,l500,l550,l1000,fin double A[20],B[20],C[20],P[20]; double a1,b1,d,e,f,p1,q,r,t,u,v,x,y,z; int i,j,k,m,n; void Print_results() { //Label: l10 double a1; int i; printf("%13.8f %13.8f", x,y); //calculate error estimation (optional) u=P[0]; v=0.0; for (i=1; i<=m; i++) { r=u*x-v*y; v=u*y+v*x; u=r+P[i]; } a1=sqrt(u*u+v*v); u=0.0; v=0.0; for (i=1; i<=m; i++) { r=u*x-v*y; v=u*y+v*x; u=r+(m-i+1)*P[i-1]; } if (a1==0.0) goto l10; a1=a1/sqrt(u*u+v*v); l10: printf(" %e\n", a1); } //Solve x^2+p1x+q=0 void Solve2() { //Label: l10 double d; d=p1*p1-4*q; if (d<0) goto l10; d=sqrt(d); y=0.0; x=(-p1+d)/2.0; Print_results(); x=(-p1-d)/2.0; Print_results(); return; l10: d=sqrt(-d)/2.0; x=-p1/2.0; y=d; Print_results(); y=-d; Print_results(); } void main() { //Enter polynomial A(i) and copy in P(i) printf("\n Input degree of equation: "); scanf("%d",&n); m=n; //save N for error estimation for (i=0; i<=n; i++) { printf(" A(%d) = ", n-i); scanf("%lf",&A[i]); P[i]=A[i]; } printf("\n ROOTS ERROR\n"); //Init section p1=0.0; q=0.0; k=100; e=0.001; //Factoring main loop l100:if (n<=2) goto l500; j=0; l200: if (j>k) goto l1000; j++; //calculate B(i) and C(i) B[0]=0.0; B[1]=0.0; C[0]=0.0; C[1]=0.0; for (i=2; i<=n+2; i++) { B[i]=A[i-2]-p1*B[i-1]-q*B[i-2]; C[i]=-B[i-1]-p1*C[i-1]-q*C[i-2]; } //calculate dp=a1 and dq=b1 x=B[n+1]; y=B[n+2]; z=C[n]; t=C[n+1]; u=C[n+2]; d=t*t-z*(u+x); if (d==0) goto l1000; a1=(z*y-x*t)/d; b1=(-x*(q*z+p1*t)-y*t)/d; //New p1 and q p1=p1+a1; q=q+b1; f=(fabs(a1)+fabs(b1))/(fabs(p1)+fabs(q)); if (f>e) goto l200; //A factor has been found Solve2(); //Update polynomial n=n-2; for (i=0; i<=n; i++) A[i]=B[i+2]; goto l100; l500: //Last factor, first or second degree if (n==2) goto l550; x=-A[1]/A[0]; y=0.0; //first degree Print_results(); printf("\n"); return; l550: p1=A[1]/A[0]; q=A[2]/A[0]; Solve2(); printf("\n"); return; l1000: printf("\n Process not convergent !\n\n"); } //end of file bairsto1.cpp