/*********************************************
*        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 <stdio.h>
#include <math.h>

//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